PDE Models 14.2

https://wolfram.com/language/core-areas/fem

Oliver Rübenkönig

PDEModels Overview

This is the most important slide

Where to find examples…

PDEModels Overview

... if your application is missing

PDEModels Overview V14.2

OpenCascadeLink (updated)

FEM

Solid Mechanics (updated)

Electromagnetics (new components)

OpenCascadeLink

1D/2D Drafting

Surgical Fixation Plate

FEM update - speed up

tEnd=π;
data=N[Table[{t,Sin[t]},{t,0,tEnd,tEnd/100}]];
ifun=Interpolation[data];

vars={u[t,x,y],t,{x,y}};
eqn=D[u[t,x,y],t]+DiffusionPDETerm[vars,1]==ifun[t];
RepeatedTiming[
NDSolveValue[{eqn,u[0,x,y]==0},u,{t,0,tEnd},{x,y}∈Disk[]];
]

OLD:

{3.50353, Null}

FEM -  NeumannBoundaryUnitNormal

NDSolve`FEM`BoundaryUnitNormal

Used when non-normal flux values can be specified

ElectricCurrentDensityValue[x>=0,{V[x,y],{x,y}},<|"CurrentDensity"->{J0x,J0y}|>]

FEM -  NeumannBoundaryUnitNormal

Model a different boundary integration term

op=D[u[t,x],t]+(InactiveDiv[(-1/20(InactiveGrad[u[t,x], {x}, Inactive])), {x}, Inactive]);

Default (RefLink[NeumannValue, paclet:ref/NeumannValue]): (OverVector[n])·(1/20 ∇u)

Would like: (OverVector[n])·(1/20 ∇u-u)

(RefLink[NeumannValue, paclet:ref/NeumannValue]) (OverVector[n])·{u} is set up:

NeumannValue[NeumannBoundaryUnitNormal[{x}].{u[t,x]},x==0];

FEM -  NeumannBoundaryUnitNormal

Poisson

solution=NDSolveValue[{-Laplacian[u[x,y],{x,y}]==1,DirichletCondition[u[x,y]==0,True]},u,{x,y}∈Disk[]];

Plot3D[solution[x,y],{x,y}∈Disk[]]

Flux through boundary

NIntegrate[NeumannBoundaryUnitNormal[{x,y}].Grad[-solution[x,y],{x,y}],{x,y}∈RegionBoundary[Disk[]]]

Flux through sub-region

NIntegrate[NeumannBoundaryUnitNormal[{x,y}].Grad[-solution[x,y],{x,y}],{x,y}∈RegionBoundary[Rectangle[{-1/2,-1/2},{1/2,1/2}]]]

Solid Mechanics

Alessandro’s Talk: Progress in Plasticity

Extended Plane Stress|Strain

Extended Plane Stress|Strain

Goal: Compute z-direction stress|strain for 2D cases

Plane Stress|Strain

Ω=Rectangle[{0,0},{40*10^-3,30*10^-3}];
pars=<|"YoungModulus"->200*10^9,"PoissonRatio"->0.3,"Thickness"->10^-4|>;

vars={{u[x,y],v[x,y]},{x,y}};

displacement=NDSolveValue[{SolidMechanicsPDEComponent[vars,pars]=={0,0},SolidFixedCondition[x==0,vars,pars],SolidDisplacementCondition[x==40*10^-3,vars,pars,<|"Displacement"->{-40/10^6,0}|>]},{u[x,y],v[x,y]},{x,y}∈Ω]

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

Plane Stress|Strain

strain=SolidMechanicsStrain[vars,pars,displacement]

Dimensions[strain]

stress=SolidMechanicsStress[vars,pars,strain]

Dimensions[stress]

estrain=PDEModels`SolidMechanicsExtendedStrain[vars,pars,displacement,strain,stress]

estrain[[3,3]]

ContourPlot[estrain[[3,3]],{x,y}∈Ω,PlotRange->All]

Extended Plane Stress|Strain

vars={{u[x,y],v[x,y],w[x,y]},{x,y,z}};

displacement=NDSolveValue[{SolidMechanicsPDEComponent[vars,pars]=={0,0},SolidFixedCondition[x==0,vars,pars],SolidDisplacementCondition[x==40*10^-3,vars,pars,<|"Displacement"->{-40/10^6,0}|>]},{u[x,y],v[x,y],w[x,y]},{x,y}∈Ω]

VectorDisplacementPlot[displacement[[1;;2]],{x,y}∈Ω]

Extended Plane Stress|Strain

strain=SolidMechanicsStrain[vars,pars,displacement]

Dimensions[strain]

ContourPlot[strain[[3,3]],{x,y}∈Ω,PlotRange->All,PlotLegends->Automatic]

stress=SolidMechanicsStress[vars,pars,strain]

stress[[3,3]]

Electromagnetics

Electromagnetics Overview - V14.1

Electrostatics - V14.0

Electric Currents - V14.1

Magnetostatics for Permanent Magnets - V14.2

Quasistatic Magnetic Fields - V14.2

Electromagnetics Overview

Maxwell’s equations

∇×E=-((∂B)/(∂t))

∇×H=J+((∂D )/(∂t))

∇·D=(Subscript[ρ, v])

∇·B=0

Solving full Maxwell’s equations is resource intensive

Simplified equations based on assumptions

Electromagnetics Overview

Magnetostatic

∇·((Subscript[μ, 0])(Subscript[μ, r])(-∇(Subscript[V, m])+M))=0

PDEComponent and Boundary Conditions

MagnetostaticPDEComponent

MagneticPotentialCondition

MagneticCurrentDensityValue

1D, 2D, 2D-axis and 3D models

Stationary models

Scalar equation - No currents

Permanent Magnent

bmesh=Import[]

outline=Show[bmesh["Edgeframe"],bmesh["Wireframe"["MeshElementStyle"->Directive[Opacity[0.2],FaceForm[LightBlue],EdgeForm[]]]],Axes->True,AxesLabel->{x,y,z}]

Mesh

regionMarker=<|"IronCoin"->1,"Iron"->2,"Magnet"->3,"Air"->4|>;
mesh=ToElementMesh[bmesh,"RegionMarker"->{{{-0.03,0.005,0.003},regionMarker["IronCoin"]},{{0.06,0.03,0.005},regionMarker["Iron"]},{{0.02,0.03,0.005},regionMarker["Magnet"]},{{0,0.05,0.05},regionMarker["Air"]}}]

parts=Map[mesh["Wireframe"[ElementMarker==#[[1]],"MeshElement"->"MeshElements","ElementMeshDirective"->Directive[EdgeForm[],FaceForm[#[[2]]]]]]&,{{1,Gray},{2,Gray},{3,Red},{4,LightBlue}}];
Show[mesh["Edgeframe"],parts[[1;;3]]]

Model

∇·((Subscript[μ, 0])(Subscript[μ, r])(-∇(Subscript[V, m])+M))=0

vars={Vm[x,y,z],{x,y,z}};

pars=<||>;

pars["RelativePermeability"]=Piecewise[{{5000,ElementMarker==regionMarker["IronCoin"]||ElementMarker==regionMarker["Iron"]}},1];

pars["Magnetization"]={Piecewise[{{800000,ElementMarker==regionMarker["Magnet"]}},0],0,0};

op=MagnetostaticPDEComponent[vars,pars];

Boundary Conditions

antisymmetryPlane=y==0;

Show[Region[Style[ImplicitRegion[antisymmetryPlane,{{x,-0.05,0.12},{y,0,0.1},{z,0,0.06}}],Red,Opacity[0.3]]],mesh["Edgeframe"]]

(Subscript[Γ, antisym])=MagneticPotentialCondition[antisymmetryPlane,vars,pars]

Solving

VmFun=NDSolveValue[{op==0,(Subscript[Γ, antisym])},Vm,{x,y,z}∈mesh]

DiscontinousInterpolatingFunction

DIFVmFun=DiscontinuousInterpolatingFunction[VmFun]

B-Field

Hfield=-Grad[DIFVmFun[x,y,z],{x,y,z}];

μr=pars["RelativePermeability"];

μ0=π/2500000;

M=pars["Magnetization"];

Bfield=EvaluateOnElementMesh[{x,y,z},μr*μ0 *(Hfield+M),mesh];
Bfield=Through[Bfield[x,y,z]];

mesh2=ElementMeshProjection[mesh,ReflectionTransform[{0,1,0}]];
mesh3=ElementMeshProjection[mesh,ReflectionTransform[{0,0,1}]];
mesh4=ElementMeshProjection[mesh3,ReflectionTransform[{0,1,0}]];
parts=Map[mesh["Wireframe"[ElementMarker==#[[1]],"MeshElement"->"MeshElements","ElementMeshDirective"->Directive[EdgeForm[],FaceForm[#[[2]]]]]]&,{{1,Gray},{2,Gray},{3,Red},{4,LightBlue}}];
Show[mesh["Edgeframe"],parts[[1;;3]],mesh2["Edgeframe"],mesh3["Edgeframe"],mesh4["Edgeframe"],SliceVectorPlot3D[Evaluate[Piecewise[{{{-1,1,1}*Bfield/.y->Abs@y,y<0}},Bfield/.y->Abs@y]], z==0,{x,-0.05,0.12},{y,-0.1,0.1},{z,0,0.06},VectorPoints->20,Boxed->False,Axes->False,VectorScaling->Automatic,VectorSizes->Large]]

Quasistatic Magnetics

σ((∂A)/(∂t))+∇×((1/(Subscript[μ, 0]))∇×A-M)= (Subscript[J, e])

PDEComponent and Boundary Conditions

MagneticPDEComponent

MagneticPotentialCondition

2D, 2D-axis models

[3D models]

Stationary models

Vector equation - Currents

Quasistatic Magnetics — Example

mesh["Wireframe"["MeshElementStyle"->FaceForm/@Join[markerColors,{LightOrange,Yellow}]]]

BH-Curve

HField=;
BField=;

ListLinePlot[Partition[Riffle[HField,BField],2],AxesLabel->{"H [A/m]","B [T]"}]

μData=BField/HField;

BHData=Transpose[{BField,μData}];

μB=Interpolation[BHData,InterpolationOrder->2
,ExtrapolationHandler->{If[#<BField[[1]],BHData[[1,2]],BHData[[-1,2]]]&,"WarningMessage"->False}]

Model

vars={{0,0,Az[x,y]},{x,y}};

BNorm=Sqrt[Total[Grad[Az[x,y],{x,y}]^2]+$MachineEpsilon];

μ0=π/2500000;
μr=Piecewise[{
{μB[BNorm]/μ0,ElementMarker==2||ElementMarker==3}
},1];

J0={0,0,Piecewise[{{10,ElementMarker==4},{-10,ElementMarker==6}},0]};

pars=<|"RelativePermeability"->μr,"ExternalCurrentSource"->J0|>;

Solving

AzSol=NDSolveValue[{MagneticPDEComponent[vars,pars]==0,MagneticPotentialCondition[x^2+y^2≥0.95,vars,pars]},Az,{x,y}∈mesh]

AzDIF=DiscontinuousInterpolatingFunction[AzSol]

B-Field

Bfield=Curl[{0,0,AzDIF[x,y]},{x,y,z}];

BfieldNormRange=MinMax[Sqrt[Total[Bfield^2]]/. getValuesOnGrid] ;
Show[DensityPlot[Sqrt[Total[Bfield^2]],{x,y}∈mesh,],VectorPlot[Bfield[[1;;2]],{x,y}∈mesh,],bmesh["Wireframe"]]

New Application Examples

Inductive Heating

Martin’s Talk: Transient Cerebral Aneurysm

Future of PDE Modeling

Plastic deformation

Axisymmetric Fluid Flow

Ricardo’s talk: Investigation of NN based PDE solving

NDSolve rocks!

Thanks

Init

Needs["NDSolve`FEM`"]
Needs["OpenCascadeLink`"]
$HistoryLength=0;

mr=;
bmesh=ToBoundaryMesh[mr];
airCoordinate1={0.,0.};
rotorCoordinate={0.,0.3125};
airCoordinate2={0,0.8};
statorCoordinate1={0,0.95};
fingerCoordinates={{0.25,0.6},{0.6,0.25},{0.6,-0.25},{0.25,-0.6},{-0.25,-0.6},{-0.6,-0.25},{-0.6,0.25},{-0.25,0.6}};
tipCoordinates={{0.18,0.44},{0.44,0.18},{0.44,-0.18},{0.18,-0.44},{-0.18,-0.44},{-0.44,-0.18},{-0.44,0.18},{-0.18,0.44}};
coilCoordinates={{0.15,0.64},{0.35,0.56},{0.56,0.35},{0.64,0.15},{0.64,-0.15},{0.56,-0.35},{0.35,-0.56},{0.15,-0.64},{-0.15,-0.64},{-0.35,-0.56},{-0.56,-0.35},{-0.64,-0.15},{-0.64,0.15},{-0.56,0.35},{-0.35,0.56},{-0.15,0.64}};
markerColors={Blue,Red,Gray,Orange};
markerCoordinates={{airCoordinate1,airCoordinate2},{rotorCoordinate},{statorCoordinate1,fingerCoordinates,tipCoordinates},{coilCoordinates}};
Show[
bmesh["Wireframe"],
Graphics[MapThread[{PointSize[0.02],#1,Point/@#2}&,{markerColors,markerCoordinates}]]];
coilMarkers={4,4,4,4,5,5,5,5,6,6,6,6,5,5,5,5};
Short[markerSpecification=Join[{{airCoordinate1,1},{airCoordinate2,1},{rotorCoordinate,2},{statorCoordinate1,3}},{#,3}&/@fingerCoordinates,{#,3}&/@tipCoordinates,MapThread[{#1,#2}&,{coilCoordinates,coilMarkers}]]];
mesh=ToElementMesh[bmesh,"RegionMarker"->markerSpecification];

getValuesOnGrid=(head:DiscontinuousInterpolatingFunction |InterpolatingFunction)[a__][d__]:>head[a]["ValuesOnGrid"];