Note
Go to the end to download the full example code.
ComputeHyperelasticLaws#
Compute hyperelastic constitutive laws.
================= Neo-Hookean =================
W = K*(I1/I3**(1/3) - 3)
dWdI1 = K/I3**(1/3)
dWdI3 = -I1*K/(3*I3**(4/3))
dW = 2 * (dWdI1 * dI1dC + dWdI3 * dI3dC)
dWdI1 = K/I3**(1/3)
d2WdI1dI3 = -K/(3*I3**(4/3))
dWdI3 = -I1*K/(3*I3**(4/3))
d2WdI3dI1 = -K/(3*I3**(4/3))
d2WdI3dI3 = 4*I1*K/(9*I3**(7/3))
d2W = 4 * (dWdI1 * d2I1dC + dWdI3 * d2I3dC) + 4 * (d2WdI1dI3 * TensorProd(dI1dC, dI3dC) + d2WdI3dI1 * TensorProd(dI3dC, dI1dC) + d2WdI3dI3 * TensorProd(dI3dC, dI3dC))
================ Mooney-Rivlin ================
W = K*(sqrt(I3) - 1)**2 + K1*(I1/I3**(1/3) - 3) + K2*(I2/I3**(2/3) - 3)
dWdI1 = K1/I3**(1/3)
dWdI2 = K2/I3**(2/3)
dWdI3 = -I1*K1/(3*I3**(4/3)) - 2*I2*K2/(3*I3**(5/3)) + K*(sqrt(I3) - 1)/sqrt(I3)
dW = 2 * (dWdI1 * dI1dC + dWdI2 * dI2dC + dWdI3 * dI3dC)
dWdI1 = K1/I3**(1/3)
d2WdI1dI3 = -K1/(3*I3**(4/3))
dWdI2 = K2/I3**(2/3)
d2WdI2dI3 = -2*K2/(3*I3**(5/3))
dWdI3 = -I1*K1/(3*I3**(4/3)) - 2*I2*K2/(3*I3**(5/3)) + K*(sqrt(I3) - 1)/sqrt(I3)
d2WdI3dI1 = -K1/(3*I3**(4/3))
d2WdI3dI2 = -2*K2/(3*I3**(5/3))
d2WdI3dI3 = 4*I1*K1/(9*I3**(7/3)) + 10*I2*K2/(9*I3**(8/3)) + K/(2*I3) - K*(sqrt(I3) - 1)/(2*I3**(3/2))
d2W = 4 * (dWdI1 * d2I1dC + dWdI2 * d2I2dC + dWdI3 * d2I3dC) + 4 * (d2WdI1dI3 * TensorProd(dI1dC, dI3dC) + d2WdI2dI3 * TensorProd(dI2dC, dI3dC) + d2WdI3dI1 * TensorProd(dI3dC, dI1dC) + d2WdI3dI2 * TensorProd(dI3dC, dI2dC) + d2WdI3dI3 * TensorProd(dI3dC, dI3dC))
============== Ciarlet-Geymonat ==============
W = K*(sqrt(I3) - log(sqrt(I3)) - 1) + K1*(I1/I3**(1/3) - 3) + K2*(I2/I3**(2/3) - 3)
dWdI1 = K1/I3**(1/3)
dWdI2 = K2/I3**(2/3)
dWdI3 = -I1*K1/(3*I3**(4/3)) - 2*I2*K2/(3*I3**(5/3)) + K*(-1/(2*I3) + 1/(2*sqrt(I3)))
dW = 2 * (dWdI1 * dI1dC + dWdI2 * dI2dC + dWdI3 * dI3dC)
dWdI1 = K1/I3**(1/3)
d2WdI1dI3 = -K1/(3*I3**(4/3))
dWdI2 = K2/I3**(2/3)
d2WdI2dI3 = -2*K2/(3*I3**(5/3))
dWdI3 = -I1*K1/(3*I3**(4/3)) - 2*I2*K2/(3*I3**(5/3)) + K*(-1/(2*I3) + 1/(2*sqrt(I3)))
d2WdI3dI1 = -K1/(3*I3**(4/3))
d2WdI3dI2 = -2*K2/(3*I3**(5/3))
d2WdI3dI3 = 4*I1*K1/(9*I3**(7/3)) + 10*I2*K2/(9*I3**(8/3)) + K*(1/(2*I3**2) - 1/(4*I3**(3/2)))
d2W = 4 * (dWdI1 * d2I1dC + dWdI2 * d2I2dC + dWdI3 * d2I3dC) + 4 * (d2WdI1dI3 * TensorProd(dI1dC, dI3dC) + d2WdI2dI3 * TensorProd(dI2dC, dI3dC) + d2WdI3dI1 * TensorProd(dI3dC, dI1dC) + d2WdI3dI2 * TensorProd(dI3dC, dI2dC) + d2WdI3dI3 * TensorProd(dI3dC, dI3dC))
=========== Saint-Venant-Kirchhoff ===========
W = I1**2*(lmbda/8 + mu/4) - I1*(3*lmbda/4 + mu/2) - I2*mu/2 + 0.5*K*(I3 - 1)**2 + 9*lmbda/8 + 3*mu/4
dWdI1 = 2*I1*(lmbda/8 + mu/4) - 3*lmbda/4 - mu/2
dWdI2 = -mu/2
dWdI3 = 0.5*K*(2*I3 - 2)
dW = 2 * (dWdI1 * dI1dC + dWdI2 * dI2dC + dWdI3 * dI3dC)
dWdI1 = 2*I1*(lmbda/8 + mu/4) - 3*lmbda/4 - mu/2
d2WdI1dI1 = lmbda/4 + mu/2
dWdI2 = -mu/2
dWdI3 = 0.5*K*(2*I3 - 2)
d2WdI3dI3 = 1.0*K
d2W = 4 * (dWdI1 * d2I1dC + dWdI2 * d2I2dC + dWdI3 * d2I3dC) + 4 * (d2WdI1dI1 * TensorProd(dI1dC, dI1dC) + d2WdI3dI3 * TensorProd(dI3dC, dI3dC))
=============== Holzapfel-Ogden ===============
W = C0*(exp(C1*(I1/I3**(1/3) - 3)) - 1) + C2*(exp(C3*(I4 - 1)**2) - 1)/(1 + exp(-ks*(I4 - 1))) + C4*(exp(C5*(I6 - 1)**2) - 1)/(1 + exp(-ks*(I6 - 1))) + C6*(exp(C7*I8**2) - 1) + bulk*(I3 - 2*log(sqrt(I3)) - 1)/4 + mu1*(I1/I3**(1/3) - 3) + mu2*(I2/I3**(2/3) - 3)
dWdI1 = C0*C1*exp(C1*(I1/I3**(1/3) - 3))/I3**(1/3) + mu1/I3**(1/3)
dWdI2 = mu2/I3**(2/3)
dWdI3 = -C0*C1*I1*exp(C1*(I1/I3**(1/3) - 3))/(3*I3**(4/3)) - I1*mu1/(3*I3**(4/3)) - 2*I2*mu2/(3*I3**(5/3)) + bulk*(1 - 1/I3)/4
dWdI4 = C2*C3*(2*I4 - 2)*exp(C3*(I4 - 1)**2)/(1 + exp(-ks*(I4 - 1))) + C2*ks*(exp(C3*(I4 - 1)**2) - 1)*exp(-ks*(I4 - 1))/(1 + exp(-ks*(I4 - 1)))**2
dWdI6 = C4*C5*(2*I6 - 2)*exp(C5*(I6 - 1)**2)/(1 + exp(-ks*(I6 - 1))) + C4*ks*(exp(C5*(I6 - 1)**2) - 1)*exp(-ks*(I6 - 1))/(1 + exp(-ks*(I6 - 1)))**2
dWdI8 = 2*C6*C7*I8*exp(C7*I8**2)
dW = 2 * (dWdI1 * dI1dC + dWdI2 * dI2dC + dWdI3 * dI3dC + dWdI4 * dI4dC + dWdI6 * dI6dC + dWdI8 * dI8dC)
dWdI1 = C0*C1*exp(C1*(I1/I3**(1/3) - 3))/I3**(1/3) + mu1/I3**(1/3)
d2WdI1dI1 = C0*C1**2*exp(C1*(I1/I3**(1/3) - 3))/I3**(2/3)
d2WdI1dI3 = -C0*C1**2*I1*exp(C1*(I1/I3**(1/3) - 3))/(3*I3**(5/3)) - C0*C1*exp(C1*(I1/I3**(1/3) - 3))/(3*I3**(4/3)) - mu1/(3*I3**(4/3))
dWdI2 = mu2/I3**(2/3)
d2WdI2dI3 = -2*mu2/(3*I3**(5/3))
dWdI3 = -C0*C1*I1*exp(C1*(I1/I3**(1/3) - 3))/(3*I3**(4/3)) - I1*mu1/(3*I3**(4/3)) - 2*I2*mu2/(3*I3**(5/3)) + bulk*(1 - 1/I3)/4
d2WdI3dI1 = -C0*C1**2*I1*exp(C1*(I1/I3**(1/3) - 3))/(3*I3**(5/3)) - C0*C1*exp(C1*(I1/I3**(1/3) - 3))/(3*I3**(4/3)) - mu1/(3*I3**(4/3))
d2WdI3dI2 = -2*mu2/(3*I3**(5/3))
d2WdI3dI3 = C0*C1**2*I1**2*exp(C1*(I1/I3**(1/3) - 3))/(9*I3**(8/3)) + 4*C0*C1*I1*exp(C1*(I1/I3**(1/3) - 3))/(9*I3**(7/3)) + 4*I1*mu1/(9*I3**(7/3)) + 10*I2*mu2/(9*I3**(8/3)) + bulk/(4*I3**2)
dWdI4 = C2*C3*(2*I4 - 2)*exp(C3*(I4 - 1)**2)/(1 + exp(-ks*(I4 - 1))) + C2*ks*(exp(C3*(I4 - 1)**2) - 1)*exp(-ks*(I4 - 1))/(1 + exp(-ks*(I4 - 1)))**2
d2WdI4dI4 = C2*C3**2*(2*I4 - 2)**2*exp(C3*(I4 - 1)**2)/(1 + exp(-ks*(I4 - 1))) + 2*C2*C3*ks*(2*I4 - 2)*exp(C3*(I4 - 1)**2)*exp(-ks*(I4 - 1))/(1 + exp(-ks*(I4 - 1)))**2 + 2*C2*C3*exp(C3*(I4 - 1)**2)/(1 + exp(-ks*(I4 - 1))) - C2*ks**2*(exp(C3*(I4 - 1)**2) - 1)*exp(-ks*(I4 - 1))/(1 + exp(-ks*(I4 - 1)))**2 + 2*C2*ks**2*(exp(C3*(I4 - 1)**2) - 1)*exp(-2*ks*(I4 - 1))/(1 + exp(-ks*(I4 - 1)))**3
dWdI6 = C4*C5*(2*I6 - 2)*exp(C5*(I6 - 1)**2)/(1 + exp(-ks*(I6 - 1))) + C4*ks*(exp(C5*(I6 - 1)**2) - 1)*exp(-ks*(I6 - 1))/(1 + exp(-ks*(I6 - 1)))**2
d2WdI6dI6 = C4*C5**2*(2*I6 - 2)**2*exp(C5*(I6 - 1)**2)/(1 + exp(-ks*(I6 - 1))) + 2*C4*C5*ks*(2*I6 - 2)*exp(C5*(I6 - 1)**2)*exp(-ks*(I6 - 1))/(1 + exp(-ks*(I6 - 1)))**2 + 2*C4*C5*exp(C5*(I6 - 1)**2)/(1 + exp(-ks*(I6 - 1))) - C4*ks**2*(exp(C5*(I6 - 1)**2) - 1)*exp(-ks*(I6 - 1))/(1 + exp(-ks*(I6 - 1)))**2 + 2*C4*ks**2*(exp(C5*(I6 - 1)**2) - 1)*exp(-2*ks*(I6 - 1))/(1 + exp(-ks*(I6 - 1)))**3
dWdI8 = 2*C6*C7*I8*exp(C7*I8**2)
d2WdI8dI8 = 4*C6*C7**2*I8**2*exp(C7*I8**2) + 2*C6*C7*exp(C7*I8**2)
d2W = 4 * (dWdI1 * d2I1dC + dWdI2 * d2I2dC + dWdI3 * d2I3dC + dWdI4 * d2I4dC + dWdI6 * d2I6dC + dWdI8 * d2I8dC) + 4 * (d2WdI1dI1 * TensorProd(dI1dC, dI1dC) + d2WdI1dI3 * TensorProd(dI1dC, dI3dC) + d2WdI2dI3 * TensorProd(dI2dC, dI3dC) + d2WdI3dI1 * TensorProd(dI3dC, dI1dC) + d2WdI3dI2 * TensorProd(dI3dC, dI2dC) + d2WdI3dI3 * TensorProd(dI3dC, dI3dC) + d2WdI4dI4 * TensorProd(dI4dC, dI4dC) + d2WdI6dI6 * TensorProd(dI6dC, dI6dC) + d2WdI8dI8 * TensorProd(dI8dC, dI8dC))
15 from EasyFEA import Display
16
17 try:
18 import sympy
19 except ModuleNotFoundError:
20 raise Exception("sympy must be installed!")
21
22
23 def Compute(W, params: list, details=True):
24 print(f"W = {W}\n")
25
26 # dW
27 dW = ""
28 for param_i in params:
29 p_i = str(param_i)
30 dWdIi = sympy.diff(W, param_i)
31 if dWdIi != 0:
32 dW += " + "
33 if details:
34 print(f"dWd{p_i} = {dWdIi}")
35 dW += f"dWd{p_i} * d{p_i}dC"
36 else:
37 dW += f"({dWdIi}) * d{p_i}dC"
38
39 dW = f"dW = 2 * ({dW})\n"
40 dW = dW.replace("+ -", "- ")
41 dW = dW.replace("( + ", "(")
42 print(dW)
43
44 # d2W
45 d2W1 = ""
46 d2W2 = ""
47
48 for param_i in params:
49 p_i = str(param_i)
50 dWdIi = sympy.diff(W, param_i)
51 if dWdIi != 0:
52 d2W1 += " + "
53 if details:
54 print(f"dWd{p_i} = {dWdIi}")
55 d2W1 += f"dWd{p_i} * d2{p_i}dC"
56 else:
57 d2W1 += f"({dWdIi}) * d2{p_i}dC"
58
59 for param_j in params:
60 p_j = str(param_j)
61 d2WdIiIj = sympy.diff(dWdIi, param_j)
62 if d2WdIiIj != 0:
63 d2W2 += " + "
64 if details:
65 print(f"d2Wd{p_i}d{p_j} = {d2WdIiIj}")
66 d2W2 += f"d2Wd{p_i}d{p_j} * TensorProd(d{p_i}dC, d{p_j}dC)"
67 else:
68 d2W2 += f"({d2WdIiIj}) * TensorProd(d{p_i}dC, d{p_j}dC)"
69
70 if d2W2 == "":
71 d2W = f"d2W = 4 * ({d2W1})"
72 else:
73 d2W = f"d2W = 4 * ({d2W1}) + 4 * ({d2W2})"
74 d2W = d2W.replace("+ -", "- ")
75 d2W = d2W.replace("( + ", "(")
76 print(d2W)
77
78
79 if __name__ == "__main__":
80 Display.Clear()
81
82 I1, I2, I3, I4, I6, I8 = sympy.symbols("I1, I2, I3, I4, I6, I8")
83
84 J1 = I1 * I3 ** (sympy.Rational(-1, 3))
85 J2 = I2 * I3 ** (sympy.Rational(-2, 3))
86 J = I3 ** (sympy.Rational(1, 2))
87
88 # -------------------------------------
89 # Neo-Hookean
90 # -------------------------------------
91
92 Display.Section("Neo-Hookean")
93
94 K = sympy.symbols("K")
95
96 W = K * (J1 - 3)
97
98 Compute(W, [I1, I2, I3])
99
100 # -------------------------------------
101 # Mooney-Rivlin
102 # -------------------------------------
103
104 Display.Section("Mooney-Rivlin")
105
106 K1, K2 = sympy.symbols("K1, K2")
107
108 W = K1 * (J1 - 3) + K2 * (J2 - 3) + K * (J - 1) ** 2
109
110 Compute(W, [I1, I2, I3])
111
112 # -------------------------------------
113 # Ciarlet-Geymonat
114 # -------------------------------------
115
116 Display.Section("Ciarlet-Geymonat")
117
118 W = K1 * (J1 - 3) + K2 * (J2 - 3) + K * (J - 1 - sympy.log(J))
119
120 Compute(W, [I1, I2, I3])
121
122 # -------------------------------------
123 # Saint-Venant-Kirchhoff
124 # -------------------------------------
125
126 Display.Section("Saint-Venant-Kirchhoff")
127
128 lmbda, mu = sympy.symbols("lmbda, mu")
129
130 # W = lmbda/8 * (I1**2 - 6*I1 + 9) + mu/4 * (I1**2 - 2*I1 - 2*I2 + 3)
131 W = (
132 (lmbda / 8 + mu / 4) * I1**2
133 - mu * I2 / 2
134 - (3 * lmbda / 4 + mu / 2) * I1
135 + 9 * lmbda / 8
136 + 3 * mu / 4
137 + 1 / 2 * K * (I3 - 1) ** 2
138 )
139
140 Compute(W, [I1, I2, I3])
141
142 # -------------------------------------
143 # Holzapfel-Ogden
144 # -------------------------------------
145
146 Display.Section("Holzapfel-Ogden")
147
148 C0, C1, C2, C3, C4, C5, C6, C7 = sympy.symbols("C0:8")
149
150 ks = sympy.symbols("ks")
151 bulk, mu1, mu2 = sympy.symbols("bulk, mu1, mu2")
152
153 chi = lambda Ii: 1 / (1 + sympy.exp(-ks * (Ii - 1)))
154
155 W = (
156 C0 * (sympy.exp(C1 * (J1 - 3)) - 1)
157 + C2 * chi(I4) * (sympy.exp(C3 * (I4 - 1) ** 2) - 1)
158 + C4 * chi(I6) * (sympy.exp(C5 * (I6 - 1) ** 2) - 1)
159 + C6 * (sympy.exp(C7 * I8**2) - 1)
160 + bulk / 4 * (J**2 - 1 - 2 * sympy.ln(J))
161 + mu1 * (J1 - 3)
162 + mu2 * (J2 - 3)
163 )
164
165 Compute(W, [I1, I2, I3, I4, I6, I8])
Total running time of the script: (0 minutes 2.130 seconds)