#### Large Deformation

Large Deformation

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.

Out[]=

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: https://www.comsol.com/model/download/837521/models.sme.large_deformation_beam.pdf. There is also a video demonstration at https://www.comsol.com/video/modeling-nonlinear-structural-mechanics-in-comsol-multiphysics-aug-16-2018 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 we are talking about the lower of the two beams.

Set up some parameters.

In[]:=

Needs["NDSolve`FEM`"]length=3.2;hight=1/10;thickness=1/10;Ω=Rectangle[{0,0},{length,hight}];mesh=ToElementMesh[Ω];bmesh=ToBoundaryMesh[mesh];

In[]:=

(*Graphics[{{Gray,Line[{{0,-0.04},{0,0.14}}],Table[Line[{{-0.1,-0.02+j},{0,-0.02+0.02+j}}],{j,0,0.125,0.025}]},{Directive[EdgeForm[Black],FaceForm[]],Ω},{Blue,Arrow[{{length,5*hight},{length,hight}}],Text["Fy",{length,5*hight},{-2,0}]},{Blue,Arrow[{{length*1.3,1/2hight},{length,1/2hight}}],Text["Fx",{length*1.3,1/2*hight},{0,-2}]}}]*)

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.

In[]:=

(*forceperarea*)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.

In[]:=

PlaneStress[{Y_,nu_},{u_,v_},X:{x_,y_}]:=Module[{pStress},(*usesglobalthickness*)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}∈Ω];Show[bmesh["Wireframe"],deform2=ElementMeshDeformation[mesh,{uif2,vif2},"ScalingFactor"->1]["Wireframe"["MeshElementStyle"EdgeForm[Red]]]]

Out[]=

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

σ xx |

σ yy |

τ xy |

1 | nu | 0 |

nu | 1 | 0 |

0 | 0 | (1-nu)/2 |

ϵ xx |

ϵ yy |

γ xy |

This is then converted to a 2 by 2 matrix

In[]:=

Normal[SymmetrizedArray[Thread[Rule[{{1,1},{2,2},{1,2}},{,,}]],Automatic,Symmetric]]

σ

xx

σ

yy

τ

xy

Out[]=

{{,},{,}}

σ

xx

τ

xy

τ

xy

σ

yy

And plugged into (with an additional thickness factor):

In[]:=

Table[Inactive[Div][-{{,},{,}}[[i,;;]],X],{i,Length[U]}]

σ

xx

τ

xy

τ

xy

σ

yy

Out[]=

{·{-,-},·{-,-}}

∇

{x,y}

σ

xx

τ

xy

∇

{x,y}

τ

xy

σ

yy

In[]:=

makeStressModel[e_]:=Block[{materialData,materialModel,exx,eyy,gxy,eVec,strain,nu,Y,stress},(*linearelasticplanestressmodel*)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]]);(*forplanestressezz=nu/(1-nu)*(exx+eyy)doesnotplayarolesincenu=0*)strain={exx,eyy,gxy};stress=(materialModel.strain);Normal[SymmetrizedArray[Thread[Rule[{{1,1},{2,2},{1,2}},stress]],Automatic,Symmetric]]]

In[]:=

(*usesglobaldefinedthickness*)makePDEModel[stressMatrix_]:=Table[Inactive[Div][-thickness*stressMatrix[[i,;;]],X],{i,Length[U]}]

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]};(*geometriclinear*)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}∈mesh];Show[bmesh["Wireframe"],deform1=ElementMeshDeformation[mesh,{uif,vif},"ScalingFactor"->1]["Wireframe"["MeshElementStyle"EdgeForm[Brown]]]]

Out[]=

In[]:=

Norm[uif["ValuesOnGrid"]-uif2["ValuesOnGrid"]]Norm[vif["ValuesOnGrid"]-vif2["ValuesOnGrid"]]

1.43378×6.49507×

-10

10

-9

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.

In[]:=

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

With the same approach as above

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

Question: Why is this approach wrong or where did I make a mistake?

A second way to specify the large strains is by

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

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

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

Still, no joy.