finite element method – Large deformation of solids

Link to notebook with this question and code

I’d like to understand how large deformations of solid mechanics work and how they are implemented. For this am looking at the following reference problem: A long slender beam is fixated at the left. A larger force is pushing in the negative x direction and a smaller force in the negative y direction.

enter image description here

To keep things a simple as possible this example uses a linear plane stress model (a linear material model) – it’s just about the large deformation.

A reference to this benchmark problem can be found here. There is also a video starting at time 28:00. Geometric nonlinearity is explained at 26:22.

The result for the plane stress case is a largely deformed beam. In the picture below (taken from the source given above) we are talking about the lower of the two beams.

enter image description here

Set up some parameters.

length = 3.2;
hight = 1/10;
thickness = 1/10;
(CapitalOmega) = Rectangle({0, 0}, {length, hight});
mesh = ToElementMesh((CapitalOmega));
bmesh = ToBoundaryMesh(mesh);

The total force is given. I am not 100% sure what total force means in Comsol speak, but I assume it means force per area. The force in the negative x direction is 10^3 times larger then the one in the negative y direction. We use Poisson’s ratio of 0 and Young’s modulus of 210 GPa.

(*force per area*)
fx = -3.844*10^6/(hight*thickness);
fy = fx*10^-3;
bc2 = {NeumannValue(fx, x == length), NeumannValue(fy, x == length)};
bc1 = DirichletCondition({u(x, y) == 0, v(x, y) == 0}, x == 0);
materialData = {nu -> 0, Y -> 210*10^9};

We start by using a linear – no large deformation – plane stress operator used on SE in several places and known to produce correct results. The only thing I changed is adding a ‘thickness’ factor to account for the thickness that is not 1.

PlaneStress({Y_, nu_}, {u_, v_}, X : {x_, y_}) := Module({pStress},
  (* uses global thickness *)
  pStress = -thickness * 
    Y/(1 - nu^2)*{{{{1, 0}, {0, (1 - nu)/2}}, {{0, nu}, {(1 - nu)/2, 
        0}}}, {{{0, (1 - nu)/2}, {nu, 0}}, {{(1 - nu)/2, 0}, {0, 1}}}};
  {Inactive(Div)(pStress((1, 1)) . Inactive(Grad)(u, X), X) + 
    Inactive(Div)(pStress((1, 2)) . Inactive(Grad)(v, X), X), 
   Inactive(Div)(pStress((2, 1)) . Inactive(Grad)(u, X), X) + 
    Inactive(Div)(pStress((2, 2)) . Inactive(Grad)(v, X), X)})

Solve the equation:

planeStress = 
  PlaneStress({Y, nu} /. materialData, {u(x, y), v(x, y)}, {x, y});
{uif2, vif2} = 
  NDSolveValue({planeStress == bc2, bc1}, {u, 
    v}, {x, y} (Element) (CapitalOmega));
 deform2 = 
  ElementMeshDeformation(mesh, {uif2, vif2}, "ScalingFactor" -> 1)(
   "Wireframe"("MeshElementStyle" -> EdgeForm(Red))))

Next, I am solving the same problem with a different setup. This setup will allow me later to easily change to a large deformation setup. We use the following:

enter image description here

This is then converted to a 2 by 2 matrix:

  Thread(Rule({{1, 1}, {2, 2}, {1, 2}}, {Subscript((Sigma), xx), 
     Subscript((Sigma), yy), Subscript((Tau), xy)})), Automatic, 

(* {{Subscript((Sigma), xx), Subscript((Tau), 
  xy)}, {Subscript((Tau), xy), Subscript((Sigma), yy)}} *)

And plugged into (with an additional thickness factor):

 Inactive(Div)(-{{Subscript((Sigma), xx), Subscript((Tau), 
      xy)}, {Subscript((Tau), xy), Subscript((Sigma), 
      yy)}}((i, ;;)), X), {i, Length(U)})

(* {Inactive(
   Div)({-Subscript((Sigma), xx), -Subscript((Tau), xy)}, {x, y}), 
 Inactive(Div)({-Subscript((Tau), xy), -Subscript((Sigma), yy)}, {x,
    y})} *)

makeStressModel(e_) := Block(
  {materialData, materialModel, exx, eyy, gxy, eVec, strain, nu, Y, 
  (* linear elastic plane stress model *)
  materialModel = 
   Y/(1 - nu^2)*{{1, nu, 0}, {nu, 1, 0}, {0, 0, (1 - nu)/2}};
  exx = e((1, 1));
  eyy = e((2, 2));
  gxy = (e((1, 2)) + e((2, 1)));
  (* for plane stress ezz =nu/(1-nu)*(exx+
  eyy) does not play a role since nu = 0 *)
  strain = {exx, eyy, gxy};
  stress = (materialModel . strain);
    Thread(Rule({{1, 1}, {2, 2}, {1, 2}}, stress)), Automatic, 

(* uses global defined thickness *)
makePDEModel(stressMatrix_) := 
 Table(Inactive(Div)(-thickness*stressMatrix((i, ;;)), X), {i, 

We verify that this gives the same results as a the previous plane stress operator:

X = {x, y};
U = {u(x, y), v(x, y)};

(* geometric linear *)
gu = Grad(U, X);
gut = Transpose(gu);
e = 1/2 (gu + gut);
stress = makeStressModel(e);

pde = makePDEModel(stress /. materialData);
{uif, vif} = 
  NDSolveValue({pde == bc2, bc1}, {u, v}, {x, y} (Element) mesh);

 deform1 = 
  ElementMeshDeformation(mesh, {uif, vif}, "ScalingFactor" -> 1)(
   "Wireframe"("MeshElementStyle" -> EdgeForm(Brown))))

enter image description here

Norm(uif("ValuesOnGrid") - uif2("ValuesOnGrid"))
Norm(vif("ValuesOnGrid") - vif2("ValuesOnGrid"))

(* 1.43378*10^-10

With the framework in place and verified, I want to look at the large deformations. I set up the deformation gradient f and set up the non-linear strain.

(* geometric nonlinear *)
f = Grad(U, {x, y}) + IdentityMatrix(2);
e1 = 1/2 (Transpose(f) . f - IdentityMatrix(2)) // Expand // Simplify;

With the same approach as above

stress = makeStressModel(e1);
pde = makePDEModel(stress /. materialData);
{uif, vif} = 
  NDSolveValue({pde == bc2, bc1}, {u, 
    v}, {x, y} (Element) (CapitalOmega));
mesh = uif("ElementMesh");
 deform2 = 
  ElementMeshDeformation(mesh, {uif, vif}, "ScalingFactor" -> 1)(
   "Wireframe"("MeshElementStyle" -> EdgeForm(Red))))

enter image description here

Show(deform1, deform2)

enter image description here

This is not correct. I would have expected a failure from the nonlinear solver for the large deformation or a much larger deformation. I am trying to understand what I am missing.
Question: Why is this approach wrong or where did I make a mistake?

A second way to specify the large strains is by

gu = Grad(U, X);
gut = Transpose(gu);
e2 = 1/2 (gu + gut + gut . gu);
FullSimplify(e1 - e2)

(* {{0, 0}, {0, 0}} *)

But it gives the same result. I have verified that the PDE coefficients parse correctly

{state} = 
  NDSolve`ProcessEquations({pde == bc2, bc1}, {u, 
    v}, {x, y} (Element) (CapitalOmega));
GetInactivePDE(state) - pde

(*{0, 0}*)

The next idea was that using the Cauchy stress is not correct. So I changed to the second Piola-Kirchhoff stress. We compute the stress as before

stress = makeStressModel(e1);

Now, find the deformation gradient, it’s inverse and transpose and J the determinant of f.

f = Grad(U, X) + IdentityMatrix(2);
invF = Inverse(f);
invFT = Transpose(invF);
piolaStress = Det(f)*invF . stress . invFT;
pde = makePDEModel(piolaStress /. materialData);
{uif, vif} = 
  NDSolveValue({pde == bc2, bc1}, {u, 
    v}, {x, y} (Element) (CapitalOmega));
deform3 = 
  ElementMeshDeformation(mesh, {uif, vif}, "ScalingFactor" -> 1)(
   "Wireframe"("MeshElementStyle" -> EdgeForm(Blue)));
Show(deform2, deform3)

enter image description here

Still, no joy.

{state} = 
  NDSolve`ProcessEquations({pde == bc2, bc1}, {u, 
    v}, {x, y} (Element) (CapitalOmega));
GetInactivePDE(state) - pde
(* {0, 0} *)

Link to notebook with this question and code