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))
=========== 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))
14 from EasyFEA import Display
15
16 try:
17 import sympy
18 except ModuleNotFoundError:
19 raise Exception("sympy must be installed!")
20
21
22 def Compute(W, params: list, details=True):
23 print(f"W = {W}\n")
24
25 # dW
26 dW = ""
27 for param_i in params:
28 p_i = str(param_i)
29 dWdIi = sympy.diff(W, param_i)
30 if dWdIi != 0:
31 dW += " + "
32 if details:
33 print(f"dWd{p_i} = {dWdIi}")
34 dW += f"dWd{p_i} * d{p_i}dC"
35 else:
36 dW += f"({dWdIi}) * d{p_i}dC"
37
38 dW = f"dW = 2 * ({dW})\n"
39 dW = dW.replace("+ -", "- ")
40 dW = dW.replace("( + ", "(")
41 print(dW)
42
43 # d2W
44 d2W1 = ""
45 d2W2 = ""
46
47 for param_i in params:
48 p_i = str(param_i)
49 dWdIi = sympy.diff(W, param_i)
50 if dWdIi != 0:
51 d2W1 += " + "
52 if details:
53 print(f"dWd{p_i} = {dWdIi}")
54 d2W1 += f"dWd{p_i} * d2{p_i}dC"
55 else:
56 d2W1 += f"({dWdIi}) * d2{p_i}dC"
57
58 for param_j in params:
59 p_j = str(param_j)
60 d2WdIiIj = sympy.diff(dWdIi, param_j)
61 if d2WdIiIj != 0:
62 d2W2 += " + "
63 if details:
64 print(f"d2Wd{p_i}d{p_j} = {d2WdIiIj}")
65 d2W2 += f"d2Wd{p_i}d{p_j} * TensorProd(d{p_i}dC, d{p_j}dC)"
66 else:
67 d2W2 += f"({d2WdIiIj}) * TensorProd(d{p_i}dC, d{p_j}dC)"
68
69 if d2W2 == "":
70 d2W = f"d2W = 4 * ({d2W1})"
71 else:
72 d2W = f"d2W = 4 * ({d2W1}) + 4 * ({d2W2})"
73 d2W = d2W.replace("+ -", "- ")
74 d2W = d2W.replace("( + ", "(")
75 print(d2W)
76
77
78 if __name__ == "__main__":
79 Display.Clear()
80
81 I1, I2, I3, I4, I6, I8 = sympy.symbols("I1, I2, I3, I4, I6, I8")
82
83 J1 = I1 * I3 ** (sympy.Rational(-1, 3))
84 J2 = I2 * I3 ** (sympy.Rational(-2, 3))
85 J = I3 ** (sympy.Rational(1, 2))
86
87 # -------------------------------------
88 # Neo-Hookean
89 # -------------------------------------
90
91 Display.Section("Neo-Hookean")
92
93 K = sympy.symbols("K")
94
95 W = K * (J1 - 3)
96
97 Compute(W, [I1, I2, I3])
98
99 # -------------------------------------
100 # Mooney-Rivlin
101 # -------------------------------------
102
103 Display.Section("Mooney-Rivlin")
104
105 K1, K2 = sympy.symbols("K1, K2")
106
107 W = K1 * (J1 - 3) + K2 * (J2 - 3) + K * (J - 1) ** 2
108
109 Compute(W, [I1, I2, I3])
110
111 # -------------------------------------
112 # Saint-Venant-Kirchhoff
113 # -------------------------------------
114
115 Display.Section("Saint-Venant-Kirchhoff")
116
117 lmbda, mu = sympy.symbols("lmbda, mu")
118
119 # W = lmbda/8 * (I1**2 - 6*I1 + 9) + mu/4 * (I1**2 - 2*I1 - 2*I2 + 3)
120 W = (
121 (lmbda / 8 + mu / 4) * I1**2
122 - mu * I2 / 2
123 - (3 * lmbda / 4 + mu / 2) * I1
124 + 9 * lmbda / 8
125 + 3 * mu / 4
126 + 1 / 2 * K * (I3 - 1) ** 2
127 )
128
129 Compute(W, [I1, I2, I3])
130
131 # -------------------------------------
132 # Holzapfel-Ogden
133 # -------------------------------------
134
135 Display.Section("Holzapfel-Ogden")
136
137 C0, C1, C2, C3, C4, C5, C6, C7 = sympy.symbols("C0:8")
138
139 ks = sympy.symbols("ks")
140 bulk, mu1, mu2 = sympy.symbols("bulk, mu1, mu2")
141
142 chi = lambda Ii: 1 / (1 + sympy.exp(-ks * (Ii - 1)))
143
144 W = (
145 C0 * (sympy.exp(C1 * (J1 - 3)) - 1)
146 + C2 * chi(I4) * (sympy.exp(C3 * (I4 - 1) ** 2) - 1)
147 + C4 * chi(I6) * (sympy.exp(C5 * (I6 - 1) ** 2) - 1)
148 + C6 * (sympy.exp(C7 * I8**2) - 1)
149 + bulk / 4 * (J**2 - 1 - 2 * sympy.ln(J))
150 + mu1 * (J1 - 3)
151 + mu2 * (J2 - 3)
152 )
153
154 Compute(W, [I1, I2, I3, I4, I6, I8])
Total running time of the script: (0 minutes 0.485 seconds)