Heat Transfer Model

The heat transfer model is a specific implementation of the Model interface dedicated to handle the dynamic heat equation.

Theory

The strong form of the dynamic heat equation can be expressed as

\[\rho c_v \dot{T} + \nabla \cdot \vec{\kappa} \nabla T = b\]

with \(T\) the scalar temperature field, \(c_v\) the specific heat capacity, \(\rho\) the mass density, \(\mat{\kappa}\) the conductivity tensor, and \(b\) the heat generation per unit of volume. The discretized weak form with a finite number of elements is

\[\forall i \quad \sum_j \left( \int_\Omega \rho c_v N_j N_i d\Omega \right) \dot{T}_j - \sum_j \left( \int_\Omega \vec{\kappa} \nabla N_j \nabla N_i d\Omega \right) T_j = - \int_{\Gamma} N_i \vec{q} \cdot \vec{n} d\Gamma + \int_\Omega b N_i d\Omega\]

with \(i\) and \(j\) the node indices, \(\vec{n}\) the normal field to the surface \(\Gamma = \partial \Omega\). To simplify, we can define the capacity and the conductivity matrices as

\[C_{ij} = \int_\Omega \rho c_v N_j N_i d\Omega \qquad \textrm{and} \qquad K_{ij} = - \int_\Omega \vec{\kappa} \nabla N_j \nabla N_i d\Omega\]

and the system to solve can be written

\[\mat{C} \cdot \vec{\dot{T}} = \vec{Q}^{\text{ext}} -\mat{K} \cdot \vec{T}~,\]

with \(\vec{Q}^{\text{ext}}\) the consistent heat generated.

Using the Heat Transfer Model

A material file name has to be provided during initialization. Currently, the HeatTransferModel object uses dynamic analysis with an explicit time integration scheme. It can simply be created like this

HeatTransferModel model(mesh, spatial_dimension);

while an existing mesh has been used (see ref{sect:common:mesh}). Then the model object can be initialized with:

model.initFull()

This function will load the material properties, and allocate / initialize the nodes and element Arrays More precisely, the heat transfer model contains 4 Arrays:

• temperature contains the nodal temperature \(T\) (zero by default after the initialization).

• temperature_rate contains the variations of temperature \(\dot{T}\) (zero by default after the initialization).

• blocked_dofs contains a Boolean value for each degree of freedom specifying whether the degree is blocked or not. A Dirichlet boundary condition (\(T_d\)) can be prescribed by setting the blocked_dofs value of a degree of freedom to true. The temperature and the temperature_rate are computed for all degrees of freedom where the blocked_dofs value is set to false. For the remaining degrees of freedom, the imposed values (zero by default after initialization) are kept.

• external_heat_rate contains the external heat generations. \(\vec{Q^{ext}}\) on the nodes.

• internal_heat_rate contains the internal heat generations. \(\vec{Q^{int}} = -\mat{K} \cdot \vec{T}\) on the nodes.

Only a single material can be specified on the domain. A material text file (e.g. material.dat) provides the material properties as follows:

model heat_transfer_model [
  capacity = %\emph{XXX}%
  density = %\emph{XXX}%
  conductivity = [%\emph{XXX}% ... %\emph{XXX}%]
]

where the capacity and density are scalars, and the conductivity is specified as a \(3\times 3\) tensor.

Explicit Dynamic

The explicit time integration scheme in Akantu uses a lumped capacity matrix \(\mat{C}\) (reducing the computational cost, see Chapter Solid Mechanics Model). This matrix is assembled by distributing the capacity of each element onto its nodes. Therefore, the resulting \(\mat{C}\) is a diagonal matrix stored in the capacity Array of the model.

model.assembleCapacityLumped();

Note

Currently, only the explicit time integration with lumped capacity matrix is implemented within Akantu.

The explicit integration scheme is Forward Euler [Cur92].

• Predictor: \(\vec{T}_{n+1} = \vec{T}_{n} + \Delta t \dot{\vec{T}}_{n}\)

• Update residual: \(\vec{R}_{n+1} = \left( \vec{Q^{ext}_{n+1}} - \vec{K}\vec{T}_{n+1} \right)\)

• Corrector : \(\dot{\vec{T}}_{n+1} = \mat{C}^{-1} \vec{R}_{n+1}\)

The explicit integration scheme is conditionally stable. The time step has to be smaller than the stable time step, and it can be obtained in Akantu as follows:

time_step = model.getStableTimeStep();

The stable time step is defined as:

(8)\[\Delta t_{\st{crit}} = 2 \Delta x^2 \frac{\rho c_v}{\mid\mid \mat{\kappa} \mid\mid^\infty}\]

where \(\Delta x\) is the characteristic length (e.g the in-radius in the case of linear triangle element), \(\rho\) is the density, \(\mat{\kappa}\) is the conductivity tensor, and \(c_v\) is the specific heat capacity. It is necessary to impose a time step which is smaller than the stable time step, for instance, by multiplying the stable time step by a safety factor smaller than one.

const Real safety_time_factor = 0.1;
Real applied_time_step = time_step * safety_time_factor;
model.setTimeStep(applied_time_step);

The following loop allows, for each time step, to update the temperature, residual and temperature_rate fields following the previously described integration scheme.

for (Int s = 1; (s-1)*applied_time_step < total_time; ++s) {
  model.solveStep();
}

An example of explicit dynamic heat propagation is presented in examples/heat_transfer/explicit_heat_transfer.cc. This example consists of a square 2D plate of \(1 \text{m}^2\) having an initial temperature of \(100 \text{K}\) everywhere but a none centered hot point maintained at \(300 \text{K}\). Fig. 24 presents the geometry of this case. The material used is a linear fictitious elastic material with a density of \(8940 \text{kg}/\text{m}^3\), a conductivity of \(401 \text{W}/\text{m}/\text{K}\) and a specific heat capacity of \(385 \text{J}/\text{K}/\text{kg}\). The time step used is \(0.12 \text{s}\).

Fig. 24 Initial temperature field

Fig. 25 Temperature field after 15000 time steps = 30 minutes. The lines represent iso-surfaces.