8
$\begingroup$

The Solid Mechanics capabilities of Mathematica make it an interesting alternative to using other software especially when considering multi-physics problems. I would like to collect here (for the community) a collection of alternative materials models (e.g. hyper elastic models) for Finite Element Analysis. On one hand as a resource for others but also as a way for people to check to see if their models are actually correct. Some of these models are already built in to Mathematica, but it would be nice if the community could collect materials commonly used materials models here.

The currently available models are listed here. There are compressible and nearly incompressible models. The current models are:

There is also a model comparison and model calibration is dead easy. There is also a section on multiplicative decomposition (e.g. adding thermal and other strains)

A list of potential materials models could include the hyperelastic models (please edit add more if interested).

Viscoelastic models

Plasticity models (see e.g. here)

$\endgroup$
4
  • $\begingroup$ This is a heads up for version the upcoming 13.3. I have renames "MaterialModel" and "MaterialModelFunction" to "SolidMechanicsMaterialModel" and "SolidMechanicsMaterialModelFunction", respectively. Both "MaterialModel" and "MaterialModelFunction" will continue to work but are regarded legacy. This reason this was do is that then in the future we could have parameters with <| "SolidMechanicsMaterialModel" -> xyz, "AnotherDomainMaterialModel"-> zyx |> in the same parameter set up. $\endgroup$
    – user21
    Commented Apr 20, 2023 at 6:42
  • $\begingroup$ On a related note, I have also changes "ModelForm" to "SolidMechanicsModelForm" because "ModelForm" should be reserved for deciding if a model is conservative of not. Again, old code will continue to work. $\endgroup$
    – user21
    Commented Apr 20, 2023 at 6:43
  • 1
    $\begingroup$ Version 13.3 also has a nearly incompressible Yeoh model. More to come in 14.0. Also, the hyperelasticity monograph has a few updates in 13.3 $\endgroup$
    – user21
    Commented Jun 26, 2023 at 12:11
  • 1
    $\begingroup$ Version 14.0 has many, many new features for hyperelasticity. I took the liberty of updating the question with links to relevant sections. I hope you don't mind. Enjoy V14! $\endgroup$
    – user21
    Commented Jan 10 at 13:30

2 Answers 2

7
$\begingroup$

With the help of @user21 we have made an implementation of thermal expansion coupled to the NeoHookean model. The NeoHookean model used is the same as that implemented in Abaqus, which is slightly different to that implemented in Mathematica and Comsol.

MaterialModelNeoHookThermal[vars_, pars_, data__] :=
     Module[{U, X, dim, lambda, mu, gradU, f, idm, c, cinv, j,  stressMatrix, fInelastic, thermalstrain},

     lambda =   PDEModels`GetSolidMechanicsMaterialParameter[vars, pars,"LameParameter"];

     [FailureQ[lambda], Return[$Failed, Module];];
     mu = PDEModels`GetSolidMechanicsMaterialParameter[vars, pars, "ShearModulus"];
     If[FailureQ[mu], Return[$Failed, Module];];

     U = vars[[1]];
     X = vars[[-1]];
     dim = Length[X];
     idm = IdentityMatrix[3];
     gradU = ConstantArray[0, {3, 3}];

     gradU[[1 ;; dim, 1 ;; dim]] = Grad[U, X];
     f = gradU + idm;

     fInelastic = linearIsothermalstrain + idm;

     f = f . Inverse[fInelastic];

     If[pars["ModelForm"] == "PlaneStress", Print["Warning - PlaneStress not yet checked"];];

     j = Det[f];
     c = Transpose[f] . f;
     cinv = Inverse[c];

     (*energy density function based on Abaqus Neo-Hook*)
     stressMatrix = mu*j^(-2/3) (idm - (1/3) Tr[c] cinv) + lambda*(j^2 - j)*cinv;
     (*returns second Piola-Kirchhoff stress*)
     Simplify[stressMatrix[[1 ;; dim, 1 ;; dim]]]
   ]

The code was tested against data from Abaqus and matched well:

mesh = ToElementMesh[Rectangle[], MaxCellMeasure -> 0.01,"MeshElementType" -> "QuadElement"];

(*Properties*)
poisson = 0.3;
e = 1;

vars = {{u[x, y], v[x, y]}, {x, y}};

pars = <|"MaterialModelFunction" ->MaterialModelNeoHookThermal, 
 "ModelForm" -> "PlaneStrain",
 "LameParameter" -> e/(3 (1 - 2*poisson)), 
 "ShearModulus" -> e/(2 (1 + poisson)),
 "ThermalStrainTemperature" -> 1,
 "ThermalStrainReferenceTemperature" -> 0,
 "Thickness" -> 1, "MassDensity" -> 1|>;




linearIsothermalstrain = \[Alpha]*(1 + poisson)*IdentityMatrix[3];
If[Length[vars[[-1]]] == 2, linearIsothermalstrain[[-1, -1]] = 0];


 pdeModel = {SolidMechanicsPDEComponent[vars, pars] == {0, 
 0}, SolidDisplacementCondition[x == 0, vars, pars, <|"Displacement" -> {0, 0}|>]};

pfun =ParametricNDSolveValue[pdeModel, {u[x, y], v[x, y]}, {x, y} \[Element] mesh, \[Alpha]];

\[Alpha]Max = 0.2;
nsteps = 20;
displacementNeoHookThermal = {0, 0};

The model needs to be solved slowly by increasing the alpha parameter step wise.

Monitor[Do[\[Alpha]New = step*\[Alpha]Max/nsteps;
PrintTemporary["\[Alpha]New: ", \[Alpha]New];
displacementNeoHookThermalAbaqus = 
 pfunAbaqus[\[Alpha]New, 
  "InitialSeeding" -> 
   Thread[Equal[vars[[1]], displacementNeoHookThermalAbaqus]]];
, {step, 1, nsteps, 1}],
step] // AbsoluteTiming;

The following is a comparison of the same simulation run in Abaqus and Mathematica.

Comparison Mathematica Abaqus

Once again many thanks to @user21 for his help in getting this running.

$\endgroup$
4
$\begingroup$

Standard reinforcing material

The standard reinforcing material is simple and popular model of fiber-reinforced material based on the Neo-Hookean model with the following strain energy density function $W$: $$ W = \frac{G}{2} \left( \mathrm{Tr} \mathbb{C} - 3 \right) + \frac{\lambda}{2} \left( \log \det \mathbb{F} \right)^2 - G \det \mathbb{F} + \frac{G}{2} k \left( \mathbf{A} \cdot \mathbb{C} \mathbf{A} - 1 \right)^2, $$ where $\mathbb{F}$ is deformation gradient, $\mathbb{C}$ is right Cauchy-Green strain tensors, $\lambda$ is Lamé parameter, $G$ is shear modulus, $k \geq 0$ is parameter measuring the strength of fibers and $\mathbf{A}$ is unit vector describing the fiber direction in the reference configuration. For $k = 0$ we would obtain the standard Neo-Hookean model.

The Cauchy stress tensor $\mathbb{T}$ is given by the formula $$ \mathbb{T} = \frac{1}{\det \mathbb{F}} \left[ \left( \lambda \log \det \mathbb{F} - G \right) \mathbb{I} + G \mathbb{B} + 2 G k \left( \bf{A} \cdot \mathbb{C} \bf{A} - 1 \right) \mathbb{F} \bf{A} \otimes \bf{A} \mathbb{F}^\top \right], $$ where $\mathbb{B}$ is the left Cauchy-Green strain tensor.

StandardReinforcingMaterial[vars_, pars_, data_] := 
 Module[{u, x, dim, idm, G, \[Lambda], k, A, F, J, c, stressMatrix},
  u = vars[[1]]; 
  x = vars[[-1]];
  G = pars["ShearModulus"];
  \[Lambda] = pars["LaméParameter"];
  k = pars["k"];
  A = pars["A"];

  dim = Length[x]; 
  idm = IdentityMatrix[dim];

  F = Grad[u, x] + idm;
  J = Det[F];
  c = Transpose[F].F;

  stressMatrix = 1/J*((\[Lambda]*Log[J] - G)*idm + G*F.Transpose[F] + 2*G*k*(Transpose[A].c.A - 1)*(F.(A\[TensorProduct]A).Transpose[F]));
  stressMatrix = Simplify[stressMatrix];
  stressMatrix
 ]

Example:

ClearAll["Global`*"]
Needs["NDSolve`FEM`"]

\[CapitalOmega] = Rectangle[{0, 0}, {1, 1}];
mesh = ToElementMesh[\[CapitalOmega]]; 
vars = {{u[x, y], v[x, y]}, {x, y}}; 
pars = <|"MaterialModelFunction" -> StandardReinforcingMaterial, "Thickness" -> 1, "ConstitutiveStressMeasure" -> "Cauchy", "ShearModulus" -> 400, "LaméParameter" -> 270, "A" -> {Sqrt[2]/2, -Sqrt[2]/2}, "k" -> 1|>;

pde := {SolidMechanicsPDEComponent[vars, pars] == SolidBoundaryLoadValue[x == 1, vars, pars, <|"Pressure" -> {p, 0}|>], SolidFixedCondition[x == 0, vars, pars]};
displacement = NDSolveValue[pde /. p -> 400, {u[x, y], v[x, y]}, {x, y} \[Element] \[CapitalOmega]];

deformedMesh = ElementMeshDeformation[mesh, displacement, "ScalingFactor" -> 1];
VectorDisplacementPlot[displacement, {x, y} \[Element] \[CapitalOmega], VectorSizes -> Full]

And this is the result:

standard reinforcing material

Conjugate stress / strain basis model of biological membrane (2D)

This is a non-standard model of two dimensional transversely isotropic material that can be used to model e.g. biological membranes. This model is rather complicated as it uses the QR decomposition of the deformation gradient and requires definitions of two conjugate bases: one for strain and one for stress. Thus I won't explain the model and will just explain the input parameters and output. The explanation of the model can be found in this article.

The following material parameters are needed:

  • bulk modulus $K$,
  • P-wave modulus $M$,
  • shear modulus $G$,
  • extent of anisotropy $n$,
  • matrix of anisotropy $\mathbb{Q}$.

The extent of anistotropy measures the strength of fibers. For $n > 1$ the fibers are stronger than the matrix material, for $n < 1$ the fibers are weaker than the matrix material. The matrix of anisotropy defines the direction of fibers. It is a matrix of rotation that gives a tangent unit vector to the direction of fibers in reference configuration when multiplied by a unit vector in the $x$ direction.

The output is the Cauchy stress tensor.

QRMembrane[vars_, pars_, data_] := 
 Module[{u, x, dim, idm, n, Q, K0, M, G, F, FA, R, sin\[Theta], cos\[Theta], a, b, \[Gamma], U, \[Delta], \[Epsilon], pi, \[Sigma], \[Tau], STilde11, STilde22, STilde12, STilde, stressMatrix},
  u = vars[[1]];
  x = vars[[-1]]; 
  K0 = pars["BulkModulus"];
  M = pars["PWaveModulus"];
  G = pars["ShearModulus"];
  n = pars["n"];
  Q = pars["Q"];

  dim = Length[x];
  idm = IdentityMatrix[dim];

  F = idm + Grad[u, x];
  FA = Q.F.Transpose[Q];

  sin\[Theta] = -FA[[2, 1]]/Sqrt[FA[[1, 1]]^2 + FA[[2, 1]]^2];
  cos\[Theta] = FA[[1, 1]]/Sqrt[FA[[1, 1]]^2 + FA[[2, 1]]^2];
  R = {{cos\[Theta], sin\[Theta]}, {-sin\[Theta], cos\[Theta]}};

  a = Sqrt[FA[[1, 1]]^2 + FA[[2, 1]]^2];
  b = (FA[[1, 1]] FA[[2, 2]] - FA[[1, 2]] FA[[2, 1]])/Sqrt[FA[[1, 1]]^2 + FA[[2, 1]]^2];
  \[Gamma] = (FA[[1, 1]] FA[[1, 2]] + FA[[2, 1]] FA[[2, 2]])/(FA[[1, 1]]^2 + FA[[2, 1]]^2);
  U = {{a, a*\[Gamma]}, {0, b}};

  \[Delta] = Log[Sqrt[a^n*b^(1/n)]];
  \[Epsilon] = Log[Sqrt[(a^n)/(b^(1/n))]];

  pi = 4.*K0*\[Delta];
  \[Sigma] = 2.*M*\[Epsilon];
  \[Tau] = G*\[Gamma];

  STilde11 = 1./2.*(n*pi + n*\[Sigma]);
  STilde22 = 1./2.*(pi/n - \[Sigma]/n);
  STilde12 = (b/a)*\[Tau];
  STilde = {{STilde11, STilde12}, {STilde12, STilde22}};

  stressMatrix = 1/Det[F]*Transpose[Q].R.STilde.Transpose[R].Q;
  stressMatrix = Simplify[stressMatrix];
  stressMatrix
 ]

Example:

ClearAll["Global`*"]
Needs["NDSolve`FEM`"]

\[CapitalOmega] = Rectangle[{0, 0}, {1, 1}];
mesh = ToElementMesh[\[CapitalOmega]]; 
vars = {{u[x, y], v[x, y]}, {x, y}}; 
pars = <|"MaterialModelFunction" -> QRMembrane, "Thickness" -> 1, "ConstitutiveStressMeasure" -> "Cauchy", "BulkModulus" -> 1270, "PWaveModulus" -> 1670, "ShearModulus" -> 400, "Q" -> {{Sqrt[2]/2, -Sqrt[2]/2}, {Sqrt[2]/2, Sqrt[2]/2}}, "n" -> 1.5|>;

pde := {SolidMechanicsPDEComponent[vars, pars] == SolidBoundaryLoadValue[x == 1, vars, pars, <|"Pressure" -> {p, 0}|>], SolidFixedCondition[x == 0, vars, pars]};
displacement = NDSolveValue[pde /. p -> 400, {u[x, y], v[x, y]}, {x, y} \[Element] \[CapitalOmega]];

deformedMesh = ElementMeshDeformation[mesh, displacement, "ScalingFactor" -> 1];
VectorDisplacementPlot[displacement, {x, y} \[Element] \[CapitalOmega], VectorSizes -> Full]

This is the result:

conjugate stress / strain model

I would like to help @user21 for his help with the code.

$\endgroup$
2
  • 1
    $\begingroup$ (+1) This is great, thanks for sharing! If you get a chance adding a simple example for both cases would be useful, I think. $\endgroup$
    – user21
    Commented Mar 31, 2023 at 8:35
  • 1
    $\begingroup$ Thank you. No problem, I will update the answer with some examples. $\endgroup$
    – lemurman
    Commented Mar 31, 2023 at 9:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.