Tag: laplacianFoam

fvOptions SemiImplicitSource

Last Updated: May 5, 2019

When we solve a partial differential equation with a source term, we must pay attention to how to treat it numerically for better convergence. The SemiImplicitSource fvOption is available in OpenFOAM so that we can handle a linearized source term in numerically stable way.

Keywords

  • Source term linearization
  • SemiImplicitSource in OpenFOAM
Source Term Linearization

If the source term under consideration is a non-linear function of a conserved variable, the linearization of it is a fundamental approach to discretizing it in the finite volume method. This topic is precisely covered in the famous book “Numerical Heat Transfer and Fluid Flow” by Suhas V. Patankar:

In Section 4.2.5, the concept of the linearization of the source term was introduced. One of the basic rules (Rule 3) required that when the source term is linearized as

\begin{equation}
S = S_C + S_P \phi_P, \tag{7.6} \label{eq:sourceTerm}
\end{equation}

the quantity \(S_P\) must not be positive. Now, we return to the topic of source-term linearization to emphasize that often source terms are the cause of divergence of iterations and that proper linearization of the source term frequently holds the key to the attainment of a converged solution.

Rule 3: Negative-slope linearization of the source term If we consider the coefficient definitions in Eqs. (3.18), it appears that, even if the neighbor coefficients are positive, the center-point coefficient \(a_P\) can become negative via the \(S_P\) term. Of course, the danger can be completely avoided by requiring that \(S_P\) will not be positive. Thus, we formulate Rule 3 as follows:

When the source term is linearized as \(\bar{S} = S_C + S_P T_P\), the coefficient \(S_P\) must always be less than or equal to zero.

Linearized Source Specification for a Scalar Field in OpenFOAM

In OpenFOAM, the linearized source terms can be specified using SemiImplicitSource fvOption. Its settings are described in the constant/fvOptions file.

  • selectionMode: [Required] Domain where the source is applied (all/cellSet/cellZone/points)
  • injectionRateSuSp: [Required] Conserved variable name and two coefficients of linearized source term

  • volumemode: [Required] Choice of how to specify the two coefficients

– absolute: values are given as [variable]

– specific: values are given as [variable]/\({\rm m^3}\)

Example 1. Heat Conduction with Heat Source

Let’s consider a simple 1D heat conduction problem in a solid rod with a thermal energy generation per unit volume and time \(\dot{q}_v {\rm [W/m^3]}\). In this case, the steady heat conduction equation in 1D is given by

\begin{align}
\alpha \frac{d^2 T}{d x^2} + \frac{\dot{q}_v}{\rho c} = 0. \tag{1} \label{eq:conductionEqn}
\end{align}

Here, we assume that the thermal diffusivity \(\alpha {\rm [m^2/s]}\) is constant.

We can obtain the general solution for this equation \eqref{eq:conductionEqn} by integrating it twice.

\begin{equation}
T(x) = -\frac{S_C}{2\alpha}x^2 + C_1 x + C_2, \tag{2} \label{eq:generalSolution}
\end{equation}

where \(C_1\) and \(C_2\) are integration constants and we put \(S_C {\rm [K/s]} := \dot{q}_v/\rho c\). If we assume that the temperatures at both ends of the rod are maintained at constant temperatures \(T(0) = T_1, T(L) = T_2\), we can determine the values of two integration constants as

\begin{align}
C_2 &= T_1, \tag{3} \\
C_1 &= \frac{T_2 – T_1}{L} + \frac{S_C L}{2\alpha}. \tag{4}
\end{align}

Then, the temperature distribution in the rod is expressed by the following equation:

\begin{equation}
T(x) = -\frac{S_C}{2\alpha}x^2 + \left( \frac{T_2 – T_1}{L} + \frac{S_C L}{2\alpha} \right)x + T_1. \tag{5} \label{eq:solution}
\end{equation}

t_heatsource_case1

t_heatsource

If we confine the heat source region to one-quarter of the whole region as shown in the following picture, the quadratic temperature distribution is limited to the region and the linear profile is obtained in the source free region.

t_heatsource_case2_ps2

t_heatsource_case2

t_heatsource_case2b

Example 2. Passive Scalar Transport

Introduction to laplacianFoam and simple validation calculation


In this blog post, I will try to give a description of the governing equation of the laplacianFoam in OpenFOAM that solves a simple Laplace equation, e.g. for thermal diffusion in a solid.

Governing Equation

The heat conduction equation is given by the following equation:

\begin{align}
\frac{\partial T}{\partial t} = \frac{1}{\rho c_p} \nabla \cdot \left(k \nabla T\right) + \frac{q}{\rho c_p}, \tag{1} \label{eq:conductionEqn}
\end{align}

where \(T\;{\rm [ K]}\) is the absolute temperature field, \(\rho\;{\rm [kg/m^3]}\) is the density field, \(q\;{\rm [W/m^3]}\) is the rate of energy generation per unit volume, \(k\;{\rm [W/(m\cdot K)]}\) is the thermal conductivity and \(c_p\;{\rm [J/(kg\cdot K)]}\) is the specific heat at constant pressure.

If the heat capacity \(\rho c_p\) is spatially uniform, the Eq. \eqref{eq:conductionEqn} can be transformed into the following form irrespective of whether the thermal conductivity \(k\) is spatially uniform or not:

\begin{align}
\frac{\partial T}{\partial t} = \nabla \cdot \left(\alpha \nabla T\right) + \frac{q}{\rho c_p}, \tag{2} \label{eq:conductionEqn2}
\end{align}

where \(\alpha = k/\rho c_p\;{\rm [m^2/s]}\) is the thermal diffusivity.

The laplacianFoam (in OpenFOAM-4.x and earlier versions) doesn’t consider the heat generation and the implemented equation is

\begin{align}
\frac{\partial T}{\partial t} = \nabla \cdot \left(\alpha \nabla T\right), \tag{3} \label{eq:laplacianFoam}
\end{align}

but the solver in the latest development version supports the fvOptions so that we can solve \eqref{eq:conductionEqn2} and specify a volumetric heat source.

We can find the Eq. \eqref{eq:conductionEqn2} solved in laplacianFoam.C.

The variable DT represents the thermal diffusivity \(\alpha\) and it is specified in the constant/transportProperties file.

Comparison with Analytical Solution

We consider the steady state problem of source-free heat conduction in a concentric cylinder whose inner and outer walls are maintained at constant temperature of 400 K and 300 K respectively. This problem can be analytically solved and the following relation holds for the temperature distribution in the radial direction \(T(r)\):

\begin{align}
\frac{T_1 – T(r)}{T_1 – T_2} = \frac{\ln{(r/r_1)}}{\ln{(r_2/r_1)}}, \tag{4} \label{eq:cylinderT}
\end{align}

where the subscripts 1 and 2 indicate the values at the inner and outer walls respectively shown in Figure 1.

conduction_pb1
Fig. 1 Problem setting

The temperature distribution obtained using laplacianFoam is shown in Figure 2 and the comparison between the numerical and analytical solutions is shown in Figure 3. The agreement with the analytical solution is good.

Fig. 2
Fig. 2 Temperature distribution
conduction_pb1_comparison
Fig. 3 Comparison of numerical and analytical solutions
Introduction of Source Term

In the next post, I’ll deal with a simple example case to introduce how to specify a heat generation source term in laplacianFoam using the fvOptions functionality (scalarSemiImplicitSource).

t_heatsource