D2Q9 Heat Transfer

Before start, it is assumed that the reader is already familiar with:

• Basic model creation in the TCLB code, see first tutorial.

• Basics of the LB physics, see obligatory SRT LBM tutorial and supplementary Schan Chen SRT tutorial.

We will cover following topics:

• Heat transfer on D2Q9 lattice using Double Distribution Function (DDF) approach.

Macroscopic Equation

We are going to simulate the heat transfer within a moving fluid.

The motion of the fluid is described by the Navier-Stokes equation: $$\begin{eqnarray} \rho\left( \frac{\partial \textbf{u} }{\partial t} + \textbf{u} \cdot \nabla \textbf{u}\right) &=& -\nabla p + \nabla \cdot (\mu[ \nabla \textbf{u} + (\nabla \textbf{u})^\top]) + \textbf{F}_b \end{eqnarray}$$

The temperature can be recovered from the (simplified) balance of Enthalpy $H=\rho c_p T$:

$$\begin{align} \frac{\partial}{\partial t} (\rho c_p T ) + \nabla \cdot (\boldsymbol{u} \rho c_p T ) &= \nabla \cdot k \nabla T + \dot{q} \label{eq:heat_transfer_eq} \end{align}$$

Where $k$ is the thermal conductivity coefficient.

The SRT-lattice Boltzmann Equation

The continous Boltzmann equation with a forcing term is known as: $$\frac{\partial f}{\partial t} + \xi \cdot \nabla f + \frac{\boldsymbol{F}}{\rho} \cdot \nabla_{\xi} f = \Omega(f)$$

While, the discrete form of the Boltzmann transport equation has been already expressed as:

$$f_i(x+c_i,t+1) =f_i(x,t) + \Omega_f(x,t) + S_i(x,t)$$

Here, we add another distribution function $h$, which will be responsible for the motion of a scalar (like temperature) field.

$$h_i(x+c_i,t+1)=h_i(x,t) + \Omega_h(x,t)$$

As before, the collision operator $\Omega$ is defined using the BGK, a.k.a Single Relaxation Time (SRT) scheme:

$$\begin{aligned} \Omega_f(x,t) &= \frac{1}{\tau_f} \left[f_i^{eq}(x,t) - f_i(x,t) \right] \\ \Omega_h(x,t) &= \frac{1}{\tau_h} \left[h_i^{eq}(x,t) - h_i(x,t) \right] \end{aligned}$$ For the D2Q9 model, we choose nine discrete directions as,

$$\boldsymbol{c} = \left( \begin{array}& 0 & 1 & 0 & -1 & 0 & 1 & -1 & -1 & 1 \\ 0 & 0 & 1 & 0 & -1 & 1 & 1 & -1 & -1 \end{array} \right)$$

Additionally, we define the equilibrium distribution function found by an expansion of a Maxwellian distribution as,

$$\begin{aligned} f_i^{eq}(x,\boldsymbol{u}^{eq},t) &= \omega_i \rho \left( 1 + \frac{\boldsymbol{c}_i \cdot \boldsymbol{u}}{c_s^2} + \frac{(\boldsymbol{c}_i \cdot \boldsymbol{u})^2}{2c_s^4} - \frac{\boldsymbol{u}^2}{2c_s^2} \right) \\ h_i^{eq}(x,\boldsymbol{u}^{eq},t) &= \omega_i \rho \left( 1 + \frac{\boldsymbol{c}_i \cdot \boldsymbol{u}}{c_s^2} + \frac{(\boldsymbol{c}_i \cdot \boldsymbol{u})^2}{2c_s^4} - \frac{\boldsymbol{u}^2}{2c_s^2} \right) \end{aligned}$$

The lattice velocities and weight coefficients are defined as usual:

$$c_s^2 = \frac{1}{3}$$

$$\omega_i=\left\{ \begin{array}{ll} \frac{4}{9} & i=0\\ \frac{1}{9} & i=1,2,3,4\\ \frac{1}{36} & i=5,6,7,8 \end{array} \right.$$

(De)Coupling of Navier-Stokes and Energy equations

Physically, equations shall be coupled:

• Energy $\rightarrow$ NS: equation of state $f(p,\rho,T) = 0$. Usually ideal gas $p(\rho, T) = \rho cs^2 = \rho R T$ is assumed for single phase LBM models.
• NS $\rightarrow$ Energy: kinetic energy + dissipation (viscous heating) and compression work.

In simplified models, the NS equation is decoupled from energy equation. Then, the equation of state has constant temperature $p(\rho, T) = \rho cs^2 = \rho R T_0$ and the sound speed is fixed as $c_s = \sqrt{R T_0} = 1/3$ . As a result, these models are incompressible. To account for thermal advection, the Boussinesq approximation is usually employed:

$$\begin{aligned} \rho(T) &\approx \rho_0 (1- \alpha_V (T-T_0)) \\ \boldsymbol{F}_{bouyancy} &= [\rho(T) - \rho_0] \boldsymbol{g} = - \boldsymbol{g} \rho_0 \alpha_V (T-T_0) \end{aligned}$$

In present tutorial we will skip the viscous dissipation and compression work and treat the energy equation as a transport of a simple passive scalar (temperature), thus there is no source term in discrete transport equation for $h$.

In this tutorial, to incorporate the force into the fluid-LB scheme we will use the He's scheme:

$$f_i(x+c_i,t+1)=f_i(x,t) + \Omega_f(x,t) + S_i(x,t)$$

Accoring to the He's idea ¹ we can approximate the gradient of unknown particle distribution function $\nabla_{\xi} f$ in continous Boltzmann equation by its dominant, equilibrium part $\nabla_{\xi} f^{eq}$

$$\frac{\boldsymbol{F}}{\rho} \cdot \nabla_{\xi} f \approx \frac{\boldsymbol{F}}{\rho} \cdot \nabla_{\xi} f^{eq} = \frac{\boldsymbol{F}}{\rho} \cdot f^{eq} \frac{\boldsymbol{c} - \boldsymbol{u}}{c_s^2}$$

Finally, the discrete forcing term can be expressed as: $$S_i(x,t) = \left(1 - \frac{\Delta t}{2 \tau_f} \right) \frac{\boldsymbol{F}}{\rho} \cdot f_i^{eq} \frac{\boldsymbol{c}_i - \boldsymbol{u}}{c_s^2}$$

The macroscopic velocity is defined as: $$\boldsymbol{u} = \boldsymbol{u}^{eq}= \frac{1}{\rho} \sum_i \boldsymbol{c}_i f_i + \frac{\boldsymbol{F}\Delta t}{2 \rho}$$

The $\boldsymbol{u}^{eq}$ is the velocity entering $f_i^{eq}$ and $h_i^{eq}$.

The origin of the $\frac{\Delta t}{2 \tau_f}$ and $\frac{1}{2}$ terms in the last two equations comes from the second order (trapezoidal) time discretization of the Bolztmann equation. We left the proof as an exercice for the reader ;)

Model Creation in TCLB

Dynamics.R

First, we have to define the variables accessible within the nodes, outputs from the simulation and batch (input) settings. These informations are stored in Dynamics.R file.

Densities

As in the first tutorial, we have to set up a folder named d2q9_srt_heat in ~TCLB/models/tutorial/ and create the generic file structure (conf.mk, Dynamics.c.Rt, Dynamics.R).

To start off the model, we look to add the nine distribution functions for the nine discrete velocities. This is implemented in the Dynamics.R file, but instead of adding each as its own field we look to stream the functions in the process of calling them at each node. To do this we instead add them as Densities:

source("lib/lattice.R")  # import R symbols

#  Density - table of variables of LB Node to stream
#  name - variable name to stream
#  dx,dy,dz - direction of streaming
#  comment - additional comment

# by hand: AddDensity( name="f[0]", dx= 0, dy= 0, group="f")
AddDensity(
    # name = paste("f",1:9-1,sep=""),  # without brackets
    name = paste("f[",1:9-1,"]",sep=""),
    dx   = d2q9[,1],
    dy   = d2q9[,2],
    comment=paste("flow LB density F",1:9-1),
    group="f"
)

AddDensity(
    name = paste("h[",1:9-1,"]",sep=""),
    dx   = d2q9[,1],
    dy   = d2q9[,2],
    comment=paste("heat LB density H",1:9-1),
    group="h"
)

Notice that the dx and dy coordinates correspond with the c matrix previously given. At this stage, it helps to assess what other values will be needed to initialise and run the LBM. To start off the method, a fluid density and initial velocities must be specified, then to perform the collision operation the relaxation time and the magnitude of the applied body force need to be defined .

Quantities

We want to be able to interrogate the macroscopic fluid velocity, density and temperature so these must be added as Quantities. For this, we define (in Dynamics.R):

#   Outputs:
AddQuantity(name="Rho")
AddQuantity(name="T")
AddQuantity(name="U", vector=T)

Settings

Here, we define simulation settings (like properties of fluid). Once the simulation is launched, the solver will look for them in the batch.xml file.

#   Inputs: Flow Properties
AddSetting(name="VelocityX", default=0.0, comment='inlet/outlet/init x-velocity component', zonal=TRUE)
AddSetting(name="GravitationY", default=0.0, comment='applied rho*GravitationY')
AddSetting(name="GravitationX", default=0.0, comment='applied rho*GravitationX')
AddSetting(name="omega_nu", comment='inverse of viscous relaxation time', default=1.0)
AddSetting(name="nu", omega_nu='1.0/(3*nu+0.5)',  default=0.16666666,  comment='kinematic viscosity')
AddSetting(name="omega_bulk", comment='inverse of bulk relaxation time', default=1.0)
AddSetting(name="bulk_visc", omega_bulk='1.0/(3*bulk_visc+0.5)',  comment='bulk viscosity')


#   Inputs: Initial Thermal Properties
AddSetting(name="InitTemperature", default=0, comment='Initial/Inflow temperature distribution', zonal=T)

#   Inputs: Fluid Thermal Properties
AddSetting(name="omega_k", default=1.0 , comment='inverse of thermal relaxation time')
AddSetting(name="k", omega_k='1.0/(3*k+0.5)', default=0.16666666, comment='thermal conductivity of fluid (W/(m·K))')
AddSetting(name="cp", , default=1.0, comment='specific heat capacity at constant pressure of fluid (J/(kg·K))')
AddSetting(name="BoussinesqCoeff", default=0.0, comment='BoussinesqCoeff=rho_0*thermal_exp_coeff')

BC - Node types

#   Boundary things
AddNodeType(name="HeaterDirichletTemperature", group="ADDITIONALS_HEAT")
AddNodeType(name="HeaterNeumannHeatFlux", group="ADDITIONALS_HEAT")

Dynamics.c.Rt

From the above, we have all the variables we need to implement the LBM. The first step is to initialise the lattice over the required domain, this is incorporated as part of the dynamics that are occurring in the model and must be incorporated into the Dynamics.c file within the Init() function.

Initialization

CudaDeviceFunction void Init() {
    real_t rho = 1.0;
    real_t T = InitTemperature;
    vector_t u; u.x = VelocityX; u.y = 0; u.z = 0;

    SetEquilibrium(f,rho,u);
    SetEquilibrium(h,rho*T,u);
}

Notice here that we are calling the SetEquilibrium() twice to calculate the equilibrium distribution for both $f_i^{eq}$ and $h_i^{eq}$.

CudaDeviceFunction void SetEquilibrium(real_t x_eq[9], real_t rhoX, vector_t u){
    // modify the SetEquilibrium() function from previous tutorials so that it can be used for both f and h populations (distributions).
}

Additionally, in all Dynamics.c files the Color() function is required. (even with ./configure --disable-graphics).

Macroscopic Quantities

CudaDeviceFunction real_t getRho() {
// This function defines the macroscopic density at the current node.
    return f[8]+f[7]+f[6]+f[5]+f[4]+f[3]+f[2]+f[1]+f[0];
}

CudaDeviceFunction vector_t getRawU() {
// This function defines the macroscopic velocity at the current node.
    real_t rho;
    // implement rho as sum(f[i]).
    // rho = ... ;
    vector_t u;
    // the sum(f[i]*ex[i]) is interpreted as x-component of momentum, i.e., rho*u_x.
    // implement u.x as sum(f[i]*ex[i]) divided by rho. Similarly for u.y. See previous tutorials.
    // u.x = ... ;
    // u.y = ... ;
    u.z = 0;
    return u;
}

CudaDeviceFunction real_t getT(){
    real_t T;
    // implement T as sum(h[i]) divided by rho.
    return T; 
}

CudaDeviceFunction vector_t getU()
{
    real_t localTemperature = getT();
    vector_t u = getRawU();
    real_t m00 = getRho();
    vector_t Force = getForce(localTemperature, m00);
    u.x += Force.x/(2*m00);
    u.y += Force.y/(2*m00);
    u.z = 0;
    return u;
}

CudaDeviceFunction float2 Color() {
    float2 ret;
    ret.x = getT();
    ret.y = ret.x;
    return ret;
}

Main loop

CudaDeviceFunction void Run() {
// This defines the dynamics that we run at each node in the domain.
    switch (NodeType & NODE_BOUNDARY) {
        case NODE_Wall:
            BounceBack();
            break;
    }

    switch (NodeType & NODE_COLLISION) {
        case NODE_BGK:
            CollisionBGK();
            break;
        case NODE_MRT:
        //case NODE_CM:
            //CollisionCM();
            break;
        }

    if ((NodeType & NODE_ADDITIONALS_HEAT) == NODE_HeaterDirichletTemperature) {
         HeatDirichletEquilibriumScheme();
    }
}

Boussinesq approximation

CudaDeviceFunction vector_t getForce(real_t localTemp, real_t rho)
{
    // Boussinesq approximation
    // rho(T) ~ rho_0*(1-thermal_exp_coeff*(Temp-Temp_0))
    // F_b = (rho(T) - rho_0)*Grav_Y = -Grav_Y*rho_0*thermal_exp_coeff*(Temp-Temp_0)
    // see chapter 8.4.1, eq 8.44, p313 from 'The Lattice Boltzmann Method: Principles and Practice'
    // by T. Krüger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva, E.M. Viggen
    real_t refTemperature = 0;

    // implement the BoussinesqForce
    // real_t BoussinesqForce = ...; // BoussinesqCoeff=rho_0*thermal_exp_coeff

    vector_t Force;
    Force.x = GravitationX*rho;
    Force.y = GravitationY*rho+BoussinesqForce;
    Force.z = 0;

    return Force;
}

SRT Collision

CudaDeviceFunction void CollisionBGK() {
// Here we perform a single relaxation time collision operation.

    real_t localTemperature = getT();
    real_t m00 = getRho();

    vector_t Force = getForce(localTemperature, m00);
    vector_t u = getRawU();
    u.x += Force.x/(2*m00);
    u.y += Force.y/(2*m00);

    real_t f_eq[9]; real_t h_eq[9];
    // calcluate f_eq and h_eq using the SetEquilibrium() function.

    real_t S[9];
    for (int i=0; i< 9; i++) {
            // Using the formulas from the theroretical section, calculate the forcing term.
            // tclb hint: use d2q9_ex[i] and d2q9_ey[i] to perform the vector operations.
            // S[i] = ... ;

            // Calculate collision for fluid populations
            // f[i] =  ... ;    
    }


    for (int i=0; i< 9; i++) {
            // Calculate collision for temperature populations
            // h[i] = ...;
    }
}

Boundary Conditions

Temperature: Dirichlet BC

CudaDeviceFunction void HeatDirichletEquilibriumScheme()
{
    // equilibrium scheme for BC - don't care and impose rho*Teq
    // see chapter 5.3.4.2, eq 5.34, p191 from 'The Lattice Boltzmann Method: Principles and Practice'
    // by T. Krüger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva, E.M. Viggen

    real_t T = InitTemperature;
    real_t m00 = getRho();
    vector_t u = getU();

    real_t d_rhoT = m00*T;
    real_t h_eq[9];
    SetEquilibrium(h_eq, d_rhoT, u);

    for (int i = 0; i < 9; i++) {
        h[i] = h_eq[i];
    }
}

Setting up a Simulation

Create a configuration file, my_heat_cube.xml in example/

<?xml version="1.0"?>
<CLBConfig version="2.0" output="output/">
      <Geometry nx="128" ny="128" >
      <BGK><Box /></BGK>

        <None name="heated_region">
            <!-- <Sphere dx="80" nx="10" dy="25" ny="10" />  -->
        </None>        

        <HeaterDirichletTemperature  name="dirichlet_hot"> 
          <Box nx="1" dx="-2"/>
        </HeaterDirichletTemperature>

        <HeaterDirichletTemperature  name="dirichlet_cold"> 
          <Box dx="1" nx="1" />
        </HeaterDirichletTemperature>

        <Wall mask="ALL" name="border" >
                <Box nx="1"/>
                <Box dx="-1"/>
                <Box ny="1"/>
                <Box dy="-1"/>
        </Wall>
  </Geometry>
    <Model>
    <Param name="k" value="0.1666666" />
    <Param name="cp" value="1" />

    <Param name="InitTemperature" value="0"/>
    <Param name="InitTemperature" value="100" zone="heated_region"/>
    <Param name="InitTemperature" value="1" zone="dirichlet_hot"/>
    <Param name="InitTemperature" value="-1" zone="dirichlet_cold"/>

    <Param name="GravitationY" value="-1.0E-05" />
    <Param name="BoussinesqCoeff" value="10"/>

    </Model>

<Solve Iterations="10" output="output/"> <VTK Iterations="1"/> </Solve>
<Failcheck Iterations="500" nx="128" ny="128" />
<Log Iterations="100"/>
<Solve Iterations="10000" output="output/"> <VTK Iterations="500"/> </Solve>
</CLBConfig>

Question

What is the difference between specifing the BC for temperature before/after the <Wall mask="ALL" name="border" >... </Wall> region? Is the following description correct? <HeaterDirichletTemperature name="dirichlet_hot"> <Box dx="-1"/> </HeaterDirichletTemperature> <HeaterDirichletTemperature name="dirichlet_cold"> <Box nx="1" /> </HeaterDirichletTemperature>

Compilation

With the model and set-up files created, we can now look to make and run d2q9srt. To do this first enter the TCLB directory and call the run file along with the input file location, after this, the analysis can be performed in Paraview:

make d2q9_srt_heat
CLB/d2q9_srt_heat/main my_heat_cube.xml
paraview output/my_heat_cube_VTK_P00_..pvti

In case the calculations are run on CPU, an mpi run can also be initiated e.g:

mpirun -np 4 CLB/d2q9srt/main my_heat_cube.xml

References

Tip

A good starting point explaining the 'zoology' of various LBM models is a modern (2017) book 'The Lattice Boltzmann Method: Principles and Practice' written by T. Krüger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva, E.M. Viggen .

---

1. Xiaoyi He, Xiaowen Shan, and Gary D. Doolen. 'Discrete Boltzmann equation model for nonideal gases' in Physical Review E - Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics (1998). ↩