Solid Mechanics

PDE Modeling

Oliver  Rübenkönig 
Wolfram Research

Previously on PDE modeling

Acoustic & HeatTransfer in 12.0

Mass Transport added in 12.1

Framework & Basic Building Blocks in 12.2

[Electromagnetics Models 13.0]

Experimental SolidMechanics 13.0

PDEModels Overview

paclet:PDEModels/tutorial/PDEModelsOverview

Field of physics - SolidMechanics

xyzPDEComponent[vars, pars] -> PDE

xyzCondition[predicate, vars, pars] -> DirchletCondition

xyzValue[predicate, vars, pars] -> NeumannValue

Overview example

Geometry

Ω=Cuboid[{0,0,0},{10,1,1}];

Variables and parameters

vars={{u[x,y,z],v[x,y,z],w[x,y,z]},{x,y,z}};
pars=<|"Material"->(Entity["titanium", Entity["Element","Titanium"], "Entity[\"Element\", \"Titanium\"]", "element"])|>;

<|"YoungModulus"->\[FormalCapitalY],"PoissonRatio"->\[FormalNu]|>;

PDE component:

op=SolidMechanicsPDEComponent[vars,pars];

MatrixForm[op]

Boundary conditions:

(Subscript[Γ, force])=SolidBoundaryLoadValue[x==10,vars,pars,<|"Force"->{0,0,(Quantity[-1000, "N", newtons, "Newtons"])}|>]

(Subscript[Γ, wall])=SolidFixedCondition[x==0,vars,pars]

Solid Mechanics

beamDisplacement=NDSolveValue[{op==(Subscript[Γ, force]),(Subscript[Γ, wall])},{u[x,y,z],v[x,y,z],w[x,y,z]},{x,y,z}∈Ω];

Head/@Head/@beamDisplacement

{(Subscript[p, x]),(Subscript[p, y]),(Subscript[p, z])} -> {(Subscript[p, x])+u(x,y,z),(Subscript[p, y])+v(x,y,z),(Subscript[p, z])+w(x,y,z)}

Visualize

Exaggeration of the deformation:

VectorDisplacementPlot3D[beamDisplacement,{x,y,z}∈Ω,]

Stress & Strain

SolidMechanicsPDEComponent[vars, pars]
->displacement

SolidMechanicsStrain[vars, pars,displacement]
->strain

SolidMechanicsStress[vars, pars,strain]
->stress

Strain

strains=SolidMechanicsStrain[vars,pars,beamDisplacement]

(({{(Subscript[ϵ, x,x]), (Subscript[γ, x,y]), (Subscript[γ, x,z])}, {(Subscript[γ, y,x]), (Subscript[ϵ, y,y]), (Subscript[γ, y,z])}, {(Subscript[γ, z,x]), (Subscript[γ, y,z]), (Subscript[ϵ, z,z])}}))

normal strains ϵ, shear strains γ, engineering strains (Subscript[γ, i,j])=2(Subscript[ϵ, i,j])

Show[,
SliceContourPlot3D[strains[[1,1]],mesh,{x,y,z}∈mesh,]]

Stress

stresses=SolidMechanicsStress[vars,pars,strains]

(({{(Subscript[σ, x,x]), (Subscript[τ, x,y]), (Subscript[τ, x,z])}, {(Subscript[τ, y,x]), (Subscript[σ, y,y]), (Subscript[τ, y,z])}, {(Subscript[τ, z,x]), (Subscript[τ, y,z]), (Subscript[σ, z,z])}}))

normal stress σ, shear stress τ

Show[,
SliceContourPlot3D[stresses[[1,2]],mesh,{x,y,z}∈mesh,]]

Solid Mechanics + Heat = Thermoelasticity

Show[
Graphics[{Text["Aluminum",{len 4/5,3hi}],Text["Carbon steel",{len 4/5,-hi}],Text["100C",{-0.01,hi}]}]
,Graphics[Line[{{0,-hi},{0,3 hi}}]],Graphics[Table[Line[{{0,y},{-0.0025,y-0.0025}}],{y,3 hi,-hi,-hi/2}]],mesh["Wireframe"["MeshElementStyle"->{Directive[EdgeForm[Green],FaceForm[Green]],Directive[EdgeForm[Red],FaceForm[Red]]}]]
,mesh["Wireframe"["MeshElement"->"BoundaryElements"]]
]

Mesh generation

len=1/10;hi=5/1000;
bmesh=ToBoundaryMesh[
"Coordinates"->{{0,0},{len,0},{len,hi},{len,2 hi},{0, 2hi},{0,hi}},"BoundaryElements"->{LineElement[{{1,2},{2,3},{3,4},{4,5},{5,6},{6,1},{3,6}}]}];
material=<|"Al"->1,"Ti"->2|>;
mesh=ToElementMesh[bmesh,"RegionMarker"->{{{len/2,3 hi/2},material["Al"]},{{len/2,hi/2},material["Ti"]}}];
mesh["Wireframe"["MeshElementStyle"->{Directive[FaceForm[Green]],Directive[FaceForm[Red]]}]]

Variables

svars={{u[x,y],v[x,y]},{x,y}};
tvars={T[x,y],{x,y}};

Material parameters

smtpars=<|"Material"->{
{y>=hi,(Entity["aluminum", Entity["Element","Aluminum"], "Entity[\"Element\", \"Aluminum\"]", "element"])},
{y<hi,(Entity["titanium", Entity["Element","Titanium"], "Entity[\"Element\", \"Titanium\"]", "element"])}},
"ThermalStrainTemperature"->T[x,y],
"Thickness"->10^-3
|>;

PDE model

pdeModel={
SolidMechanicsPDEComponent[svars,smtpars]=={0,0},
SolidFixedCondition[x==0,svars,smtpars],
HeatTransferPDEComponent[tvars,smtpars]==HeatTransferValue[y==0||y==2 hi,tvars,smtpars,<|"HeatTransferCoefficient"->(Quantity[5, "W""/("(("m")^2)"K"")", StringTemplate[`1` per `2`][watts, StringTemplate[`1` squared][meter] kelvin], (("Watts")/("Kelvins" (("Meters")^2)))]),"AmbientTemperature"->(Quantity[15, "°C", degrees Celsius, "DegreesCelsius"])|>],
HeatTemperatureCondition[x==0,tvars,smtpars,<|"SurfaceTemperature"->(Quantity[100, "°C", degrees Celsius, "DegreesCelsius"])|>]
};

Solve

AbsoluteTiming[{uif,vif,Tif}=NDSolveValue[pdeModel, {u[x,y],v[x,y],T[x,y]},{x,y}∈mesh];]

Post-process

VectorDisplacementPlot[{uif,vif},{x,y}∈Rectangle[{0,0},{len,2hi}]]

ContourPlot[Tif,{x,y}∈mesh,]

Strain

strain=SolidMechanicsStrain[svars,smtpars,{uif,vif}];

ContourPlot[strain[[1,1]],{x,y}∈mesh,]

Stress

newPars=smtpars;
newPars["ThermalStrainTemperature"]=Tif;
stress=SolidMechanicsStress[svars,newPars,strain];

vms=PDEModels`VonMisesStress[svars,smtpars,stress];

ContourPlot[vms,{x,y}∈mesh,]

Eigenmodes

Load a STEP file:

shape=OpenCascadeShapeImport[FileNameJoin[{,"linkrods.step"}]];

Convert to a boundary element mesh:

bmesh=OpenCascadeShapeSurfaceMeshToBoundaryMesh[shape];
bmesh["Wireframe"[]]

Generate a mesh:

mesh=ToElementMesh[bmesh,"MeshOrder"->1];

Setup the solid mechanics PDE:

emvars={{u[t,x,y,z],v[t,x,y,z],w[t,x,y,z]},t,{x,y,z}};
empars=<|"Material"->(Entity["titanium", Entity["Element","Titanium"], "Entity[\"Element\", \"Titanium\"]", "element"]),"AnalysisType"->"Eigenmode"|>;
emop=SolidMechanicsPDEComponent[emvars,empars];

Compute the eigenvalues and eigenvectors:

{evals,evecs}=NDEigensystem[{emop=={0,0,0}},emvars[[1]],t,{x,y,z}∈mesh,12];

Compute the eigenfrequencies:

eigenfrequencies=Sqrt[evals]/(2π)

Visualize the deformation modes at their frequencies:

GraphicsGrid[Partition[MapThread[Show[
mesh["Edgeframe"[PlotLabel->Style["freq="<>ToString[#1]<>" Hz",8]]],
deformedMesh2=ElementMeshDeformation[mesh,#2,"ScalingFactor"->35];
deformedMesh2["Edgeframe"],
deformedMesh2["Wireframe"[
"ElementMeshDirective"->Directive[EdgeForm[],FaceForm[LightGray]]]]
]&,{eigenfrequencies[[7;;]],evecs[[7;;]]}],3]]

Nonlinear Elastic Material

Region and vars:

Ω=Rectangle[{0,0},{1,1}];
vars={{u[x,y],v[x,y]},{x,y}};

Hypoelastic (not Hyperelastic!) material model as a function:

materialModel[vars_,pars_,data_]:=Module[
{n,e0,s0,ν,dim,idm,ea,evol,eTemp,ee,se,K,stressMatrix,strainMeasure},
(* model parameters *)
n=10;e0=1/1000;s0=10;ν=3/10;
strainMeasure="StrainMeasure"/.data;
(*model*)
dim=Length[vars[[1]]];
idm=IdentityMatrix[dim];
ea=strainMeasure+$MachineEpsilon;
evol=Tr[ea];
eTemp=(ea-idm*evol/3);
ee=Sqrt[2/3 Total[eTemp^2,2]];
se=Piecewise[{
{s0*(Sqrt[(1+n^2)/(n-1)^2-(n/(n-1)-ee/e0)^2]-1/(n-1)),ee<=e0}},s0*((ee/e0)^(1/n))];
K=(n*s0/e0)/(3*(1-2 ν));
stressMatrix=K*evol*idm+2/3*se*eTemp/ee;
stressMatrix
]

Parameters:

pars=<|"MaterialModelFunction"->materialModel,"ModelForm"->"PlaneStrain","Thickness"->1|>;

PDE:

pde={SolidMechanicsPDEComponent[vars,pars]==SolidBoundaryLoadValue[x==1,vars,pars,<|"Pressure"->{p,0}|>]
,DirichletCondition[{u[x,y]==0},x==0]
,DirichletCondition[{v[x,y]==0},y==0]
};

Parametric in the pressure:

pfun=ParametricNDSolveValue[pde,vars[[1,All,0]],{x,y}∈Ω,p];

Pressure 1-20:

data=Table[{pfun[p][[1]][1,0],p},{p,1,20,1}];

Force vs. Displacement:

ListLinePlot[data,]

A few more things

Analysis types

Stationary

Transient

Eigenmode

Frequency Response

Parametric

Material Models

Linear elastic isotropic/orthotropic/anisotropic

?Nonlinear elastic (Hypoelastic)

Large deformation (e.g. Hyperelastic) -> future version

Geometric models

OpenCascadeLink:

STEP import

Shout out

Mauricio Fernández

Craig Carter

Thank you for your attention!

paclet:guide/PDEModeling

paclet:PDEModels/tutorial/PDEModelsOverview

ruebenko@wolfram.com

Initializations

Needs["NDSolve`FEM`"]
Needs["OpenCascadeLink`"]