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

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

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.

Pictures:

Code (purely elastic):

Needs["NDSolveFEM"]
(*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["NDSolveFEM"]
(*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)

• Could you explain a bit what a user material subroutine is? Then I can perhaps help. – user21 May 6 '15 at 19:52
• You can set the coefficient matrices explicitly, if that is what you're after. See reference.wolfram.com/language/FEMDocumentation/tutorial/…. (You should respond to user21. He can really help usually.) – Michael E2 May 8 '15 at 20:35
• @user21 I will edit the question and add a 1D example. at Michael: not, that's not what I mean, I will explain it in the 1D example. – Mauricio Fernández May 10 '15 at 13:00
• I had a look at this and to me it still seams that these are just a method to specify your own PDE coefficients, which you can do in Mathematica.... – user21 Feb 4 '16 at 5:58
• @user21 yes, partly, but their use if also aimed at evolving internal variables and for fixed time step adapting the global jacobian accordingly. I made an edit almost at the beginning of the question and uploaded a hand written pdf with the main idea. If you want you can take a look at it. If this is too time consuming, I can understand that closing the question may be better. But thanks! – Mauricio Fernández Feb 4 '16 at 12:29

Note that the PDE coefficients can depend on space and time. Let me illustrate that with a 2D stationary example.

Here is a plain stress operator:

ps = {Inactive[
Div][{{0, -((Y*\[Nu])/(1 - \[Nu]^2))}, {-(Y*(1 - \[Nu]))/(2*(1 \
- \[Nu]^2)), 0}}.Inactive[Grad][v[x, y], {x, y}], {x, y}] +
Inactive[
Div][{{-(Y/(1 - \[Nu]^2)),
0}, {0, -(Y*(1 - \[Nu]))/(2*(1 - \[Nu]^2))}}.Inactive[Grad][
u[x, y], {x, y}], {x, y}],
Inactive[
Div][{{0, -(Y*(1 - \[Nu]))/(2*(1 - \[Nu]^2))}, {-((Y*\[Nu])/(1 \
- \[Nu]^2)), 0}}.Inactive[Grad][u[x, y], {x, y}], {x, y}] +
Inactive[
Div][{{-(Y*(1 - \[Nu]))/(2*(1 - \[Nu]^2)),
v[x, y], {x, y}], {x, y}]};


First we consider the case where Young's modulus is constant.

op = ps /. {Y -> 10^3, \[Nu] -> 33/100};
{uif, vif} = NDSolveValue[{op == {0, NeumannValue[-1, x == 5]},
DirichletCondition[{u[x, y] == 0., v[x, y] == 0.}, x == 0]
}, {u, v}, {x, 0, 5}, {y, 0, 1}];


The left hand side of the bar is clamped, the right hand side has a load is -1 unit in the y-direction. Visualize the deformation of the object:

mr = MeshRegion[uif["ElementMesh"]];
c = MeshCoordinates[mr];
Show[
BoundaryMesh[mr],
HighlightMesh[
MeshRegion[c + Transpose[{uif @@@ c, vif @@@ c}],
MeshCells[mr, {2, All}]], Style[1, Red]]
]


Now, in the second case we change the material model such that Young's modulus the space dependent:

op = ps /. {Y -> (10^3 + 10^3*x), \[Nu] -> 33/100};
{uif, vif} = NDSolveValue[{op == {0, NeumannValue[-1, x == 5]},
DirichletCondition[{u[x, y] == 0., v[x, y] == 0.}, x == 0]
}, {u, v}, {x, 0, 5}, {y, 0, 1}];


The rest of the setup remains the same.

Visualize the deformation:

mr = MeshRegion[uif["ElementMesh"]];
c = MeshCoordinates[mr];
Show[
BoundaryMesh[mr],
HighlightMesh[
MeshRegion[c + Transpose[{uif @@@ c, vif @@@ c}],
MeshCells[mr, {2, All}]], Style[1, Red]]
]


So, as long as the material model remains linear you can code use it. The change is simply done be specifying the PDE / model you want.

Hope this helps.

• hmmm... that's not my problem. It is clear that the coefficients of a linear PDE can vary. In my example, you are just changing Young's modulus E or the cross section area A in space. The problem above contains a second unknown, the plastic deformation epsilon_p. Its evolution is triggered only if the unknown displacement u induces stresses reaching the yield stress (phi=0). Otherwise, the unkown epsilon_p does not evolve or stops evoling. I'll try to write the complete system clearly and update the pdf and the pictures. – Mauricio Fernández May 11 '15 at 15:08
• One think that might be possible is to use WhenEvent then. – user21 May 11 '15 at 15:16
• Can you use that in the Mathematica FEM? May be in that case the complete system (see above in the updated pdf) can be just written down as a separated evolution equation. – Mauricio Fernández May 11 '15 at 15:28
• @MauricioLobos, version 10.3 WhenEvent and the FEM work together. I am still not 100% sure this is what you are after, but I wanted to mention that if it is, this should work now. – user21 Jan 14 '16 at 4:06
• Thanks for the information. I gave it a try now but I just dont know how to use it properly. Still thank you very much for the help. – Mauricio Fernández Jan 14 '16 at 14:44