Modeling Multi-Physics with PDEs

Yi-Lin Chiu,

Software Developer, WRI.

Join the Conversation #WolframTechConf

Partial differential equations (PDEs) Model

Resources: paclet:PDEModels/tutorial/PDEModelsOverview

Workflow

Sequential PDE Modeling Example

1. Electrostatics Model:

∇·(-G·∇V)=0

2. Heat Transfer Model:

ρ (Subscript[C, p])((∂T)/(∂t))+∇·(-k ∇T)=G((\[LeftBracketingBar]∇V\[RightBracketingBar])^( 2))

(({{V, :, Electric potential}, {T, :, Temperature}})

Model Setup: Electrostatics Model

Model Setup: Electrostatics Model

Domain

Ω=ToElementMesh[];

Electrostatics Model

∇·(-G·∇V)=0

ElectricModel=(InactiveDiv[(G (InactiveGrad[V[x,y,z], {x,y,z}, Inactive])), {x,y,z}, Inactive]);

Boundary Conditions

(Subscript[Γ, volt])=
{DirichletCondition[V[x,y,z]==(Subscript[V, o]),leftEndQ[x,y,z]],
DirichletCondition[V[x,y,z]==0,rightEndQ[x,y,z]]};

colors={,,};
Legended[Ω["Wireframe"["MeshElement"->"BoundaryElements","MeshElementStyle"->(Directive[FaceForm[#1],EdgeForm[]]&)/@colors,ImageSize->{750,450}]],SwatchLegend[{,Orange,},{Text[Style[" (Subscript[V, o]) [V]",24,]],Text[Style[" 0 [V]",24,Orange]]}]]

Solve: Electrostatics Model

Solve: Electrostatics Model

∇·(-G·∇V)=0

Solve

Vfun=NDSolveValue[{ElectricModel==0,(Subscript[Γ, volt])},V,{x,y,z}∈Ω]

Post-Processing: Electrostatics Model

Post-Processing: Electrostatics Model

VRange=MinMax[Vfun["ValuesOnGrid"]];
legendBar=BarLegend[{"SunsetColors",VRange},LegendLabel->Text[Style["[V]",24,Opacity[0.6]]],LabelStyle->{FontSize->20},LegendMarkerSize->400];
options={AspectRatio->Automatic,PerformanceGoal->"Quality",PlotPoints->50,Mesh->None,PlotTheme->"Detailed",PlotLegends->None,AxesLabel->{x,y,z},LabelStyle->Directive[20],ColorFunctionScaling->False,ImageSize->750,PlotLabel->Text[Style["Electric Potential Field: V(x,y,z)",32]],Ticks->{Automatic,Automatic,{-r,r}}};
Legended[RegionPlot3D[Ω,ColorFunction->Function[{x,y,z},
ColorData[{"SunsetColors",VRange}][Vfun[x,y,z]]],Evaluate[options]],legendBar]

Post-Processing: Electrostatics Model

Electric Heat Generation

Q=G((\[LeftBracketingBar]∇V\[RightBracketingBar])^( 2))

Q=G*((Norm[Grad[Vfun[x,y,z],{x,y,z}]])^2);

legendBar=BarLegend[{"TemperatureMap",{0,(10^7)}},LegendLabel->Text[Style["[W/\!\(\*SuperscriptBox[\(m\), \(3\)]\)]",Opacity[0.6],24]],LabelStyle->{FontSize->15},LegendMarkerSize->200];
options={LabelStyle->Directive[20],ImageSize->750,PlotRange->{0,1.2*(10^7)},ColorFunction->ColorData[{"TemperatureMap",{0,(10^7)}}],ContourStyle->None,ColorFunctionScaling->False,Contours->50,PlotLabel->Text[Style["Q(x,y)",32]],PlotPoints->100,FrameLabel->{"x","y"},AspectRatio->Automatic};
crossSection=RegionUnion[Rectangle[{0,-r},{L,r}],Rectangle[{R+r,2 R-r},{L,2 R+r}],Rectangle[{R+r,4 R-r},{L+R+r,4 R+r}],Rectangle[{L+2 R,-r},{L+2 R+2 r,3 R}],Annulus[{L,R},{R-r,R+r},{((3 π)/2),((5 π)/2)}],Annulus[{R+r,3 R},{R-r,R+r},{(π/2),((3 π)/2)}],Annulus[{L+R+r,3 R},{R-r,R+r},{0,(π/2)}]];
bmeshWireframe=ToBoundaryMesh[crossSection]["Wireframe"];
Show[Legended[ContourPlot[Q/.{z->0},{x,y}∈crossSection,Evaluate[options]],legendBar],bmeshWireframe]

Total Heat Generation

(Subscript[q, gain])=(∫_( Ω)) Q(x,y,z)\[DifferentialD]Ω

(Subscript[q, gain])=NIntegrate[Q,{x,y,z}∈Ω]

Model Setup: Heat Transfer Model

Model Setup: Heat Transfer Model

Heat Transfer Model

Q=G((\[LeftBracketingBar]∇V\[RightBracketingBar])^( 2))

Q=G*((Norm[Grad[Vfun[x,y,z],{x,y,z}]])^2);

ρ (Subscript[C, p])((∂T)/(∂t))+∇·(-k ∇T)=Q

HeatModel3D=
TimeHeatTransferModel[T[t,x,y,z],t,{x,y,z},k,ρ,Cp,"NoFlow",Q];

Initial Conditions

ic={T[0,x,y,z]==(Subscript[T, amb])};

Boundary Conditions

(Subscript[Γ, temp])=DirichletCondition[T[t,x,y,z]==(Subscript[T, amb]),leftEndQ[x,y,z]||rightEndQ[x,y,z]];

(Subscript[Γ, radiation])=NeumannValue[ε σ ((Subscript[T, amb]^4)-((T[t,x,y,z])^4)),];

(Subscript[Γ, convection])=NeumannValue[h ((Subscript[T, amb])-T[t,x,y,z]),];

colors={Orange,Blue,Blue};
Legended[Ω["Wireframe"[ImageSize->{750,450},"MeshElement"->"BoundaryElements","MeshElementStyle"->(Directive[FaceForm[#1],EdgeForm[]]&)/@colors]],SwatchLegend[{Blue,Orange},{Text[Style[" 300 [K]",24,Blue]],Text[Style["  Radiation + Convection",24,Orange]]}]]

Solve: Heat Transfer Model

ρ (Subscript[C, p])((∂T)/(∂t))+∇·(-k ∇T)=Q

Solve

Tfun=
NDSolveValue[{HeatModel3D==(Subscript[Γ, radiation])+(Subscript[Γ, convection]),(Subscript[Γ, temp]),ic},T,{t,0,300},{x,y,z}∈Ω];

Post-Processing: Heat Transfer Model

TRange=MinMax[Tfun["ValuesOnGrid"]];
legendBar=BarLegend[{"TemperatureMap",TRange},LegendLabel->Text[Style["[K]",Opacity[0.6],24]],LabelStyle->{FontSize->20},LegendMarkerSize->400];
options={AspectRatio->Automatic,PerformanceGoal->"Quality",PlotPoints->50,Mesh->None,PlotTheme->"Detailed",PlotLegends->None,AxesLabel->{x,y,z},LabelStyle->Directive[20],ColorFunctionScaling->False,ImageSize->{800,600},PlotLabel->Text[Style["Temperature Field: T(x,y,z)",32]],Ticks->{Automatic,Automatic,{-r,r}}};
nframes=50;
frames=Table[Labeled[Legended[RegionPlot3D[Ω,ColorFunction->Function[{x,y,z},
ColorData[{"TemperatureMap",TRange}][Tfun[t,x,y,z]]],Evaluate[options]],legendBar],Text[Style[Row[{"Time"," = ",Round[t]," s"}],Opacity[0.6],24]],Spacings->{0, 0}],{t,0,300,(300/nframes)}];
ListAnimate[frames,SaveDefinitions->True]

Post-Processing: Heat Transfer Model

Heat Loss

({{(Subscript[q, radiation]), =, (∫_( Γ)) ε σ ((T^( 4))-((Subscript[T, amb])^4))\[DifferentialD]Γ}, {(Subscript[q, convection]), =, (∫_( Γ)) h (T-(Subscript[T, amb]))\[DifferentialD]Γ}, {(Subscript[q, conduction]), =, (∫_( Γ)) k∇T\[DifferentialD]Γ}})

(Subscript[q, radiation])=
NIntegrate[Piecewise[{{ε σ (((Tfun[300,x,y,z])^4)-(Subscript[T, amb]^4)),(Subscript[Γ, side])},{0,True}}],{x,y,z}∈bmesh]

(Subscript[q, convection])=

(Subscript[q, conduction])=+

Sequential PDE Modeling

({{Electrostatics Model, ⟶, Heat Transfer Model}})

({{(Electric potential: V) ,                                              , (Temperature:T)}})

({{Model A, (Overscript[⟷, ??]), Model B}})

Workflow

Fully-coupled PDE Modeling Example

Heat Transfer Model:

ρ (Subscript[C, p])V·∇T+∇·(-k ∇T)=0

Fluid Dynamics Model:

({{ρ V·∇V+∇·(-μ ∇V)+∇p, =, ρ·g·α·(T-(Subscript[T, in]))}, {∇·V, =, 0}})

(({{V, :, Flow Velocity}, {T, :, Temperature}})

Graphics[{Arrowheads[0.05`],{Opacity[0.5`],LightBlue,Rectangle[{0.,0},{0.035`,0.045`}]},{Thick,Circle[{0.025`,0.015`},0.0025`],Circle[{0.01`,0.015`},0.0025`,{(π/2),((3 π)/2)}]},{Thick,Circle[{0.025`,0.03`},0.0025`],Circle[{0.01`,0.03`},0.0025`,{(π/2),((3 π)/2)}]},{Translate[{LightOrange,Disk[{0,0.015`},0.0025`]},Table[{i,0},{i,0.01,0.035`,0.015`}]]},{Translate[{LightOrange,Disk[{0,0.03`},0.0025`]},Table[{i,0},{i,0.01,0.035`,0.015`}]]},{Thick,Dashed,Line[{{0.01`,0},{0.0175`,0},{0.0175`,0.045`},{0.01`,0.045`},{0.01`,0.0325`}}],Line[{{0.01`,0.0275`},{0.01`,0.0175`}}],Line[{{0.01`,0.0125`},{0.01`,0}}],Circle[{0.01`,0.015`},0.0025`,{-(π/2),(π/2)}],Circle[{0.01`,0.03`},0.0025`,{-(π/2),(π/2)}]},{Red,Arrowheads[0.03`],Thick,Arrow[{{0.0125`,0.024`},{0.011`,0.0275`}}],Arrow[{{0.0125`,0.021`},{0.011`,0.0175`}}],Text[Style["tubes",24],{0.014`,0.0225`}]},{Thick,Red,Translate[Arrow[{{0.,0.037`},{0.,0.044`}}],Table[{i,0},{i,0.005,0.035,0.00875`}]],Text[Style["outflow",24],{0.027,0.035`}]},{Thick,Blue,Translate[Arrow[{{0.,0.001`},{0.,0.008`}}],Table[{i,0},{i,0.005,0.035`,0.00875`}]],Text[Style["inflow",24],{0.027,0.01`}]}},ImageSize->350,PlotLabel->Text[Style["Heat Exchanger",32]],Frame->True,FrameTicks->None,LabelStyle->Directive[20]]

Model Setup

Domain

Ω=;

Heat transfer model + Fluid dynamics model

HeatModel2D=
HeatTransferModel[T[x,y],{x,y},(Subscript[k, water]),(Subscript[ρ, water]),(Subscript[Cp, water]),{u[x,y],v[x,y]},"NoSource"];

FlowModel2D=
FluidDynamicsModel[{u[x,y],v[x,y]},p[x,y],{x,y},(Subscript[μ, water]),(Subscript[ρ, water]),{0,α*g*(Subscript[ρ, water])*(T[x,y]-(Subscript[T, in]))}];

Boundary Conditions

(Subscript[Γ, inlet])=DirichletCondition[{u[x,y]==0,v[x,y]==(Subscript[v, in])},y==0];

(Subscript[Γ, outlet])=;

(Subscript[Γ, wall])=;

(Subscript[Γ, sym])=;

(Subscript[Γ, temp])={,};

BCs={(Subscript[Γ, temp]),(Subscript[Γ, inlet]),(Subscript[Γ, outlet]),(Subscript[Γ, wall]),(Subscript[Γ, sym])};

Solve

Coupled PDE models

CoupledPDE=
{FlowModel2D=={0,0,0},HeatModel2D==0,BCs};

Solve

{ufun,vfun,pfun,Tfun}=
NDSolveValue[CoupledPDE,{u,v,p,T},{x,y}∈Ω,];

Post-Processing

TRange=MinMax[Tfun["ValuesOnGrid"]];
legendBar=;
options={};
boundaryHighlight=;
g1=Show[boundaryHighlight,VectorPlot[{ufun[x,y],vfun[x,y]},{x,y}∈ufun["ElementMesh"],AspectRatio->Automatic]];
g2=Show[ContourPlot[Tfun[x,y],{x,y}∈ufun["ElementMesh"],Evaluate[options]],boundaryHighlight];
g3=Show[Legended[ContourPlot[Tfun[x,y],],legendBar],VectorPlot[{ufun[x,y],vfun[x,y]},{x,y}∈ufun["ElementMesh"],PlotTheme->"Scientific"],boundaryHighlight];
GraphicsRow[{g1,g2,g3},ImageSize->{1000,500},Spacings->0]

Post-Processing

Mean Temperature at Outlet

(Subscript[T, mean])=(1/(\[LeftBracketingBar]Γ\[RightBracketingBar]))(∫_( Γ)) T(x,y)\[DifferentialD]Γ

(Subscript[T, out])=NIntegrate[Tfun[x,H],{x,0,W}]/W

Resources: paclet:PDEModels/tutorial/PDEModelsOverview

Thank you for your attention!

paclet:PDEModels/tutorial/PDEModelsOverview

yilinc@wolfram.com

Initialization Code

Example: Tungsten Wire (electrostatics)

Load the finite element package.

Needs["NDSolve`FEM`"]

Define the electrostatics model function.

ElectrostaticsModel[V_,X_List,G_]:=Module[{a,factor},
factor=-G*IdentityMatrix[Length[X]];
a=PiecewiseExpand[Piecewise[{{factor,True}}]];
Inactive[Div][G Inactive[Grad][V,X],X]]

Specify geometry parameters.

L=0.1;
r=0.0025;
R=0.0075;

Specify model parameters.

G=1.79*(10^7);
(Subscript[V, o])=0.2;

Specify the path of the wire.

arc1=Table[R*{Cos[θ],Sin[θ],0.},{θ,-π/2,π/2,π/30}];
arc2=Reverse[-arc1];
arc3=Table[R*{Cos[θ],Sin[θ],0.},{θ,π/2,0,-π/30}];
path=;

Set up the functions leftEndQ and rightEndQ to represent both ends of the wire.

tolerance=0.01r;
leftEndQ[x_,y_,z_]:=Abs[x]≤tolerance&&(y^2)+(z^2)≤(r^2)
rightEndQ[x_,y_,z_]:=Abs[y+r]≤tolerance&&(((x-(L+2R+r)))^2)+(z^2)≤(r^2)
(Subscript[Γ, side])=Not[leftEndQ[x,y,z]||rightEndQ[x,y,z]];

Example: Tungsten Wire (heat transfer)

Define the heat transfer model function.

ClearAll[HeatTransferModel]
HeatTransferModel[T_,X_List,k_,ρ_,Cp_,Velocity_,Source_]:=Module[{V,Q,a=k},
V=If[Velocity==="NoFlow",0,ρ*Cp*Velocity.Inactive[Grad][T,X]];
Q=If[Source==="NoSource",0,Source];
If[ FreeQ[a, _?VectorQ],a=a*IdentityMatrix[Length[X]]];
If[ VectorQ[a],a=DiagonalMatrix[a]];
(* Note the - sign in the operator *)
a=PiecewiseExpand[Piecewise[{{-a,True}}]];
Inactive[Div][a.Inactive[Grad][T,X],X]+V-Q]

TimeHeatTransferModel[T_,TimeVar_,X_List,k_,ρ_,Cp_,Velocity_,Source_]:=ρ*Cp*D[T,{TimeVar,1}]+HeatTransferModel[T,X,k,ρ,Cp,Velocity,Source]

Specify model parameters.

ρ=1.93*(10^4);
Cp=134;
k=175;
h=20;
ε=1;
(Subscript[T, amb])=300;

Extract the value of Stefan-Boltzmann Constant.

σ=First[UnitConvert[Quantity["StefanBoltzmannConstant"]]];

Construct the boundary element mesh for the boundary integration.

bmesh=ToBoundaryMesh[Ω];

Example: Heat Exchanger

Define the fluid dynamics model function.

FluidDynamicsModel[V_List,p_,X_List,μ_,ρ_,VolumeForce_]:=Module[{F},
F=If[VolumeForce==="NoForce",Table[0,Length[X]],VolumeForce];
Append[
Table[ρ V.Inactive[Grad][V[[i]],X]+Inactive[Div][-μ Inactive[Grad][V[[i]],X],X]+D[p,X[[i]]],{i,Length[X]}]-F,
Div[V,X]]]

Specify geometry parameters.

W=0.005;
H=0.045;
(Subscript[r, cyl])=0.0025;

Specify model parameters.

(Subscript[ρ, water])=QuantityMagnitude[ThermodynamicData["Water","Density"]];
(Subscript[k, water])=QuantityMagnitude[ThermodynamicData["Water","ThermalConductivity"]];
(Subscript[μ, water])=QuantityMagnitude[ThermodynamicData["Water","Viscosity"]];
(Subscript[Cp, water])=QuantityMagnitude[ThermodynamicData["Water","IsobaricHeatCapacity"]];
(Subscript[v, in])=5*10^-3;
α=1.8*10^-4;
(Subscript[T, in])=300;
(Subscript[T, wire])=375;
g=9.81;

Set up the functions wire1 and wire2 to represent the boundary of two wires.

wire1=(x^2)+(((y-0.015))^2)≤((Subscript[r, cyl])^2);
wire2=(x^2)+(((y-0.03))^2)≤((Subscript[r, cyl])^2);

Code for post-processing visualization.

(*TRange=MinMax[Tfun["ValuesOnGrid"]];
legendBar=;
options={};
boundaryHighlight=;
g1=Show[boundaryHighlight,VectorPlot[{ufun[x,y],vfun[x,y]},{x,y}∈ufun["ElementMesh"],AspectRatio->Automatic]];
g2=Show[ContourPlot[Tfun[x,y],{x,y}∈ufun["ElementMesh"],Evaluate[options]],boundaryHighlight];
g3=Show[Legended[ContourPlot[Tfun[x,y],],legendBar],VectorPlot[{ufun[x,y],vfun[x,y]},{x,y}∈ufun["ElementMesh"],PlotTheme->"Scientific"],boundaryHighlight];
GraphicsRow[{g1,g2,g3},ImageSize->{1000,500},Spacings->0]*)