Homec4science
Diffusion akantu (features/fe_engine_python_interface_extension) Browse

rAKA/doc/manual/manual-heattransfermodel.tex
6e52064950bafeatures/fe_engine_pyth…

Last Change
Raw File

manual-heattransfermodel.tex
View Options

\chapter{Heat Transfer Model}
The heat transfer model is a specific implementation of the \code{Model} interface
dedicated to handle the dynamic heat equation.
\section{Theory}
The strong form of the dynamic heat equation
can be expressed as
\begin{equation}
\rho c_v \dot{T} + \nabla \cdot \vec{\kappa} \nabla T = b
\end{equation}
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
\begin{equation}
\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
\end{equation}
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
\begin{equation}
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
\end{equation}
and the system to solve can be written
\begin{equation}
\mat{C} \cdot \vec{\dot{T}} = \vec{Q}^{\text{ext}} -\mat{K} \cdot \vec{T}~,
\end{equation}
with $\vec{Q}^{\text{ext}}$ the consistent heat generated.
\section{Using the Heat Transfer Model}
A material file name has to be provided during initialization.
Currently, the \code{HeatTransferModel} object uses dynamic analysis
with an explicit time integration scheme. It can simply be created
like this
\begin{cpp}
HeatTransferModel model(mesh, spatial_dimension);
\end{cpp}
while an existing mesh has been used (see \ref{sect:common:mesh}).
Then the model object can be initialized with:
\begin{cpp}
model.initFull()
\end{cpp}
This function will load the material
properties, and allocate / initialize the nodes and element \code{Array}s
More precisely, the heat transfer model contains 4 \code{Arrays}:
\begin{description}
\item[temperature] contains the nodal temperature $T$ (zero by default after the
initialization).
\item[temperature\_rate] contains the variations of temperature $\dot{T}$
(zero by default after the
initialization).
\item[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
\textbf{blocked\_dofs} value of a degree of freedom to \code{true}.
The \textbf{temperature} and the \textbf{temperature\_rate} are
computed for all degrees of freedom where the \textbf{blocked\_dofs}
value is set to \code{false}. For the remaining degrees of freedom,
the imposed values (zero by default after initialization) are kept.
\item[external\_heat\_rate] contains the external heat generations. $\vec{Q^{ext}}$
on the nodes.
\item[internal\_heat\_rate] contains the internal heat generations. $\vec{Q^{int}} = -\mat{K} \cdot \vec{T}$ on the nodes.
\end{description}
Only a single material can be specified on the domain.
A material text file (\eg material.dat) provides the material properties as follows:
\begin{cpp}
model heat_transfer_model [
capacity = %\emph{XXX}%
density = %\emph{XXX}%
conductivity = [%\emph{XXX}% ... %\emph{XXX}%]
]
\end{cpp}
where the \code{capacity} and \code{density} are scalars, and the \code{conductivity} is specified as a $3\times 3$ tensor.
\subsection{Explicit Dynamic}
The explicit time integration scheme in \akantu uses a lumped capacity
matrix $\mat{C}$ (reducing the computational cost, see Chapter \ref{sect:smm}).
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 \code{capacity} \code{Array} of the model.
\begin{cpp}
model.assembleCapacityLumped();
\end{cpp}
\index{HeatTransferModel!assembleCapacityLumped}
\note{Currently, only the explicit time integration with lumped capacity matrix
is implemented within \akantu.} The explicit integration scheme is \index{Forward Euler} \emph{Forward Euler}
\cite{curnier92a}.
\begin{itemize}
\item Predictor: $\vec{T}_{n+1} = \vec{T}_{n} + \Delta t \dot{\vec{T}}_{n}$
\item Update residual: $\vec{R}_{n+1} = \left( \vec{Q^{ext}_{n+1}} - \vec{K}\vec{T}_{n+1} \right)$
\item Corrector : $\dot{\vec{T}}_{n+1} = \mat{C}^{-1} \vec{R}_{n+1}$
\end{itemize}
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:
\begin{cpp}
time_step = model.getStableTimeStep();
\end{cpp}
\index{HeatTransferModel!StableTimeStep}
The stable time step is defined as:
\begin{equation}\label{eqn:htm:explicit:stabletime}
\Delta t_{\st{crit}} = 2 \Delta x^2 \frac{\rho c_v}{\mid\mid \mat{\kappa} \mid\mid^\infty}
\end{equation}
where $\Delta x$ is the characteristic length (\eg the inradius 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.
\begin{cpp}
const Real safety_time_factor = 0.1;
Real applied_time_step = time_step * safety_time_factor;
model.setTimeStep(applied_time_step);
\end{cpp}
The following loop allows, for each time step, to update the \code{temperature}, \code{residual} and
\code{temperature\_rate} fields following the previously described integration scheme.
\begin{cpp}
for (UInt s = 1; (s-1)*applied_time_step < total_time; ++s) {
model.solveStep();
}
\end{cpp}
\index{HeatTransferModel!solveStep()}
An example of explicit dynamic heat propagation is presented in \\
\shellcode{\examplesdir/heat\_transfer/explicit\_heat\_transfer.cc}. \\
This example consists of a square 2D plate of \SI{1}{\metre^2}
having an initial temperature of \SI{100}{\kelvin} everywhere but a
none centered hot point maintained at
\SI{300}{\kelvin}. Figure~\ref{fig:htm:explicit:dynamic} presents
the geometry of this case. The material used is a linear fictitious
elastic material with a density of
\SI{8940}{\kilo\gram\per\metre^3}, a conductivity of
\SI{401}{\watt\per\metre\per\kelvin} and a specific heat capacity of
\SI{385}{\joule\per\kelvin\per\kilogram}. The time step used is
\SI{0.12}{\second}.
\begin{figure}[!htb]
\centering
\includegraphics[width=.4\textwidth]{figures/hot-point-1}
\hfill
\includegraphics[width=.4\textwidth]{figures/hot-point-2}
\caption{Initial temperature field (left) and after 15000 time steps = 30 minutes (right). The lines represent iso-surfaces.}
\label{fig:htm:explicit:dynamic}
\end{figure}
%%% Local Variables:
%%% mode: latex
%%% TeX-master: "manual"
%%% End:
c4science · Help