Heat conduction PDE

The classical heat conduction equation is a partial differential equation (PDE) that describes the conduction of heat through a material or system. It's important to note that the heat conduction equation can be modified or extended to include additional terms, such as the convection term, which accounts for heat transfer due to domain movements. This results in the heat convection-diffusion equation, which can be seen a generalization of the heat conduction equation. With this equation both, conduction and convection can be modeled. In the following sections the derivation how to obtain these two strong forms will be described.

Governing equations

The derivation of the classical heat conduction equation can be derived by considering Fourier's law (\ref{eq:fourierslaw}), calorics law (\ref{eq:caloriclaw}) and the first law of thermodynamics (\ref{eq:firstlawthermo}).

$$\begin{equation} \label{eq:fourierslaw} \dot{\textbf{q}} = -k(T) \nabla T \end{equation}$$

$$\begin{equation} \label{eq:caloriclaw} u = c_p \cdot T \end{equation}$$

$$\begin{equation} \label{eq:firstlawthermo} \frac{dU}{dt} = \dot{Q}(t) + P(t) \end{equation}$$

Mathematically this PDE is equivalent to the ElectricFlowPDE, both of them are Poisson problems.

The heat conduction problem can be stated as:

\begin{equation} \label{eq:heatconductionpde} \rho c_{\textrm{m}} \frac{\partial T}{\partial t} -\nabla \cdot ({k(T)\nabla{T}}) = \dot{{q}}, \end{equation} where $\rho$ is the mass density, $c_{\textrm{m}}$ the specific heat capacity (per unit mass), $k$ the thermal conductivity and $\dot{q}$ the volumetric heat source.

For moving domains an additional term with a velocity $\textbf{v}$ enters to the heat conduction PDE (\ref{eq:heatconductionpde}), which is also known as the convection term. In the next section it is shown how this term is obtained.

Heat convection diffusion Equation

The starting point to derive the convection diffusion equation builds the transport theorem due to Reynolds which is a 3D generalization of the leibniz integration rule. The theorem reads as follows:

$$\begin{equation} \label{eq:reynolstheorem} \frac{d}{dt} \int_{\Omega} \Phi d\Omega = \int_{\Omega} \frac{\partial}{\partial t}\Phi d\Omega + \int_{\partial \Omega} \Phi \textbf{v} \cdot \textbf{n} d\Gamma \end{equation}$$

Setting $\Phi = \rho u$, leads to the following:

$$\begin{equation} \label{eq:2} \frac{d}{dt} \int_{\Omega} \rho u d\Omega = \int_{\Omega} \frac{\partial}{\partial t}\rho u d\Omega + \int_{\partial \Omega} \rho u \textbf{v} \cdot \textbf{n} d\Gamma \end{equation}$$

The left hand side of equation \ref{eq:2} can be transformed with the 1st law of thermodynamics \ref{eq:firstlawthermo} and using fouriers law (\ref{eq:fourierslaw}). The right hand side of the equation can be reformulated by using calorics law (\ref{eq:caloriclaw}).

$$\begin{equation} \label{eq:3} -\int_{\partial \Omega} k(T) \nabla T \cdot \textbf{n} d\Gamma + \int_{\Omega} \dot{q}_d d\Omega = \int_{\Omega} \frac{\partial}{\partial t}\rho c_m T d\Omega + \int_{\partial \Omega} \rho c_m T \textbf{v} \cdot \textbf{n} d\Gamma \end{equation}$$

using Gauß integral theorem results in the final equation for the convective heat equation.

$$\begin{equation} \label{eq:4} \frac{\partial}{\partial t}\rho c_m T + \nabla \cdot (\rho c_m T \textbf{v})- \nabla \cdot (k(T) \nabla T) = \dot{q}_d \end{equation}$$

if $k(T) = k,\,\,$ $\rho(\textbf{x})=\rho$ and $c_m(p,T) = c_m$, the equation above can be simplified to:

$$\begin{equation} \label{eq:heatconvectiondiffusionpde} \frac{\partial T}{\partial t} + \nabla \cdot (T \textbf{v})- c^2\nabla \cdot\nabla T = \frac{1}{\rho c_m}\dot{q}_d \end{equation}$$

with $c^2$ as the diffusion coefficient.

Boundary conditions

As described above, there are two types of BC's, Neumann and Dirichlet. However, since openCFS is physics-based, there are several variations and combinations of these. They can be defined in the simulation-xml file in the BC-section

        <bcsAndLoads>
            <temperature name="" value=""/>
            <heatFlux name="" value=""/>
            <heatSource name="" value=""/>
            <heatSourceDensity>
<!--                 via: -->
<!--                    * coupled to another PDE -->
<!--                    * input from external simulation -->
<!--                    * scattered data input from csv file -->
<!--                    * sequence step (forward coupling from another PDE) -->
            </heatSourceDensity>
            <heatTransport name=""  volumeRegion="" heatTransferCoefficient="" bulkTemperature="20.0"/>
            <periodic slave="" master="" quantity=""/>
            <elecPowerDensity>
<!--                 via: -->
<!--                    * coupling to other PDE -->
            </elecPowerDensity>
        </bcsAndLoads>

• Temperature <temperature name="" value=""/>: Inhomogeneous Dirichlet value $T_{\textrm{e}}$
• Heat Flux <heatFlux name="" value=""/>: Prescribes the heat flux density (in $\frac{W}{m^2}$) across a certain boundary $\dot{\mathbf{q}}\cdot \mathbf{n}$
• Heat Source Density <heatSourceDensity/>: Prescribes the volumetric heat source densities $\dot{q}$
• Heat Transport <heatTransport name="" volumeRegion="" heatTransferCoefficient="" bulkTemperature="20.0"/>: Describes the heat transport between two media, e.g. iron-air; name defines the surface through which the heat is transported, volumeRegion defines the positive direction for the heat flux, e.g. we define the direcion iron->air as positive, then the volumeRegion has to be the iron-region. heatTransferCoefficient defines the transfer coefficient and bulkTemperature defines the ambient temperature of, e.g., the surrounding air domain
• Electric Power Density <elecPowerDensity>: Defines the electric power density as a volume source from a coupled computation as a right-hand-side quantity ($\dot{q}$)

Material

The material parameters $\rho$, $c_{\textrm{m}}$ and $k$ are defined in the material-xml file, e.g. for iron

  <material name="iron">
    <heatConduction>
      <density>
        <linear>
          <real> 7870 </real>
        </linear>
      </density>
      <heatCapacity>
        <linear>
          <real>  444 </real>
        </linear>
      </heatCapacity>
      <heatConductivity>
        <linear>
          <isotropic>
            <real> 80.4 </real>
          </isotropic>
        </linear>
      </heatConductivity>
    </heatConduction>
  </material>

Convection term and stabilisation

The original heat conduction problem can be modified by the convection term. Physically, the convection term accounts for the transport of matter (convection). The heat convection-diffusion equation is

\begin{equation} \rho c_{\textrm{m}} \frac{\partial T}{\partial t} - \nabla \cdot ({k\nabla{T}} - \rho c_{\textrm{m}}\mathbf{v} T) = \dot{{q}}, \end{equation} where $\mathbf{v}$ is the velocity field.

One can define the convective term in a certain region by the velocityId and the corresponding definition of the velocity field in the <velocityList>. Make sure that the velocity field you specify is divergence-free.

    <regionList>
      <region name="V_body"  velocityId="velocity_body"/>
    </regionList>
    <velocityList>
      <velocity name="velocity_body">
       <comp dof="y" value="10"/>
      </velocity>
    </velocityList>

The solution of a convection-diffusion equation is unstable because of the numerical discretisation. When the Péclet number exceeds a critical value ($\mathrm{Pe}>1$), the spurious oscillations result in space. The Péclet number is defined as follows

$$\mathrm{Pe} = \frac{|\mathbf{v}| l^e}{2 k},$$

where $l^e$ denotes the length of an element in the velocity direction and $k$ is a material parameter for a heat conduction problem.

We have multiple options to tackle this instability. In the XML model, one adds a tag stabilisation= to the <velocity/> element.

  <velocityList>
      <velocity name="vel" stabilisation="ArtificialDiffusion" >
          <comp dof="x" value="1" />
          <comp dof="y" value="1" />
          <comp dof="z" value="0" />
      </velocity>
  </velocityList>

Artificial diffusion

Artificial diffusion stabilisation is one of the simplest options. We modify the diffusion part of the equation, so we get an extra term in stiffness matrix in week formulation. Such an addition forces Pe not to exceed 1. The method is referred to as inconsistent, because it adds a certain amount of diffusion independently of how close the numerical solution is to the exact solution.

$$\begin{equation} K_{ji}^e = \int_{\Omega} \rho c_p \mathbf{v} \cdot \nabla (N_j) (N_i + \color{red}{\tau^e \mathbf{v} \cdot \nabla N_i}) + \lambda \nabla(N_j) \cdot \nabla(N_i), \end{equation}$$

$$\begin{gather*} \tau^e = \frac{\alpha l^e}{2 |\mathbf{v}|}, \\ l^e = \frac{2}{\Sigma_n |\frac{\mathbf{v}}{|\mathbf{v}|} \cdot \nabla{w_n}|}, \\ \alpha = \coth{\gamma} - \frac{1}{\gamma},\\ \gamma = \mathrm{Pe} = \frac{|\mathbf{v}| l^e}{2 k},\\ k = \frac{\lambda}{\rho c_p}, \\ \end{gather*}$$

where $\tau^e$ is a stabilisation parameter.

SUPG

The SUPG stabilisation is still in the process of implementation in openCFS.

Another approach to stabilise a solution is Streamline Upwind / Petrov-Galerkin method. This is a consistant method, which means it gives less numerical diffusion the closer the numerical solution comes to the exact solution. In other words, no diffusion is added in the regions where the mesh is fine enough.

The idea is to change the diffusion in the streamline direciton. For this purpose we modify the test function. The test function is \begin{equation} \hat{T}' = T' + P, \end{equation} where $P$ is a discontinuous streamline upwind contribution

$$P = \tau^e \mathbf{v} \cdot \nabla T',$$

where $\tau^e$ is the same as in Artificial diffusion stabilisation.

The test function in FEM space for 1D case is \begin{equation} \hat{T}' = N_i + \frac{\alpha l^e}{2} \frac{d N_i}{d x} , \end{equation}

The week formulation multiplied by a test function

$$\begin{equation} \rho c_p \int_{\Omega} \frac{\partial T}{\partial t} \hat{T}' + \rho c_p \int_{\Omega} \mathbf{v} \cdot \nabla ( T) \hat{T}' - \lambda \int_{\Omega} \nabla \cdot (\nabla(T)) \hat{T}' = \int_{\Omega} \dot{q}_d \hat{T}'. \end{equation}$$

From now on we use only elementwise integration as the test function is discontinuous $\Omega = \Omega^e $.
Integrate by parts and insert BCs

$$\begin{equation} \rho c_p \int_{\Omega} \frac{\partial T}{\partial t} \hat{T}'+\rho c_p \int_{\Omega} \mathbf{v} \cdot \nabla ( T) \hat{T}'+\lambda \int_{\Omega} \nabla(T) \cdot \nabla(\hat{T}') - \int_{\delta \Omega} \lambda \nabla(T) \cdot \mathbf n \hat{T}'=\int{\Omega} \dot{q}d \hat{T}' \end{equation}$$

The element stiffness matrix

$$\begin{equation*} K_{ji}^e = \int_{\Omega} \rho c_p \mathbf{v} \cdot \nabla (N_j) (N_i) + \color{red}{\rho c_p \tau^e \mathbf{v} \cdot \nabla (N_j) \mathbf{v} \cdot \nabla N_i} + \lambda \nabla(N_j) \cdot \nabla(N_i) + \color{red}{\lambda \nabla(N_j) \cdot \tau^e \mathbf{v} \nabla ( \cdot \nabla N_i)} \end{equation*}$$

The rhs

$$\begin{equation*} F_i^e = \int_{\Omega} \dot{q}_d (N_i + \color{red}{\tau^e \mathbf{v} \cdot \nabla N_i}) = \int{\Omega} \dot{q}_d N_i + \color{red}{\tau^e \dot{q}_d \mathbf{v} \cdot \nabla N_i} \end{equation*}$$

The rhs from BCs

$$\begin{equation*} F_i^{e\, BC} = \int{\delta \Omega}f (N_i + \color{red}{\tau^e \mathbf{v} \cdot \nabla N_i}) \end{equation*}$$

The mass matrix

$$\begin{equation*} M_{ji}^e = \rho c_p \int_{\Omega} N_j (N_i + \color{red}{\frac{(\coth{\frac{|\mathbf{v}| l^e}{2 k}} - \frac{1}{\frac{|\mathbf{v}| l^e}{2 k}}) l^e}{2 |\mathbf{v}|} \mathbf{v} \cdot \nabla N_i}) = \rho c_p \int_{\Omega} N_j N_i + \color{red}{\tau^e N_j \mathbf{v} \cdot \nabla N_i} \end{equation*}$$

Good pdf with the formulation in section 6.9.1.

Analysis Types

Since in general we are dealing with a time-dependent PDE, we can distinguish three analysis types:

• Static: $\frac{\partial T}{\partial t}=0$, in this case, the PDE is a purely elliptic PDE (Poisson problem)

      <analysis>
         <static/>
      </analysis>
  

• Transient: $\frac{\partial T}{\partial t}\neq 0$

    <analysis>
      <transient>
        <numSteps>10</numSteps>
        <deltaT>0.1</deltaT>
      </transient>
    </analysis>

• Harmonic: $\partial / \partial t (\cdot) = j \omega (\cdot)$

  <analysis>
      <harmonic>
          <numFreq>10</numFreq>
          <startFreq>10</startFreq>
          <stopFreq>3000</stopFreq>
          <sampling>linear</sampling>
      </harmonic>
  </analysis>
  

Postprocessing Results

Node Results

• Temperature (primary solution)

  <nodeResult type="heatTemperature">
  

• Temperature Time-Derivative

  <nodeResult type="heatTemperatureD1">
  

Element Results

• Heat Flux Density \begin{equation} {\mathbf q} = - k \cdot \nabla T. \end{equation}

  <elemResult type="heatFluxDensity">
  

Surface Element Results

• Heat Flux Density Surface Normal \begin{equation} \dot{q}_{\textrm{n}} = \dot{\mathbf q} \cdot \mathbf n \end{equation}

  <surfElemeResult type="heatFluxIntensity">
  

Surface Region Results

• Heat Flux \begin{equation} \dot{Q} = \int_{{A}} \dot{q}_{\textrm{n}} \, ds \end{equation}

  <surfRegionResult type="heatFlux">
  

Tutorial

Analysis Workflow