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