Modeling 
with 
Partial Differential Equations

Oliver   Rübenkönig 
Wolfram Research

(RowDefault[, (Spacer1), , (Spacer1), ])  #WolframTechConf

PDE Models

paclet:PDEModels/tutorial/PDEModelsOverview

Previously on PDE modeling

Acoustics & Heat Transfer in 12.0

Mass Transport in 12.1

Framework & Basic Building Blocks in 12.2

[Electromagnetics Models 13.0]

Experimental SolidMechanics 13.0

RegionSymmetry & more SolidMechancis 13.1

More RegionSymmetry 13.2

Essentials

(InactiveDiv[(κ.(InactiveGrad[u[x], {x}, Inactive])), {x}, Inactive])

Variables

Always: Dependent (u)

Always: Independent (x)

Possibly: Time/Frequency (t/ω)

Parameters

Often: Material (κ)

Field of physics

xyzPDEComponent[vars, pars] -> PDE

xyzCondition[predicate, vars, pars] -> DirchletCondition

xyzValue[predicate, vars, pars] -> NeumannValue

Overview example

({{, ({{Heat equation}, {Heat Flux}, {Fixed Temperature}})}})

Geometry

Ω3D=RegionDifference[Cylinder[{{0,0,0},{0,0,0.14}},0.1],Cylinder[{{0,0,0},{0,0,0.14}},0.02]];

Mesh

mesh=ToElementMesh[Ω3D];
mesh["Wireframe"]

Variables and parameters

vars3D={T[x,y,z],{x,y,z}};
(*pars3D=<|"Material"->(Entity["aluminum", Entity["Element","Aluminum"], "Entity[\"Element\", \"Aluminum\"]", "element"])|>;*)
pars3D=<|"ThermalConductivity"->(Quantity[235., "W""/(""m" "K"")", watts per meter kelvin, (("Watts")/("Kelvins" "Meters"))])|>;

PDE component

op=HeatTransferPDEComponent[vars3D,pars3D]

Boundary Condition 1

heatFlux=<|"HeatFlux"->(Quantity[500000, "W""/"(("m")^2), watts per meter squared, (("Watts")/(("Meters")^2))])|>;

(Subscript[Γ, flux])=HeatFluxValue[Sqrt[x^2+y^2]<=0.02 ,vars3D,pars3D,heatFlux]

Boundary Condition 2

refTemp=<|"SurfaceTemperature"->(Quantity[273.15, "K", kelvins, "Kelvins"])|>;

(Subscript[Γ, temp])=HeatTemperatureCondition[z==0||z==0.14,vars3D,pars3D,refTemp]

Solve

AbsoluteTiming[
solution3D=NDSolveValue[{op==(Subscript[Γ, flux]),(Subscript[Γ, temp])},T,{x,y,z}∈mesh];]

Visualize

Show[mesh["Edgeframe"],
SliceContourPlot3D[solution3D[x,y,z],{x,y,z}∈mesh,ColorFunction->"TemperatureMap",Boxed->False,Axes->False,PlotLegends->BarLegend[Automatic,25,LegendLabel->"T [K]"]]]

New in Version 13.1

RegionSymmetry

Geometry

Ω=Rectangle[{0.02,0},{0.1,0.14}];

Variables

vars={T[r,z],{r,z}};
pars=Join[pars3D,<|"RegionSymmetry"-> "Axisymmetric"|>]

PDE component

op=HeatTransferPDEComponent[vars,pars]

Div[Grad[-235 T[r,z],{r,θ,z},"Cylindrical"],{r,θ,z},"Cylindrical"]

Activate[op]-%

Boundary   Conditions

bc1=HeatFluxValue[r==0.02,vars,pars,heatFlux]

bc2=HeatTemperatureCondition[z==0||z==0.14,vars,pars,refTemp]

Solve

AbsoluteTiming[solution=NDSolveValue[{op==bc1,bc2},T,{r,z}∈Ω];]

Visualize

GraphicsRow[{ContourPlot[solution[r,z],{r,z}∈Ω,PlotRange->All,ColorFunction->"TemperatureMap",PlotLegends->BarLegend[Automatic,25,LegendLabel->"T [K]"],PlotLabel->"Axisymmetric",ImageSize->Small],
ContourPlot[solution3D[r,0,z],{r,z}∈Ω,PlotRange->All,ColorFunction->"TemperatureMap",PlotLegends->BarLegend[Automatic,25,LegendLabel->"T [K]"],PlotLabel->"3D",ImageSize->Small]}]

Visualize in 3D

RegionPlot3D[(0.02^2)≤(x^2)+(y^2)≤(0.1^2)&&0≤PlanarAngle[{0,0}->{{0.02,0},{x,y}}]≤((4 π)/3),{x,-0.1,0.1},{y,-0.1,0.1},{z,0,0.14},]

RegionSymmetry

PDETerms

DiffusionPDETerm, ConservativeConvectionPDETerm, DerivativePDETerm

13.2 LaplacianPDETerm

PDEComponents

HeatTransferPDEComponent, MassTransportPDEComponent

13.2 AcousticPDEComponent, PoissonPDEComponent, HelmholtzPDEComponent, WavePDEComponent

13.3 SolidMechanicsPDEComponent

PDE Models

Disk Brake, Tubular Reactor, Spherical Capacitor, Solute Centrifugation

13.2 Helmholtz Resonator

SolidMechanics

Hyperelastic Material

St. Venant-Kirchhoff

Neo-Hookean

Rayleigh Damping

Hyperelastic material

Extremely large reversible deformation

Hyperelastic material

Set up

vars={{u[x,y],v[x,y]},{x,y}};
pars=<|"MaterialModel"->"NeoHookeanIsotropic","LameParameter"->10^6,"ShearModulus"->5000,"ModelForm"->"PlaneStress","Thickness"->0.01|>;

Set up

neoHooke[f_]:={SolidMechanicsPDEComponent[vars,pars]==SolidBoundaryLoadValue[x==1,vars,pars,<|"Force"->{Quantity[f,"Newtons"],0}|>],SolidFixedCondition[x==0,vars,pars]}

Side note

SolidMechanicsPDEComponent[vars,pars]

Solve PDE

displacement2D=NDSolveValue[
neoHooke[300],{u[x,y],v[x,y]},{x,y}∈Rectangle[{0,0},{1,1}]];

Visualize

VectorDisplacementPlot[displacement2D,{x,y}∈Rectangle[{0,0},{1,1}],]

450N - Ops

displacement2D=NDSolveValue[
neoHooke[450],{u[x,y],v[x,y]},{x,y}∈Rectangle[{0,0},{1,1}]];

Solver restart - ParametricNDSolve

pfunNeoHooke=ParametricNDSolveValue[
neoHooke[f],{u[x,y],v[x,y]},{x,y}∈Rectangle[{0,0},{1,1}],f];

Solver restart

fMax=500;
nsteps=20;
displacementNeoHooke={0,0};

Monitor[Do[
fNew=step*fMax/nsteps;
displacementNeoHooke=pfunNeoHooke[fNew,
"InitialSeeding"->{u[x,y]==displacementNeoHooke[[1]],v[x,y]==displacementNeoHooke[[2]]}];
,{step,1,nsteps,1}],
step];

Visualization

VectorDisplacementPlot[displacementNeoHooke,{x,y}∈Rectangle[{0,0},{1,1}],]

Rayleigh Damping

M(Overscript[u, ¨])+Ku=F

M(Overscript[u, ¨])+C(OverDot[u])+Ku=F

C=α M+β K

Set up

vars={{u[t,x,y],v[t,x,y]},t,{x,y}};
pars=<|"Material"->(Entity["copper", Entity["Element","Copper"], "Entity[\"Element\", \"Copper\"]", "element"]),"MassDamping"->40,"StiffnessDamping"->0.0004,"Thickness"->0.001,"ModelForm"->"PlaneStress"|>;

PDE

pde={SolidMechanicsPDEComponent[vars,pars]==SolidBoundaryLoadValue[x>=1/10,vars,pars,<|"Force"->{0,(Quantity[-20, "N", newtons, "Newtons"])}|>],
SolidFixedCondition[x<=-1/10,vars,pars],{u[0,x,y]==0,v[0,x,y]==0,Derivative[1,0,0][u][0,x,y]==0,Derivative[1,0,0][v][0,x,y]==0}};

Solve

AbsoluteTiming[Monitor[displacement=NDSolveValue[pde,{u,v},{x,y}∈,{t,0,0.05},EvaluationMonitor:>(monitor=Row[{"t = ",CForm[t]}])]
,monitor];]

Visualize

Plot[Evaluate[displacement[[2]][t,0.17,0.01]],{t,0,0.05}]

Visualize

Needs["NDSolve`FEM`"]
mesh=displacement[[1]]["ElementMesh"];
frames=Table[Show[ToBoundaryMesh[mesh]["Wireframe"],ElementMeshDeformation[{t,mesh},displacement,"ScalingFactor"->1]["Wireframe"[
"ElementMeshDirective"->Directive[EdgeForm[],FaceForm[LightOrange]]]],Boxed->True,BoxStyle->Transparent],
{t,0,0.05,0.05/25}];
frames=Rasterize[#1,"Image",ImageSize->Medium]&/@frames;
ListAnimate[frames]

Summary

PDEModels Overview

RegionSymmetry

Hyperelastic Material

Rayleigh Damping

Thanks

Humphrey Hardy

Christopher Purcell

Stan Connell

Harald Biller

Questions?

Initialization

Needs["NDSolve`FEM`"]

m=NDSolve`FEM`ToElementMesh[RegionDifference[Cylinder[{{0,0,0},{0,0,0.14}},0.1],Cylinder[{{0,0,0},{0,0,0.14}},0.02]]];
Show[
m["Edgeframe"],
m["Wireframe"["MeshElementStyle"->Directive[Opacity[0.1],EdgeForm[],FaceForm[Gray]]]],
m["Wireframe"[ElementMarker==2||ElementMarker==3,"MeshElementStyle"->Directive[Opacity[0.4],EdgeForm[],FaceForm[Blue]]]],
m["Wireframe"[ElementMarker==4,"MeshElementStyle"->Directive[Opacity[0.4],EdgeForm[],FaceForm[Red]]]]
]

Graphics[{
{{FaceForm[LightGray],EdgeForm[Black],Rectangle[{0,0},{1,1}]},Text["Ω",{0.9,0.1}]},
{Black,Line[{#,#-({1,+1}*0.025)}]&/@({0,#}&/@Range[-0.01,1.01,(1.01-(-0.01))/25]),Line[{{0,-0.01},{0,0.21}}]}
,{{Blue,Arrow[{{1,#},{1+1/8,#}}]&/@Range[0,1,1/9]},Text["f",{1+1/5,1/2}]}
}]

Needs["OpenCascadeLink`"]
Show[
Graphics3D[{
{Opacity[0.2],Cylinder[{{0,0,0},{0,0,1}}]}
,{Opacity[0.4],Polygon[{{0,0,0},{1/Sqrt[2],1/Sqrt[2],0},{1/Sqrt[2],1/Sqrt[2],1},{0,0,1}}]}
,{Arrowheads[Small],Arrow[{{0,0,0},{1/Sqrt[2],1/Sqrt[2],0}}],Text["r",{1/Sqrt[2],1/Sqrt[2],0}*0.75,{-1,1}]}
,{Arrowheads[Small],Arrow[BSplineCurve[{{1/2,0,0},{1/Sqrt[2],1/Sqrt[2],0}/2,{0,1/2,0}}]],Text["θ",{1/2,0,0},{-1,-1}]}
},Boxed->False]
,OpenCascadeAxis3D[]]