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⎦
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, I6 and I8
150 # -------------------------------------
151
152 T1x, T1y, T1z = sympy.symbols("T1x, T1y, T1z")
153 T2x, T2y, T2z = sympy.symbols("T2x, T2y, T2z")
154
155 T1 = sympy.Matrix([[T1x], [T1y], [T1z]])
156 T2 = sympy.Matrix([[T2x], [T2y], [T2z]])
157
158 I8 = T1.transpose() * (C * T2)
159
160 Compute(I8, "Ii")
Total running time of the script: (0 minutes 0.714 seconds)