Model Energy Transport Dynamics

Solve Navier-Stokes equations coupled with a heat transfer equation over a rectangular region with cutout

Consider a multi-physics problem that couples the Navier–Stokes equation with a heat equation. Here

is the vector-valued velocity field,

is the pressure and

the identity matrix.

and

are the density and dynamic viscosity, respectively.

is the heat capacity of the liquid and

the temperature:

∂u ∂t +∇(ν∇u)+ρu∇u+Ip-{0,T}	=	0
∇u	=	0
∂T ∂t +∇(-ϵ∇T)+u∇T	=	0

Set up a viscous Navier–Stokes equation coupled to a heat equation making use of a Boussinesq approximation. Use the material parameters specified:

In[1]:=

op=

(1,0,0)

[t,x,y]+

∇

{x,y}

·-ν

∇

{x,y}

u[t,x,y]+{u[t,x,y],v[t,x,y]}.

∇

{x,y}

u[t,x,y]+

(0,1,0)

[t,x,y],

(1,0,0)

[t,x,y]+

∇

{x,y}

·-ν

∇

{x,y}

v[t,x,y]+{u[t,x,y],v[t,x,y]}.

∇

{x,y}

v[t,x,y]+

(0,0,1)

[t,x,y]-T[t,x,y],

(0,1,0)

[t,x,y]+

(0,0,1)

[t,x,y],

(1,0,0)

[t,x,y]+

∇

{x,y}

·-ϵ

∇

{x,y}

T[t,x,y]+{u[t,x,y],v[t,x,y]}.

∇

{x,y}

T[t,x,y];

Specify a Prandtl number (

Pa)

and a Rayleigh number (

) and computer the dynamic viscosity (ν) and the heat capacity of the liquid (ϵ):

In[2]:=

parameters={νSqrt[Pr/Ra],ϵ1/Sqrt[Pr*Ra]}/.{Pr7.1,Ra2*10^5}

Out[2]=

{ν0.00595819,ϵ0.000839181}

Update the set of equations with the material parameters specified:

In[3]:=

op=op/.parameters;

Specify a rectangular region with cutouts:

In[4]:=

sizes={length2,height1,thickness1/6,offset1/4};Ω=RegionDifference[RegionDifference[Rectangle[{0,0},{length,height}],Rectangle[{length/2-offset-thickness,height/3},{length/2-offset,height}]],Rectangle[{length/2+offset,0},{length/2+offset+thickness,2height/3}]]/.sizes;

Visualize the region:

In[5]:=

RegionPlotΩ,



Out[5]=

Set up the no-slip boundary conditions for the velocities, pressure points and temperature differences between the left and right walls:

In[6]:=

wall=DirichletCondition[{u[t,x,y]0,v[t,x,y]0},True];reference=DirichletCondition[p[t,x,y]0,x0&&y0];temperatures={DirichletCondition[T[t,x,y]1,x0],DirichletCondition[T[t,x,y]0,xlength]};

List together all the boundary conditions. Replace the region's size values:

In[7]:=

bcs={wall,reference,temperatures}/.sizes;bcs//Flatten//Column

Out[7]=

DirichletCondition[{u[t,x,y]0,v[t,x,y]0},True]

DirichletCondition[p[t,x,y]0,x0&&y0]

DirichletCondition[T[t,x,y]1,x0]

DirichletCondition[T[t,x,y]0,x2]

Set up initial conditions such that the system is initially at rest:

In[8]:=

ic={u[0,x,y]0,v[0,x,y]0,p[0,x,y]0,T[0,x,y]0};

Monitor the time-integration progress and the total time it takes to solve the PDE while using a specified mesh spacing and interpolating the velocities "

" and "

" and the temperature "

" to second order and the pressure

to first order:

In[9]:=

Monitor[AbsoluteTiming[{xVel,yVel,pressure,temperature}=NDSolveValue[{op{0,0,0,0},bcs,ic},{u,v,p,T},{x,y}∈Ω,{t,0,300},Method{"PDEDiscretization"{"MethodOfLines","SpatialDiscretization"{"FiniteElement","MeshOptions"{"MaxCellMeasure"0.00125},"InterpolationOrder"{u2,v2,p1,T2}}}},EvaluationMonitor(currentTime=Row[{"t = ",CForm[t]}])];],currentTime]

Out[9]=

{54.7317,Null}

Construct a visualization of the boundary of the region:

In[10]:=

bmesh=HighlightMeshBoundaryDiscretizeRegionΩ,

...

,{Style[1,Black],Style[2,None]},

...



Out[10]=

Visualize the pressure distribution at time

t=150

seconds:

In[11]:=

ShowContourPlotpressure[150,x,y],{x,y}∈pressure["ElementMesh"],

,bmesh,

...



Out[11]=

Visualize the temperature distribution at time

t=150

seconds:

In[12]:=

ShowContourPlottemperature[150,x,y],{x,y}∈temperature["ElementMesh"],

,bmesh,

...



Out[12]=

Visualize the velocity field at time

t=150

seconds:

In[13]:=

Showbmesh,VectorPlot{xVel[150,x,y],yVel[150,x,y]},{x,y}∈xVel["ElementMesh"],

,

...



Out[13]=

Animate the change in temperature and the velocity streamlines:

In[14]:=

frames=Once@TableShowbmesh,ContourPlottemperature[t,x,y],{x,y}∈temperature["ElementMesh"],

,StreamPlot{xVel[t,x,y],yVel[t,x,y]},{x,y}∈xVel["ElementMesh"],

,{t,0,300,2.5};

In[15]:=

FrameListVideo[Downsample[frames,5],FrameRate5]

Out[15]=

Publisher Information

Contributed by: Wolfram Staff

Model Energy Transport Dynamics

Related Symbols

Publisher Information