In[]:=
Needs["NDSolve`FEM`"];
In[]:=
ResourceFunction["FEMAddOnsInstall"][]
Out[]=
PacletObject
Name: FEMAddOns
Version: 1.4.5

In[]:=
Needs["FEMAddOns`"]
In[]:=
ClearAll[length,height,offsetx1,offsetx2,offsety];​​sizes={length->21/10,height->15/10,​​offsetx1->1/10,offsetx2->1/10,​​offsety1->1/10,offsety2->2/10​​};
In[]:=
​​Geometry
Out[]=
Geometry
In[]:=
Ω1=Rectangle[{0,0},{length,height}]/.sizes;​​Ω2=Rectangle[{0,0}+{0.2,1.5-0.049},{0.69,0.049}+{0.2,1.5-0.049}];​​​​ΩD1=RegionDifference[Ω1,Ω2]
Out[]=
Polygon
Number of points: 8
Embedding dimension: 2

In[]:=
Ω3=Rectangle[{0,0}+{1.2,1.5-0.049},{0.69,0.049}+{1.2,1.5-0.049}];​​ΩD2=RegionDifference[ΩD1,Ω3]
Out[]=
Polygon
Number of points: 12
Embedding dimension: 2

In[]:=
RegionPlot[ΩD2,AspectRatio->Automatic]
In[]:=
Mesh​​
Out[]=
Mesh
In[]:=
Ω=DiscretizeRegion[ΩD2,MaxCellMeasure0.01,AccuracyGoal4]
In[]:=
meshes=ToBoundaryMesh/@{Ω1,Ω2,Ω3};
In[]:=
bmesh=FEMUtils`BoundaryElementMeshJoin@@meshes;
In[]:=
boxCoordinate={1.05,0.8};​​PCMCoordinate1={0.5,1.46};​​PCMCoordinate2={1.5,1.46};​​markerColors={Blue,Orange,Green};​​markerCoordinate={{boxCoordinate},{PCMCoordinate1},{PCMCoordinate2}};
In[]:=
markerspec=MapIndexed[{First@#1,First@#2}&,markerCoordinate]
Out[]=
{{{1.05,0.8},1},{{0.5,1.46},2},{{1.5,1.46},3}}
In[]:=
​​mesh=ToElementMesh[bmesh,"RegionMarker"->markerspec,"MaxCellMeasure"->5*10^-3];
In[]:=
mesh["Wireframe"[]]
In[]:=
mesh["Wireframe"["MeshElementStyle"->Map[Directive[FaceForm[#],EdgeForm[]]&,markerColors]]]
In[]:=
​
In[]:=
Show[{bmesh["Wireframe"],Graphics[MapThread[{PointSize[0.02],#1,Point/@#2}&,{markerColors,markerCoordinate}]]}]
In[]:=
MaterialProperties
Out[]=
MaterialProperties
In[]:=
region=<|"air"->1,"pcm"->2|>
Out[]=
air1,pcm2
In[]:=
vars={T[t,x,y],t,{x,y}};​​HeatTransferPDEComponent[vars,<|"MassDensity"->ρ,"SpecificHeatCapacity"->Cp,"ThermalConductivity"->k|>]
Out[]=
∇
{x,y}
·(-k).
∇
{x,y}
T[t,x,y]+Cpρ
(1,0,0)
T
[t,x,y]
In[]:=
airpara=
ρ
air
->QuantityMagnitude[ThermodynamicData["Air","Density"]],​​
Cp
air
->QuantityMagnitude[ThermodynamicData["Air","IsobaricHeatCapacity"]],​​
k
air
->QuantityMagnitude[ThermodynamicData["Air","ThermalConductivity"]]
Out[]=
{
ρ
air
1.20458,
Cp
air
1006.14,
k
air
0.0258738}
In[]:=
HS30Npara=
Cp
pcmS
->2.1,
Cp
pcmL
->2.7,​​
k
pcmS
->2.34,
k
pcmL
->0.6,​​
ρ
pcmS
->1430,latent->200,​​ tpmin->-34,tpmax->-26;
In[]:=
PUFpara={ρ->50,cp->1500,k->0.025};
In[]:=
​
In[]:=
pars=<||>;​​pars["MassDensity"]=If[ElementMarker==region["air"],
ρ
air
,
ρ
pcmS
];​​pars["SpecificHeatCapacity"]=If[ElementMarker==region["air"],
k
air
,
k
pcmS
];​​pars["ThermalConductivity"]=If[ElementMarker==region["air"],
k
air
,
k
pcmS
*IdentityMatrix[2]];
In[]:=
op=HeatTransferPDEComponent[vars,pars]
Out[]=
∇
{x,y}
·-IfElementMarker1,
k
air
,
k
pcmS
IdentityMatrix[2].
∇
{x,y}
T[t,x,y]+IfElementMarker1,
k
air
,
k
pcmS
IfElementMarker1,
ρ
air
,
ρ
pcmS

(1,0,0)
T
[t,x,y]
In[]:=
​​InitialandBoundaryCondition
Out[]=
andBoundaryConditionInitial
In[]:=
Tpcm=-40;​​Tair=-40;
In[]:=
ic={u[0,x,y]==0,​​v[0,x,y]==0,​​p[0,x,y]==0,​​T[0,x,y]==Piecewise[{{Tpcm,ElementMarker==region["pcm"]},{Tair,ElementMarker==region["air"]}}]};
In[]:=
​​
Γ
convectiveLeft
=HeatTransferValue[x0,{T[t,x,y],t,{x,y}},<||>,<|"HeatTransferCoefficient"0.7,"AmbientTemperature"40|>];​​
Γ
convectiveRight
=HeatTransferValue[xlength/.sizes,{T[t,x,y],t,{x,y}},<||>,<|"HeatTransferCoefficient"0.7,"AmbientTemperature"40|>];​​
Γ
convectiveTop
=HeatTransferValue[yheight/.sizes,{T[t,x,y],t,{x,y}},<||>,<|"HeatTransferCoefficient"0.7,"AmbientTemperature"40|>];​​
Γ
convectiveBottom
==HeatTransferValue[y0/.sizes,{T[t,x,y],t,{x,y}},<||>,<|"HeatTransferCoefficient"0.7,"AmbientTemperature"40|>];
In[]:=
bcstemperatures={​​DirichletCondition[T[t,x,y]==-40,{x,y}∈Ω2||{x,y}∈Ω3]​​}​​
Out[]=
{DirichletCondition[T[t,x,y]-40,{x,y}∈Rectangle[{0.2,1.451},{0.89,1.5}]||{x,y}∈Rectangle[{1.2,1.451},{1.89,1.5}]]}
In[]:=
ClearAll[ν,ϵ,Pr,Ra];​​g=9.8;​​α=0.0034;​​ν=1.48*10^-5;​​β=2.17*10^-5;​​L=(4*2.1*1.5)/(2*(2.1+1.5));​​Temp=40;​​Subscript[T,∞]=0;​​Ra->((g*α)/(ν*β))*(Temp-Subscript[T,∞])*L^3;​​parameters={μ->Sqrt[Pr/Ra],ϵ->1/Sqrt[Pr*Ra]}/.{Pr->0.7,Ra->((g*α)/(ν*β))*(Temp-Subscript[T,∞])*L^3};​​
In[]:=
op=​​​​
(1,0,0)
u
[t,x,y]​​+Inactive[Div][(-νInactive[Grad][u[t,x,y],{x,y}]),{x,y}]​​+{u[t,x,y],v[t,x,y]}.Inactive[Grad][u[t,x,y],{x,y}]​​+
(0,1,0)
p
[t,x,y],​​​​
(1,0,0)
v
[t,x,y]​​+Inactive[Div][(-νInactive[Grad][v[t,x,y],{x,y}]),{x,y}]​​+{u[t,x,y],v[t,x,y]}.Inactive[Grad][v[t,x,y],{x,y}]​​+
(0,0,1)
p
[t,x,y]-T[t,x,y],​​​​
(0,1,0)
u
[t,x,y]+
(0,0,1)
v
[t,x,y],​​​​
(1,0,0)
T
[t,x,y]​​+Inactive[Div][(-ϵInactive[Grad][T[t,x,y],{x,y}]),{x,y}]​​+{u[t,x,y],v[t,x,y]}.Inactive[Grad][T[t,x,y],{x,y}]​​​​/.parameters;