Note
Go to the end to download the full example code.
ComputeInvariants#
Compute invariants used in hyperelastic constitutive laws.
⎡cxx cxy cxz⎤
⎢ ⎥
⎢cxy cyy cyz⎥
⎢ ⎥
⎣cxz cyz czz⎦
⎡ cxx ⎤
⎢ ⎥
⎢ cyy ⎥
⎢ ⎥
⎢ czz ⎥
⎢ ⎥
⎢√2⋅cyz⎥
⎢ ⎥
⎢√2⋅cxz⎥
⎢ ⎥
⎣√2⋅cxy⎦
===================== I1 =====================
cxx + cyy + czz
dI1dC
⎡1⎤
⎢ ⎥
⎢1⎥
⎢ ⎥
⎢1⎥
⎢ ⎥
⎢0⎥
⎢ ⎥
⎢0⎥
⎢ ⎥
⎣0⎦
d2I1dC
⎡0 0 0 0 0 0⎤
⎢ ⎥
⎢0 0 0 0 0 0⎥
⎢ ⎥
⎢0 0 0 0 0 0⎥
⎢ ⎥
⎢0 0 0 0 0 0⎥
⎢ ⎥
⎢0 0 0 0 0 0⎥
⎢ ⎥
⎣0 0 0 0 0 0⎦
===================== I2 =====================
2 2 2
1.0⋅cxx⋅cyy + 1.0⋅cxx⋅czz - 1.0⋅cxy - 1.0⋅cxz + 1.0⋅cyy⋅czz - 1.0⋅cyz
dI2dC
⎡1.0⋅cyy + 1.0⋅czz⎤
⎢ ⎥
⎢1.0⋅cxx + 1.0⋅czz⎥
⎢ ⎥
⎢1.0⋅cxx + 1.0⋅cyy⎥
⎢ ⎥
⎢ -1.0⋅√2⋅cyz ⎥
⎢ ⎥
⎢ -1.0⋅√2⋅cxz ⎥
⎢ ⎥
⎣ -1.0⋅√2⋅cxy ⎦
d2I2dC
⎡ 0 1.0 1.0 0 0 0 ⎤
⎢ ⎥
⎢1.0 0 1.0 0 0 0 ⎥
⎢ ⎥
⎢1.0 1.0 0 0 0 0 ⎥
⎢ ⎥
⎢ 0 0 0 -1.0 0 0 ⎥
⎢ ⎥
⎢ 0 0 0 0 -1.0 0 ⎥
⎢ ⎥
⎣ 0 0 0 0 0 -1.0⎦
===================== I3 =====================
2 2 2
cxx⋅cyy⋅czz - cxx⋅cyz - cxy ⋅czz + 2⋅cxy⋅cxz⋅cyz - cxz ⋅cyy
dI3dC
⎡ 2 ⎤
⎢ cyy⋅czz - cyz ⎥
⎢ ⎥
⎢ 2 ⎥
⎢ cxx⋅czz - cxz ⎥
⎢ ⎥
⎢ 2 ⎥
⎢ cxx⋅cyy - cxy ⎥
⎢ ⎥
⎢√2⋅(-cxx⋅cyz + cxy⋅cxz)⎥
⎢ ⎥
⎢√2⋅(cxy⋅cyz - cxz⋅cyy) ⎥
⎢ ⎥
⎣√2⋅(-cxy⋅czz + cxz⋅cyz)⎦
d2I3dC
⎡ 0 czz cyy -√2⋅cyz 0 0 ⎤
⎢ ⎥
⎢ czz 0 cxx 0 -√2⋅cxz 0 ⎥
⎢ ⎥
⎢ cyy cxx 0 0 0 -√2⋅cxy⎥
⎢ ⎥
⎢-√2⋅cyz 0 0 -cxx cxy cxz ⎥
⎢ ⎥
⎢ 0 -√2⋅cxz 0 cxy -cyy cyz ⎥
⎢ ⎥
⎣ 0 0 -√2⋅cxy cxz cyz -czz ⎦
===================== Ii =====================
[T1x⋅T2x⋅cxx + T1x⋅T2y⋅cxy + T1x⋅T2z⋅cxz + T1y⋅T2x⋅cxy + T1y⋅T2y⋅cyy + T1y⋅T2z ↪
↪ ⋅cyz + T1z⋅T2x⋅cxz + T1z⋅T2y⋅cyz + T1z⋅T2z⋅czz]
dIidC
⎡ T1x⋅T2x ⎤
⎢ ⎥
⎢ T1y⋅T2y ⎥
⎢ ⎥
⎢ T1z⋅T2z ⎥
⎢ ⎥
⎢ ⎛T1y⋅T2z T1z⋅T2y⎞⎥
⎢√2⋅⎜─────── + ───────⎟⎥
⎢ ⎝ 2 2 ⎠⎥
⎢ ⎥
⎢ ⎛T1x⋅T2z T1z⋅T2x⎞⎥
⎢√2⋅⎜─────── + ───────⎟⎥
⎢ ⎝ 2 2 ⎠⎥
⎢ ⎥
⎢ ⎛T1x⋅T2y T1y⋅T2x⎞⎥
⎢√2⋅⎜─────── + ───────⎟⎥
⎣ ⎝ 2 2 ⎠⎦
d2IidC
⎡0 0 0 0 0 0⎤
⎢ ⎥
⎢0 0 0 0 0 0⎥
⎢ ⎥
⎢0 0 0 0 0 0⎥
⎢ ⎥
⎢0 0 0 0 0 0⎥
⎢ ⎥
⎢0 0 0 0 0 0⎥
⎢ ⎥
⎣0 0 0 0 0 0⎦
13 from EasyFEA import Display
14
15 try:
16 import sympy
17 except ModuleNotFoundError:
18 raise Exception("sympy must be installed!")
19
20 from typing import Optional
21
22
23 def __Project_Mandel(A, orderA: int = 4):
24 assert orderA in [2, 4]
25
26 # for xx, yy, zz, yz, xz, zy
27 e = sympy.Array([[0, 5, 4], [5, 1, 3], [4, 3, 2]])
28
29 def kron(a, b):
30 return 1 if a == b else 0
31
32 if orderA == 2:
33 # Aij -> AI
34 A_I = sympy.zeros(6, 1)
35
36 for i in range(3):
37 for j in range(3):
38 A_I[e[i, j]] = sympy.sqrt((2 - kron(i, j))) * A[i, j]
39
40 res = A_I
41
42 elif orderA == 4:
43 # Aijkl -> AIJ
44
45 A_IJ = sympy.zeros(6, 6)
46
47 for i in range(3):
48 for j in range(3):
49 for k in range(3):
50 for l in range(3):
51 A_IJ[e[i, j], e[k, l]] = (
52 sympy.sqrt((2 - kron(i, j)) * (2 - kron(k, l)))
53 * A[i, j, k, l]
54 )
55
56 res = A_IJ
57
58 else:
59 raise Exception("Not implemented.")
60
61 return res
62
63
64 def __MyDiff(func, list_func: list, order: Optional[int] = None):
65 assert isinstance(list_func, list)
66
67 diff = sympy.diff(func, *list_func).as_mutable()
68
69 ndim = len(diff.shape)
70
71 e = sympy.Array([[0, 5, 4], [5, 1, 3], [4, 3, 2]])
72
73 # The use of a symmetrical matrix C on sympy is problematic here.
74 # The off-diagonal terms are incorrectly weighted and need to be corrected.
75 # The following function is used to correct the coefficients.
76 if ndim == 2:
77 for i in range(diff.shape[0]):
78 for j in range(diff.shape[1]):
79 if e[i, j] > 2:
80 diff[i, j] /= 2
81 elif ndim == 4:
82 for i in range(diff.shape[0]):
83 for j in range(diff.shape[1]):
84 for k in range(diff.shape[2]):
85 for l in range(diff.shape[3]):
86 if e[i, j] > 2 and e[k, l] > 2:
87 diff[i, j, k, l] /= 4
88 elif e[i, j] > 2 or e[k, l] > 2:
89 diff[i, j, k, l] /= 2
90
91 if order is not None:
92 diff = __Project_Mandel(diff, order)
93
94 diff = sympy.sympify(diff)
95
96 return diff
97
98
99 def Compute(func, name: str, usepprint=True):
100 Display.Section(name)
101
102 myprint = sympy.pprint if usepprint else print
103
104 func = sympy.sympify(sympy.expand(func))
105
106 myprint(func)
107
108 myprint(f"\nd{name}dC")
109 myprint(__MyDiff(func, [C], 2))
110
111 myprint(f"\nd2{name}dC")
112 myprint(__MyDiff(func, [C, C], 4))
113
114
115 if __name__ == "__main__":
116 Display.Clear()
117
118 cxx, cyy, czz, cyz, cxz, cxy = sympy.symbols("cxx, cyy, czz, cyz, cxz, cxy")
119 C = sympy.Matrix([[cxx, cxy, cxz], [cxy, cyy, cyz], [cxz, cyz, czz]])
120
121 sympy.pprint(C)
122 print()
123 sympy.pprint(__Project_Mandel(C, 2))
124
125 # -------------------------------------
126 # I1
127 # -------------------------------------
128
129 I1 = sympy.trace(C)
130
131 Compute(I1, "I1")
132
133 # -------------------------------------
134 # I2
135 # -------------------------------------
136
137 I2 = 1 / 2 * (sympy.trace(C) ** 2 - sympy.trace(C * C))
138
139 Compute(I2, "I2")
140
141 # -------------------------------------
142 # I3
143 # -------------------------------------
144
145 I3 = sympy.det(C)
146
147 Compute(I3, "I3")
148
149 # -------------------------------------
150 # I4, I6 and I8
151 # -------------------------------------
152
153 T1x, T1y, T1z = sympy.symbols("T1x, T1y, T1z")
154 T2x, T2y, T2z = sympy.symbols("T2x, T2y, T2z")
155
156 T1 = sympy.Matrix([[T1x], [T1y], [T1z]])
157 T2 = sympy.Matrix([[T2x], [T2y], [T2z]])
158
159 I8 = T1.transpose() * (C * T2)
160
161 Compute(I8, "Ii")
Total running time of the script: (0 minutes 0.366 seconds)