You are using a browser not supported by the Wolfram Cloud
Supported browsers include recent versions of Chrome, Edge, Firefox and Safari.
I understand and wish to continue anyway »
WOLFRAM NOTEBOOK
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
In[]:=
stress=makeStressModel[e1];pde=makePDEModel[stress/.materialData];{uif,vif}=NDSolveValue[{pde==bc2,bc1},{u,v},{x,y}∈Ω];mesh=uif["ElementMesh"];Show[ToBoundaryMesh[mesh]["Wireframe"],deform2=ElementMeshDeformation[mesh,{uif,vif},"ScalingFactor"->1]["Wireframe"["MeshElementStyle"EdgeForm[Red]]]]
Out[]=
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.