In some cases, it might be useful to advect scalar fields with the velocity obtained from the solution of the Navier-Stokes equation. For example, in study of mass transfer or heat transfer problems where getting analytical expression for advection velocity is not possible, the transport (advection-diffusion) equation needs to be solved along with the Navier-Stokes equation to get the scalar concentration or temperature distribution in the flow field.
In the input file, the extra field variables that are being advected need to be defined after the variables representing the velocity components. The pressure needs to be at the end of the list. For example, for a 2D simulation the expansions and variables would be defined as
1<EXPANSIONS> 2 <E COMPOSITE="C[0]" NUMMODES="5" FIELDS="u,v,c1,c2,p" TYPE="MODIFIED" /> 3</EXPANSIONS> 4 5<VARIABLES> 6 <V ID="0"> u </V> 7 <V ID="1"> v </V> 8 <V ID="2"> c1 </V> 9 <V ID="3"> c2 </V> 10 <V ID="4"> p </V> 11</VARIABLES>
where u, v are the velocity components, c1 and c2 are extra fields that are being advected and p is the pressure.
In addition, diffusion coefficients for each extra variable can be specified by adding a function
DiffusionCoefficient
1<FUNCTION NAME="DiffusionCoefficient"> 2 <E VAR="c1" VALUE="0.1" /> 3 <E VAR="c2" VALUE="0.01" /> 4</FUNCTION>
Boundary conditions for the extra fields are set up in the same way as the velocity and pressure
1<BOUNDARYCONDITIONS> 2 <REGION REF="0"> 3 <D VAR="u" VALUE="0" /> 4 <D VAR="v" VALUE="0" /> 5 <D VAR="c1" VALUE="1" /> 6 <D VAR="c2" VALUE="0" /> 7 <N VAR="p" USERDEFINEDTYPE="H" VALUE="0" /> 8 </REGION> 9</BOUNDARYCONDITIONS>
It should be noted that if the diffusion coefficient is too small, the transport equation becomes advection dominated. This leads to small grid spacing required to resolve all physical scales associated with the transport equation (the ratio of resolution required for transport to Navier Stokes equation scales with Sc3∕4, where Sc is the Schmidt number = kinematic viscosity/diffusion coefficient). Hence, small diffusion coefficient might lead to spurious oscillations if the mesh spacing is not small enough.
Next, we consider an example of an active scalar in the case of Rayleigh-Bénard convection (RBC) where the temperature (scalar) field is coupled with the momentum equations. RBC describes the motion of a layer of fluid sandwiched between two infinite parallel plates, separated by a distance, d, heated from the bottom and cooled from the top - a canonical fluid configuration of natural convection. The governing equations of RBC can be described by the non-dimensional Navier-Stokes equations with Boussinesq approximations given as,
with the following boundary conditions,
(11.32) |
where the Ra,Pr,θ,u,p refers to the Rayleigh number, Prandtl number, non-dimensional temperature and velocities respectively. To setup a simulation for a 2D simulation of RBC, the expansions and variables are prescribed as,
1<EXPANSIONS> 2 <E COMPOSITE="C[0]" NUMMODES="5" FIELDS="u,v,theta,p" TYPE="MODIFIED" /> 3</EXPANSIONS> 4 5<VARIABLES> 6 <V ID="0"> u </V> 7 <V ID="1"> v </V> 8 <V ID="2"> theta </V> 9 <V ID="3"> p </V> 10</VARIABLES>
The boundary conditions are set up in the same way as listing 11.9. Note that the boundary
conditions on the left/right edges are periodic, and its setup is shown in section 3.4.6.3.
Next, the diffusion coefficient for the momentum equations (equation 11.31a) is
prescribed directly by setting Kinvis
(within PARAMETERS
) to the Prandtl number, i.e
Kinvis = Pr
. The diffusion coefficient for the non-dimensional temperature equation
(equation 11.31b) in this case is 1, which is prescribed in DiffusionCoefficient
as,
1<FUNCTION NAME="DiffusionCoefficient"> 2 <E VAR="theta" VALUE="1" /> 3</FUNCTION>
The RaPrθj term refers to the temperature (scalar) field acting as a body force on the
v-momentum equations (equation 11.31a), which is prescribed using the BodyForce
function,
1<FUNCTION NAME="BodyForce"> 2 <E VAR="u" VALUE="0" EVARS="u v theta" /> 3 <E VAR="v" VALUE="Ra * Pr * theta" EVARS="u v theta" /> 4 <E VAR="theta" VALUE="0" EVARS="u v theta" /> 5</FUNCTION>
The onset of convection occurs above the critical Rayleigh number, Rac = 1708. The solution for Ra = 2000,Pr = 1, consisting of two pairs of counter-rotating convection rolls are shown below. Note that the wavelength of a pair of rolls is close to the critical wavelength i.e λc ≈ 2d.