.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/Hyperelasticity/ComputeHyperelasticLaws.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_examples_Hyperelasticity_ComputeHyperelasticLaws.py: ComputeHyperelasticLaws ======================= Compute hyperelastic constitutive laws. .. GENERATED FROM PYTHON SOURCE LINES 11-153 .. rst-class:: sphx-glr-script-out .. code-block:: none ================= 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)) | .. code-block:: Python :lineno-start: 12 from EasyFEA import Display try: import sympy except ModuleNotFoundError: raise Exception("sympy must be installed!") def Compute(W, params: list, details=True): print(f"W = {W}\n") # dW dW = "" for param_i in params: p_i = str(param_i) dWdIi = sympy.diff(W, param_i) if dWdIi != 0: dW += " + " if details: print(f"dWd{p_i} = {dWdIi}") dW += f"dWd{p_i} * d{p_i}dC" else: dW += f"({dWdIi}) * d{p_i}dC" dW = f"dW = 2 * ({dW})\n" dW = dW.replace("+ -", "- ") dW = dW.replace("( + ", "(") print(dW) # d2W d2W1 = "" d2W2 = "" for param_i in params: p_i = str(param_i) dWdIi = sympy.diff(W, param_i) if dWdIi != 0: d2W1 += " + " if details: print(f"dWd{p_i} = {dWdIi}") d2W1 += f"dWd{p_i} * d2{p_i}dC" else: d2W1 += f"({dWdIi}) * d2{p_i}dC" for param_j in params: p_j = str(param_j) d2WdIiIj = sympy.diff(dWdIi, param_j) if d2WdIiIj != 0: d2W2 += " + " if details: print(f"d2Wd{p_i}d{p_j} = {d2WdIiIj}") d2W2 += f"d2Wd{p_i}d{p_j} * TensorProd(d{p_i}dC, d{p_j}dC)" else: d2W2 += f"({d2WdIiIj}) * TensorProd(d{p_i}dC, d{p_j}dC)" if d2W2 == "": d2W = f"d2W = 4 * ({d2W1})" else: d2W = f"d2W = 4 * ({d2W1}) + 4 * ({d2W2})" d2W = d2W.replace("+ -", "- ") d2W = d2W.replace("( + ", "(") print(d2W) if __name__ == "__main__": Display.Clear() I1, I2, I3, I4, I6, I8 = sympy.symbols("I1, I2, I3, I4, I6, I8") J1 = I1 * I3 ** (sympy.Rational(-1, 3)) J2 = I2 * I3 ** (sympy.Rational(-2, 3)) J = I3 ** (sympy.Rational(1, 2)) # ------------------------------------- # Neo-Hookean # ------------------------------------- Display.Section("Neo-Hookean") K = sympy.symbols("K") W = K * (J1 - 3) Compute(W, [I1, I2, I3]) # ------------------------------------- # Mooney-Rivlin # ------------------------------------- Display.Section("Mooney-Rivlin") K1, K2 = sympy.symbols("K1, K2") W = K1 * (J1 - 3) + K2 * (J2 - 3) + K * (J - 1) ** 2 Compute(W, [I1, I2, I3]) # ------------------------------------- # Saint-Venant-Kirchhoff # ------------------------------------- Display.Section("Saint-Venant-Kirchhoff") lmbda, mu = sympy.symbols("lmbda, mu") # W = lmbda/8 * (I1**2 - 6*I1 + 9) + mu/4 * (I1**2 - 2*I1 - 2*I2 + 3) W = ( (lmbda / 8 + mu / 4) * I1**2 - mu * I2 / 2 - (3 * lmbda / 4 + mu / 2) * I1 + 9 * lmbda / 8 + 3 * mu / 4 + 1 / 2 * K * (I3 - 1) ** 2 ) Compute(W, [I1, I2, I3]) # ------------------------------------- # Holzapfel-Ogden # ------------------------------------- Display.Section("Holzapfel-Ogden") C0, C1, C2, C3, C4, C5, C6, C7 = sympy.symbols("C0:8") ks = sympy.symbols("ks") bulk, mu1, mu2 = sympy.symbols("bulk, mu1, mu2") chi = lambda Ii: 1 / (1 + sympy.exp(-ks * (Ii - 1))) W = ( C0 * (sympy.exp(C1 * (J1 - 3)) - 1) + C2 * chi(I4) * (sympy.exp(C3 * (I4 - 1) ** 2) - 1) + C4 * chi(I6) * (sympy.exp(C5 * (I6 - 1) ** 2) - 1) + C6 * (sympy.exp(C7 * I8**2) - 1) + bulk / 4 * (J**2 - 1 - 2 * sympy.ln(J)) + mu1 * (J1 - 3) + mu2 * (J2 - 3) ) Compute(W, [I1, I2, I3, I4, I6, I8]) .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.374 seconds) .. _sphx_glr_download_examples_Hyperelasticity_ComputeHyperelasticLaws.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: ComputeHyperelasticLaws.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: ComputeHyperelasticLaws.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: ComputeHyperelasticLaws.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_