New in:
12.1
| Modified in:
| Obsolete in:
| Excised in:
Created by: ruebenko on 02-04-2020 16:30:48
Categorization

Synonyms

Keywords

Details


Steady-state solution of the viscous heat equations

Michele Simoncelli, Nicola Marzari, and Andrea Cepellotti
In this notebook we use the finite-element method to solve the viscous heat equations, two differential equations for heat conduction in crystals that generalize Fourier's law and explain why heat propagation can become fluid like, rather than diffusive, in electronic or phononic devices. A theoretical derivation of these equations can be found in Ref. [1].
Load the FEM package.
In[202]:=
Needs["NDSolve`FEM`"]

Domain and mesh setup

Definition of the 2D domain, which corresponds to the in-plane (x,y) section of a 3D domain infinite in the z direction. Dimensions are in micrometer, see [1] for details.
Set up geometric parameters.
In[203]:=
yjumpum=3.0;​​xjumpum=3.0;​​lengthBoxum=15.0;​​widthBoxum=9.0;
Create a region and a mesh from the region.
In[207]:=
Ω=RegionUnion[Rectangle[{0,yjumpum},{xjumpum,widthBoxum-yjumpum}],Rectangle[{xjumpum,0},{lengthBoxum,widthBoxum}]];​​emesh=ToElementMesh[Ω,"MaxCellMeasure"{"Length"lengthBoxum/40.0},MeshQualityGoal1];​​emesh["Wireframe"]
Out[209]=

Parameters entering into the viscous heat equations

We are assuming that the sample is made of polycrystalline graphite with a grain size equal to 10 μm. Under these conditions, Ref. [1] gives the following numerical values for the parameters entering into the viscous heat equations:
Set up material parameters.
In[210]:=
μ
a
=0.000907647;(*inpicog/(micrometer*ns),thisaccountsforagrainsizeequalto10μm*)​​
μ
b
=0.000485155;(*inpicog/(micrometer*ns),thisaccountsforagrainsizeequalto10μm*)​​
μ
c
=0.000401477;(*inpicog/(micrometer*ns),thisaccountsforagrainsizeequalto10μm*)​​W=3.252493000;(*inmicrometer/ns*)​​
C
V
=0.000184524;(*Specificheat,inpicog/(micrometer*ns^2*K)*)​​A=0.000320256;(*Specificmomentum,inpg/micrometer^3*)​​
D
u
=0.979531000;(*Umklappmomentumdissipationtimenscale,inns^-1,thisaccountsforagrainsizeequalto10μm*)​​κ=0.002721340;(*Thermalconductivity,inpicog*micrometer/(K*ns^3),thisaccountsforagrainsizeequalto10μm*)​​
T
eq
=70.000000000;(*Equilibriumtemperaturearoundwhichatemperatureperturbationisapplied,inK*)​​δT=10.0(*temperatureperturbation,inK*);

Equations setup

The viscous heat equations are given by Eq. (10) and Eq. (11) of Ref. [1]. Here we consider the steady-state and a sample infinite in the z direction, therefore the equations that we have to solve assume the following form:
Set up the viscous heat equation.
In[220]:=
EqUx=
(1,0)
T
[x,y]*W*
C
V
*A
T
eq
-
μ
a
*
(2,0)
Ux
[x,y]-
μ
b
*
(0,2)
Ux
[x,y]-
μ
c
*
(1,1)
Uy
[x,y]+(
D
u
*A)*Ux[x,y];​​EqUy=
(0,1)
T
[x,y]*W*
C
V
*A
T
eq
-
μ
b
*
(2,0)
Uy
[x,y]-
μ
a
*
(0,2)
Uy
[x,y]-
μ
c
*
(1,1)
Ux
[x,y]+(
D
u
*A)*Uy[x,y];​​EqT=W*
T
eq
*A*
C
V
*
(1,0)
Ux
[x,y]+
(0,1)
Uy
[x,y]-κ*
(2,0)
T
[x,y]-κ*
(0,2)
T
[x,y];​​Operator={EqUx,EqUy,EqT};
We use no-slip condition for the drift velocity (Ux=0 and Uy=0 at all the boundaries). For the temperature T(x=0μm)=80 K, T(x=15μm)=60 K are set.
Set up Dirichlet boundary conditions.
In[224]:=
bcs={DirichletCondition[{Ux[x,y]0.,Uy[x,y]0.},True],​​DirichletCondition[T[x,y]
T
eq
-δT,(xlengthBoxum)&&(0.0≤y≤widthBoxum)],​​DirichletCondition[T[x,y]
T
eq
+δT,(x0.0)&&(yjumpum≤y≤(widthBoxum-yjumpum))]​​};
All the boundaries with x≠0 or x≠15 μm are adiabatic which is modeled with a Neumann zero boundary condition. Since the Neumann zero boundary condition is the default condition nothing needs to be specified.
Set up the PDE.
In[225]:=
pde=Operator{0,0,0};

Equations solving

To solve the equations
NDSolve
is used.
NDSolve
needs some help to figure out what the dependent variable ordering is. This is specified with the option
DependentVariables
. More information on the dependent variable ordering can be found in the documentation section
OrderingofDependentVariableNames
.
Solve the equations.
In[226]:=
{Ux0,Uy0,T0}=NDSolveValue[{pde,bcs},{Ux,Uy,T},{x,y}∈emesh,DependentVariables{Ux,Uy,T}];

Analysis of the viscous heat equations’ solution

Plot of the temperature profile obtained from the solution of the viscous heat equations.
Plot the temperature profile.
In[227]:=
DensityPlot[T0[x,y],{x,0,lengthBoxum},{y,0,widthBoxum},​​ColorFunction"LightTemperatureMap",​​PlotLegendsBarLegend[Automatic,LegendLabelPlaced["Temperature (K)",Below]],​​AspectRatiowidthBoxum/lengthBoxum,​​ImageSize{500,300},​​AxesTrue,​​AxesOrigin{0,0},​​PlotRange{{0,lengthBoxum},{0,widthBoxum}},​​FrameLabel{"x (μm)","y (μm)"},​​LabelStyleDirective[12,Black,Plain],​​PlotLabel"Solution of the viscous heat equations: temperature field"]
Out[227]=
Temperature (K)
Plot the drift-velocity profile resulting from the solution of the viscous heat equations.
Plot the velocity profile.
In[228]:=
StreamDensityPlot[{{Ux0[x,y],Uy0[x,y]},
0.5
(
2
(Ux0[x,y])
+
2
(Uy0[x,y])
)
},{x,y}∈Ω,​​StreamPointsFine,​​StreamStyleWhite,​​PlotLegendsBarLegend[Automatic,LegendLabelPlaced["drift velocity module (μm/ns)",Below]],​​ColorFunctionColorData["BlueGreenYellow"],​​AspectRatiowidthBoxumlengthBoxum,ImageSize{500,300},​​FrameTrue,​​FrameLabel{"x (μm)","y (μm)"},​​AxesTrue,​​AxesOrigin{0,0},​​PlotRange{{0,lengthBoxum},{0,widthBoxum}},​​LabelStyleDirective[12,Black,Plain],​​FrameStyleBlack,​​PlotLabel"Solution of the viscous heat equations: drift-velocity field"]

References

More About
XXXX
Related Tutorials
XXXX
Related Wolfram Training Courses
XXXX