FEM: user material subroutines (UMAT in ABAQUS), mechanics

Question

Just a quick question (mainly for those who are also familiar with ABAQUS or other comparable FEM programs and non-linear modelling of mechanical behavior, e.g., plasticity). Is there any module in the FEM environment in Mathematica in which you can set a general material behavior like in the UMATs of ABAQUS? I am just thinking of a standard von Mises plasticity model with ideal plastic behavior and homogeneous isotropic elastic behavior for small deformations. Or would you just write the evolution equations of the state variables in the field equations for the FEM solver and that's it? Thank you!

EDIT (2016-02-04):

If this is getting to far or I just express my question way to unclear, I would understand if closing the question is better for the community. These

Google drive link to hand written pdf

are the main ideas I know of on how a UMAT works, e.g., in ABAQUS. They are used in order to evolve internal variables locally at each integration point of each element independently of the rest. The internal variables do not exist at the FEM nodes but only locally on the integration points. For a sequential global Newton interation, the tangent at each integration point must be calculated in the UMAT (user defined material, this comes from mechanics), which are then used in the chosen quadrature method for computing the integrals and the global Jacobian. The example given below could be treated with a UMAT by considering the field $\varepsilon_p$ as an internal variable and careful treating of the corresponding cases.

EDIT1 (2015-05-10):

Ok, I should have given a minimal problem, sorry for that. Here is the simplest 1D problem from mechanics I can come up with, see pictures below or the online pdf (link below). In ABAQUS for example the in general non-linear material model is called up as an external program. How would you treat this problem with the FEM of Mathematica?

EDIT2 (2015-05-11):

I just added the section "Complete system", in order to clarify the problem.

EDIT3 (2016-01-14):

I just edited the pdf and introduced the complete solution for a simpler case. @user21: I am tried the WhenEvent but I am not sure on how to combine it with the field equations (see code below, I actually already had problems without any evolution of $\varepsilon_p$). Sorry for leaving this thing hanging around for so long. I will try to work on the description of an UMAT in the next weeks.

pdf: Google drive

Pictures:

Code (purely elastic):

Needs["NDSolve`FEM`"]
(*Domain*)
xmin = 0;
xmax = 10;
reg = ImplicitRegion[xmin <= x <= xmax, {x}];
(*Field equation*)
sig[x_, t_] := Ef[x, t]*eps[x, t];
eps[x_, t_] := D[u[x, t], x];
phi[x_, t_] := Abs[sig[x, t]] - sig0;
feq = D[Af[x, t]*sig[x, t], x] + nf[x, t];
(*Paremeters*)
E0 = 3;
Ef[x_, t_] := E0;
A0 = 7;
Af[x_, t_] := A0;
sig0 = 2;
(*Inhomogeneity*)
nf[x_, t_] := 0;
F[t_] := A0*sig0*3*t;
(*Conditions*)
iconds = u[x, 0] == 0;
bconds = {
   DirichletCondition[u[x, t] == 0, x == xmin]
   (*,DirichletCondition[u[x,t]==1*t,x==xmax]*)
   };
nconds = NeumannValue[-F[t], x == xmax];
(*Fem solution*)
ufem = NDSolveValue[{feq == nconds, iconds, bconds}, u, 
   Element[x, reg], {t, 0, 1}];
(*Plot*)
Plot[Table[u[x, ti] /. u -> ufem, {ti, 0, 1, 0.2}], {x, xmin, xmax}, 
 AxesLabel -> {"x", "u(x,t)"}]

Evolution of displacement

Defective Code (attempt for use of WhenEvent):

Needs["NDSolve`FEM`"]
(*Domain*)
xmin = 0;
xmax = 10;
reg = ImplicitRegion[xmin <= x <= xmax, {x}];
(*Field equation*)
sig[x_, t_] := Ef[x, t]*(eps[x, t] - epsp[x, t]);
eps[x_, t_] := D[u[x, t], x];
phi[x_, t_] := Abs[sig[x, t]] - sig0;
feq = {
   D[Af[x, t]*sig[x, t], x] + nf[x, t]
   , D[epsp[x, t], t]
   (*,WhenEvent[phi[x,t]==0,D[epsp[x,t],t]->D[eps[x,t],
   t]]*)
   };
(*Paremeters*)
E0 = 3;
Ef[x_, t_] := E0;
A0 = 7;
Af[x_, t_] := A0;
sig0 = 2;
(*Inhomogeneity*)
nf[x_, t_] := 0;
F[t_] := 7*sig0*3*t;
(*Conditions*)
iconds = {u[x, 0] == 0, epsp[x, 0] == 0};
bconds = {
   DirichletCondition[u[x, t] == 0, x == xmin]
   (*,DirichletCondition[u[x,t]==1*t,x==xmax]*)
   };
nconds = NeumannValue[-F[t], x == xmax];
(*Fem solution*)
fem = NDSolveValue[{feq == {nconds, 0}, iconds, bconds}, {u, epsp}, 
   Element[x, reg], {t, 0, 1}];
(*Plot*)
Plot[Table[u[x, ti] /. u -> fem[[1]], {ti, 0, 1, 0.2}], {x, xmin, 
  xmax}, AxesLabel -> {"x", "u(x,t)"}]
Plot[Table[epsp[x, ti] /. epsp -> fem[[2]], {ti, 0, 1, 0.2}], {x, 
  xmin, xmax}, AxesLabel -> {"x", "epsp(x,t)"}]
Plot[epsp[x, 1] /. epsp -> fem[[2]], {x, xmin, xmax}, 
 PlotRange -> {-1, 1}, AxesLabel -> {"x", "epsp(x,t)"}]

Weird results for $\varepsilon_p$ (it should stay at zero, since no evolution is given in the code above, WhenEvent is not active since I dont know how to properly use it in this example)

Answer 1

In version 13.1 there is support for Hyperelastic material models. Currently implemented are a St. Venant-Kirchhoff model and a neo-Hookean model. The monograph above also shows how to add your own hyperelststic material model. Here is a simple example of how to use the predefined material model:

Set up the variables:

Needs["NDSolve`FEM`"]
vars = {{u[x, y, z], v[x, y, z], w[x, y, z]}, {x, y, z}};

Create a mesh:

cubeMesh = 
  ToElementMesh[Cuboid[], MaxCellMeasure -> 0.0025, 
   "MeshElementType" -> TetrahedronElement];

Set up the model:

parsStVKI = <|"MaterialModel" -> "VenantKirchhoffIsotropic", 
   "LameParameter" -> 100, "ShearModulus" -> 1|>;

Because the hyperelastic material models can be highly nonlinear (that's not a wolfram language problem but a general problem) we set up a parametric function to be able to restart the solver with updated initial seeds:

pfunStVKI = 
  ParametricNDSolveValue[{SolidMechanicsPDEComponent[vars, 
      parsStVKI] == {0, 0, 0}, 
    SolidFixedCondition[x == 0, vars, parsStVKI], 
    SolidDisplacementCondition[x == 1, vars, 
     parsStVKI, <|"Displacement" -> {p, None, None}|>]}, 
   vars[[1]], {x, y, z} \[Element] cubeMesh, p];

Set up parameters for the restart:

pMax = 0.2;
nsteps = 40;
displacementStVKI = {0, 0, 0};

Monitor the solution process:

Monitor[Do[pNew = step*pMax/nsteps;
   displacementStVKI = 
    pfunStVKI[pNew, 
     "InitialSeeding" -> 
      Thread[Equal[vars[[1]], displacementStVKI]]];, {step, 1, nsteps,
     1}], step];

Compute the deformation. Note that the scaling factor is set to 1 - this is the true deformation of the object:

deformedMeshStVKI = 
  ElementMeshDeformation[cubeMesh, displacementStVKI, 
   "ScalingFactor" -> 1];

Visualize the deformation:

Show[cubeMesh["Edgeframe"], 
 deformedMeshStVKI[
  "Wireframe"[
   "ElementMeshDirective" -> Directive[Opacity[0.2], Red]]]]

One caveat is that this solver is not the fastest. The main reason is that the wolfram language is missing some symbolic processing capabilities that I'd need. Once they are available I'll update the code.

Now, concerning plasticity, you'd need to couple the elasticity with an WhenEvent. I have not yet done that but once I get there I'll update the post (and the documentation)

Now, if you want to go ahead and create your own material model you can do so by specifying:

MaterialModel[varsIn_, parsIn_, data__] :=
 Module[
    {X, dim, lambda, mu, EE, ee, rules, stressMatrix, strainMeasure},
  Print["Hello World"];
    
  lambda = 
   PDEModels`GetSolidMechanicsMaterialParameter[varsIn, parsIn, 
    "LameParameter"];
    If[FailureQ[lambda], Return[$Failed, Module];];
  
    mu = PDEModels`GetSolidMechanicsMaterialParameter[varsIn, parsIn,  "ShearModulus"];
    If[FailureQ[mu], Return[$Failed, Module];];
  
    X = vars[[-1]];
    dim = Length[X];
    EE = Array[ee, {dim, dim}];
    stressMatrix = (lambda * Tr[EE] * IdentityMatrix[dim]) + 
    2 * mu * EE;
    strainMeasure = "StrainMeasure" /. data;
    rules = Thread[Flatten[EE] -> Flatten[strainMeasure]];
    stressMatrix = stressMatrix /. rules;
    stressMatrix
  ]

The you'd set up the parameters as follows:

parsMyModel = <|"MaterialModelFunction" -> MaterialModel,
   "LameParameter" -> 100, "ShearModulus" -> 1,
   "StrainMeasure" -> "GreenLagrange",
   "ConstitutiveStressMeasure" -> "SecondPiolaKirchhoff",
   "EquilibriumStressMeasure" -> "FirstPiolaKirchhoff"|>;

The "MaterialModelFunction" specifies the material model. The material parameters are the same as a above. As a "StrainMeasure" we use the "GreenLagrange" and then we need to tell SolidMechanicsPDEComponent what it is that we have specified. "ConstitutiveStressMeasure" is the stress measure used for the material law. "EquilibriumStressMeasure" is the stress measure of the nonlinear equilibrium equation. The "OutputStressMeasure", by default, is the "Cauchy" stress, but can be changed. The conversions between the stresses happens automatically.

The we proceed as before:

pfunMyModel = 
  ParametricNDSolveValue[{SolidMechanicsPDEComponent[vars, 
      parsMyModel] == {0, 0, 0}, 
    SolidFixedCondition[x == 0, vars, parsMyModel], 
    SolidDisplacementCondition[x == 1, vars, 
     parsMyModel, <|"Displacement" -> {p, None, None}|>]}, 
   vars[[1]], {x, y, z} \[Element] cubeMesh, p];

(* Hello World *)

And

pfunMyModel[0]

Will now return the interpolating functions for your material law.

And thanks a lot for all your help with the solid mechanics!