What’s New in PDE Modeling

Oliver Rübenkönig, WRI

#WolframTechConf

PDEModels Overview

This is the most important slide

Where to find examples…

PDEModels Overview

... if your application is missing

PDEModels Overview V14.3 - V15.0

PDETerms

Heat Transfer

Solid Mechanics

Electromagnetics

Fluid Dynamics

LinearSolve

ExtrapolationMethod

How do you like your eggs

Introduction to PDE modeling

Visit Ricardo’s talk: Thursday, November 6

PDETerms

Coordinate Charts

Coordinate Charts

op=DiffusionPDETerm[{u[r,z],{r,z}},1,<|"CoordinateChart"->"Polar"|>]

Activate[op]

-Laplacian[u[r,z],{r,z},"Polar"]

Heat Transfer

Thickness

Makes it easier to model in 2D

Solid Mechanics

Plasticity

EquivalentStrain

Spring/Damper boundary conditions

NormalForce / NormalPressure

Plasticity

New Plasticity Monograph

Visit Alessandro’s talk: Thursday, November 6

EquivalentStrain -  “VonMisesStress” for strain

(Subscript[ϵ, eqv])=(Sqrt[(1/2)(((((Subscript[ϵ, xx])-(Subscript[ϵ, yy])))^2)+((((Subscript[ϵ, yy])-(Subscript[ϵ, zz])))^2)+((((Subscript[ϵ, zz])-(Subscript[ϵ, xx])))^2)+6((Subscript[γ, xy]^2)+(Subscript[γ, xz]^2)+(Subscript[γ, yz]^2)))])

EquivalentStrain -  Example

beam=;

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"])|>;

displacement=NDSolveValue[{SolidMechanicsPDEComponent[vars,pars]==SolidBoundaryLoadValue[x==10,vars,pars,<|"Force"->{0,0,-2*(10^7)}|>],SolidFixedCondition[x==0,vars,pars]}
,{u[x,y,z],v[x,y,z],w[x,y,z]},{x,y,z}∈beam];

strains=SolidMechanicsStrain[vars,pars,displacement]

equivalentStrain=PDEModels`EquivalentStrain[vars,pars,strains]

VectorDisplacementPlot3D[{displacement,equivalentStrain},{x,y,z}∈beam,VectorSizes->2,PlotPoints->40]

Spring/Damper boundary conditions

Ω=Rectangle[{0,0},{0.5,0.01}];
vars={{u[x,y],v[x,y]},{x,y}};
pars=<|"YoungModulus"->200*10^9,"PoissonRatio"->1/3,"Thickness"->1,"SolidMechanicsModelForm"->"PlaneStrain"|>;

op=SolidMechanicsPDEComponent[vars,pars];
fixed=SolidFixedCondition[0<=x<=0.1&&y==0,vars,pars];
load=SolidBoundaryLoadValue[0.15<=x<=0.35&&y==0.01,vars,pars,<|"Pressure"->{0,-10^4}|>];

displacement=NDSolveValue[{op==load,fixed},vars[[1]],{x,y}∈Ω];

VectorDisplacementPlot[displacement,{x,y}∈Ω,]

springBC=SolidBoundaryLoadValue[0.4<=x<=0.5&&y==0,vars,pars,<|"SpringConstant"->Quantity[10,"Kilonewtons"/"Millimeters"]|>];

displacement=NDSolveValue[{op==load+springBC,fixed},vars[[1]],{x,y}∈Ω];

VectorDisplacementPlot[displacement,{x,y}∈Ω,]

Electromagnetics

MagneticFluxDensity - B-Field

MagneticFieldIntensity - H-Field

Work for both scalar and vector potentials

MagneticFluxDensity

mesh=;

Show[HighlightMesh[mesh,],PlotRange->{{-0.08,0.08},{-0.15,0.15}}]

Magnet transversely magnetized in the direction of the y axis:

vars={Vm[x,y],{x,y}};
pars=<|"Magnetization"->{0,Piecewise[{{400000,RegionMember[Rectangle[{-0.05,-0.1},{0.05,0.1}]][{x,y}]}},0]}|>;

VmFun=NDSolveValue[{MagnetostaticPDEComponent[vars,pars]==0,MagneticPotentialCondition[x==-1||x==1||y==-1||y==1,vars,pars]},Vm[x,y],{x,y}∈mesh]

fluxDensity=PDEModels`MagneticFluxDensity[vars,pars,VmFun]

Show[Graphics[{Gray,Rectangle[{-0.05,-0.1},{0.05,0.1}]}],StreamPlot[fluxDensity,{x,-0.3,0.3},{y,-0.3,0.3},StreamColorFunction->"Rainbow"]]

FluidFlow

FluidViscousity

FluidViscousStress

Verification Tests

Swirl Flow

Martin’s talk: Wednesday, November 5

Swirl Flow

Axisymmetric -> ∂/∂θ=0 and by setting (Subscript[v, ϕ])=0

Swirl Flow -> ∂/∂θ=0

Swirl Flow

bw=0.02;bh=0.05;
radius=0.001;height=0.015;
beaker=Rectangle[{0,0},{bw,bh}];
shaft=Rectangle[{0,height},{radius,bh}];
rotorThickness=2radius;
rotor1=Rectangle[{0,height},{1/2*bw,height+rotorThickness}];
rotor2=Rectangle[{0,height+2*rotorThickness},{1/3*bw,height+3rotorThickness}];
rotor3=Rectangle[{0,height+4*rotorThickness},{1/2*bw,height+5rotorThickness}];
mixer=RegionUnion[shaft,rotor1,rotor2,rotor3];
Ω=RegionDifference[beaker,mixer];
Show[
RegionPlot[Ω,FrameTicks->None],
HighlightMesh[mixer,{Style[1,Red],Style[2,LightGray]}],AspectRatio->Automatic,ImageSize->Small]

mesh=ToElementMesh[Ω,"MaxCellMeasure"->0.0000005]

Swirl Flow

Velocity in the θ direction is not zero, and the dependent variable is specified.

Specify the variables:

vars={{u[r,z],v[r,z],w[r,z],p[r,z]},{r,θ,z}};

Set the material properties and the coordinate chart.

pars=<|"DynamicViscosity"->10^-3,"MassDensity"->10^3,"CoordinateChart"->"Cylindrical"|>;

Create the operator:

op=FluidFlowPDEComponent[vars,pars];

Swirl Flow

Parametric boundary conditions with the angular velocity ω:

mixerBC=DirichletCondition[{u[r,z]==0,v[r,z]==r*ω,w[r,z]==0},height<=z<bh&&radius<=r<bw];

Set up no-slip boundary conditions at the beaker wall:

wall=DirichletCondition[{u[r,z]==0,v[r,z]==0,w[r,z]==0},r==bw||z==0];

On the axis of symmetry, the u and v velocity are set to 0:

symmetryAxis=DirichletCondition[{u[r,z]==0,v[r,z]==0},r==0];

The top of the beaker is open, and there is no fluid velocity in the z direction:

top=DirichletCondition[{w[r,z]==0},z==bh];

A reference pressure condition is set at the top-right corner:

refPressure=DirichletCondition[{p[r,z]==0},r==bw&&z==bh];

Swirl Flow

The parametric function is created:

fun=ParametricNDSolveValue[{Activate[op]=={0,0,0,0},mixerBC,wall,symmetryAxis,top,refPressure},{u,v,w,p},{r,z}∈mesh,ω,Method->{"PDEDiscretization"->{"FiniteElement"
,"InterpolationOrder"->{u->2,v->2,w->2,p->1}
}}]

Compute the solution for an angular velocity of π [(QuantityTF[, "m""/""s", meters per second, (("Meters")/("Seconds"))])]:

{ufun,vfun,wfun,pfun}=fun[π]

Swirl Flow

Compute the magnitude of the velocity:

velMag=FunctionInterpolation[Sqrt[ufun[r,z]^2+vfun[r,z]^2+wfun[r,z]^2],{r,z}∈mesh]

Visualize the magnitude of the velocity and contours of the v component of the velocity:

Show[
HighlightMesh[mixer,{Style[1,Red],Style[2,LightGray]}],
DensityPlot[velMag[r,z],{r,z}∈mesh,],ContourPlot[vfun[r,z],{r,z}∈mesh,PlotRange->All,Contours->15,ContourShading->None],AspectRatio->Automatic,ImageSize->Small]

Swirl Flow

{minMag,maxMag}=MinMax[velMag["ValuesOnGrid"]];
legendBar=BarLegend[{"TemperatureMap",{minMag,maxMag}},];
frames=Table[{ufun,vfun,wfun,pfun}=fun[i*π];velMag=EvaluateOnElementMesh[{r,z},Sqrt[ufun[r,z]^2+vfun[r,z]^2+wfun[r,z]^2],mesh];Show[HighlightMesh[mixer,{Style[1,Red],Style[2,LightGray]}],Legended[DensityPlot[velMag[r,z],{r,z}∈mesh,],legendBar],],{i,0.1,1,0.1}];
frames=Rasterize[#,"Image",ImageResolution->120]&/@frames;
ListAnimate[frames]

LinearSolve

Method->”MUMPS”

“MUMPS”

Real & Complex valued machine precision

sa=SparseArray[{{1.2,2.3,0.0},{0.0,4.5,0.0},{0.0,2.1,6.5}}];
b={8.9,7.1,3.2};
r2=LinearSolve[sa,b,Method->"Mumps"];
Norm[Dot[sa.r2]-b]

“MUMPS” - ARM Hardware

V15.0: MUMPS default FEM solver on ARM hardware

V14.3

NDSolveValue[{SolidMechanicsPDEComponent[vars,pars]==SolidBoundaryLoadValue[x==1,vars,pars,<|"Pressure"->{9*10^8,0}|>],SolidFixedCondition[x==0,vars,pars]},{u[x,y],v[x,y]},{x,y}∈Ω
,Method->{"PDEDiscretization"->{"FiniteElement","PDESolveOptions"->{"LinearSolver"->"MUMPS"}}}
];

V15.0

NDSolveValue[{SolidMechanicsPDEComponent[vars,pars]==SolidBoundaryLoadValue[x==1,vars,pars,<|"Pressure"->{9*10^8,0}|>],SolidFixedCondition[x==0,vars,pars]},{u[x,y],v[x,y]},{x,y}∈Ω
];

ExtrapolationMethod

ifun=Interpolation[{1,2,3,5,8,5}]

ifun[0]

ifun=Interpolation[{1,2,3,5,8,5},ExtrapolationMethod->Indeterminate]

ifun[0]

New Application Examples

Schrödinger-Poisson

Schrödinger-Newton

Cooling a Room with HVAC

Schrödinger-Newton

Spine Surgery Rod

Slope Stability

Biomechanics of Human Tendon

3D Print and Mechanical Design

Thanks

Watch out for the yearly FEM/PDE conference

Init

Needs["NDSolve`FEM`"]