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