-
Notifications
You must be signed in to change notification settings - Fork 2
/
mohr.nim
142 lines (117 loc) · 3.93 KB
/
mohr.nim
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
# From: https://csmbrannon.net/2011/04/25/stress-state-analysis-python-script/
import std/[algorithm, os, math, strutils, strformat], manu
# ------------------------------------------------------------
# Helper functions for uniform printing throughout the script.
# ------------------------------------------------------------
proc headerprint(s: string) =
# Prints a centered string to divide output sections.
const width = 64
const ch = '='
let numspaces = width - len(s)
let before = int(ceil(numspaces / 2))
let after = int(floor(numspaces / 2))
echo("\n", ch.repeat(before), s, ch.repeat(after), "\n")
proc valprint(s: string, value: float) =
# Ensure uniform formatting of scalar value outputs.
echo(s.align(30), ": ", value.formatEng)
proc matprint(s: string, value: Matrix) =
# Ensure uniform formatting of matrix value outputs.
echo(s, ":\n", value)
proc usage() =
# When the user needs help, print the script usage.
headerprint(" Analyze Stress State ")
let appname = getAppFilename().extractFilename()
let s = &""" For a given stress state this script computes many
useful quantities that help to analyze the stress state.
Currently, the following values are output:
Isotropic Matrix
Deviatoric Matrix
Principal Stresses
Maximum Shear
Mean Stress
Equivalent Stress
Invariant I1
Invariant J2
Invariant J3
Lode Coordinates
Triaxiality
Command line syntax option 1:
> ./{appname} sig11 sig22 sig33
Command line syntax option 2:
> ./{appname} sig11 sig22 sig33 sig12 sig13 sig23"""
quit(s)
# ------------------------------
# The main section of the script
# ------------------------------
proc main() =
# If the script is run with an incorrect number of arguments
# or if the user is asking for help, print the usage information.
let params = commandLineParams()
if "--help" in params or "-h" in params or
params.len != 3 and params.len != 6:
usage()
# load stress components from the command line in a temporary
# container
var dum = newSeq[float](6)
for idx in 0 ..< params.len:
try:
dum[idx] = parseFloat(params[idx])
except ValueError:
quit(&"Argument '{params[idx]}' is not a valid float")
# load the stresses into our matrix and compute the
# deviatoric and isotropic stress matricies
let
sigma = matrix(@[
@[dum[0], dum[3], dum[4]],
@[dum[3], dum[1], dum[5]],
@[dum[4], dum[5], dum[2]]
])
sigmaIso = 1.0/3.0*trace(sigma)*eye64(sigma.m)
sigmaDev = sigma - sigmaIso
# compute principal stresses
var eigVals = eig(sigma).getRealEigenvalues
sort(eigVals)
reverse(eigVals)
# compute max shear stress
let maxShear = (eigVals[0]-eigVals[2])/2.0
# compute the stress invariants
let
I1 = trace(sigma)
J2 = 1.0/2.0*trace(sigmaDev*sigmaDev)
J3 = 1.0/3.0*trace(sigmaDev*sigmaDev*sigmaDev)
# compute other common stress measures
let
meanStress = 1.0/3.0*I1
eqvStress = sqrt(3.0*J2)
# compute lode coordinates
let
lodeR = sqrt(2.0*J2)
lodeZ = I1/sqrt(3.0)
temp = 3.0*sqrt(6.0)*det(sigmaDev/lodeR)
lodeTheta = 1.0/3.0*arcsin(temp)
# compute the stress triaxiality
let triaxiality = meanStress/eqvStress
# Print out what we've found
headerprint(" Stress State Analysis ")
matprint("Input Stress", sigma)
headerprint(" Component Matricies ")
matprint("Isotropic Stress", sigmaIso)
matprint("Deviatoric Stress", sigmaDev)
headerprint(" Scalar Values ")
valprint("P1", eigVals[0])
valprint("P2", eigVals[1])
valprint("P3", eigVals[2])
valprint("Max Shear", maxShear)
valprint("Mean Stress", meanStress)
valprint("Equivalent Stress", eqvStress)
valprint("I1", I1)
valprint("J2", J2)
valprint("J3", J3)
valprint("Lode z", lodeZ)
valprint("Lode r", lodeR)
valprint("Lode theta (rad)", lodeTheta)
valprint("Lode theta (deg)", radToDeg(lodeTheta))
valprint("Triaxiality", triaxiality)
headerprint(" End Output ")
main()
# End of script