Warning

Before we start the exercise 1, please download (test_laplacianFoam) the modified Laplacian solver first. Then put it in the shared folder and compile (wmake) it in the docker container. This modified solver writes out detailed information on matrix coefficients into the log file, so that it is possible to check the derived values against the computed values.

Exercise 1

Theory

Goals

• explore FVM theory

• link FVM theory and OpenFOAM solver implementation

• step by step exploring details of the solver and solving process

Step 1, Geometric and physical modeling

(27)\[\rho c_p\frac{\partial T}{\partial t} - \nabla \cdot k \nabla T = 0\]

(28)\[\frac{\partial T}{\partial t} - \nabla \cdot D \nabla T =0\]

\(k\) is the thermal conductivity and \(D\) is the thermal diffusivity \(\frac{k}{\rho c_p}\) . We here assume that \(D\) is a constant value. see also Theory.

Below is a simple 2D setup for a heat diffusion test case. We will use it to learn about the details of how openFOAM’s FV method works.

Fig. 34 Model geometry, boundary conditions and initial condition.

What we need to solve ?

Temperature field \(T\) of the model region at time \(t\), therefore \(T\) is the only unknown variable of our problem.

Step 2, Domain discretization

Mesh structure and topology

Create a 2D regular box mesh using blockMesh utility (Regular box case). Therefore, the mesh topology shown below is the same as OpenFOAM polyMesh.

Regular mesh

Fig. 35 2D regular mesh structure and topology information of OpenFOAM polyMesh (See 1. Read and plot mesh information in the notebook).

Unstructured mesh

Fig. 36 2D unstructured mesh topology information (Regular box case).

Tip

All internal faces have and only have two attributes, owner and neighbour, denote the cell indices who share the face. While the boundary faces only have owner attribute. The normal vector of face always points from the owner cell to the neighbour cell. The numbering of the vertices that make up the face is again done using the right-hand rule and the face normal always points from the cell with the lower index to the cell with the larger index.

You can find details on the mesh description in the OpenFoam manual. Have a look at it and check the structure of the files in constant\polymesh.

Step 3, Spatial discretization: The diffusion term

The equation discretization step is performed over each element of the computational domain to yield an algebraic relation that connects the value of a variable in an element to the values of the variable in the neighboring elements [Moukalled et al., 2016].

Goal

Transform partial differential equation (28) into a set of algebraic equations: \(\mathbf{A}[T] = \mathbf{b}\).

1. Integrate equation (28) over element \(C\) (the green cell in Fig. 35) that enables recovering its integral balance form as,

(29)\[\iiint\limits_{V_C} \frac{\partial T}{\partial t} dV - \iiint\limits_{V_C} \nabla \cdot D\nabla T dV = 0\]

Note that the transient term \(\iiint\limits_{V_C} \frac{\partial T}{\partial t} dV\) will be processed in the later section. Let’s set it to zero for the moment and assume steady-state.

2. Transform the volume integral of the heat flux into a surface integral by applying the divergence theorem,

(30)\[- \iiint\limits_{V_C} \nabla \cdot D\nabla T dV = -\oint\limits_{\partial V_C} (D\nabla T)\cdot d\vec{S} = 0\]

The equation (30) is actually a heat balance over cell \(C\). It is basically the integral form of the original partial differential equation and involves no approximation.

The integrant is the diffusive heat flux,

(31)\[\vec{J}^{T,D} \equiv -D\nabla T\]

3. Transform the surface integral in equation (30) as a summation over the control volume faces (still no approximation),

(32)\[\sum\limits_{f\sim faces(V_C)} \iint\limits_{f}\vec{J}^{T,D}_f \cdot d\vec{S}_f =0\]

Here \(f\) denotes the boundary face of cell \(V_C\).

Now, the key problem is how to calculate flux integral on a face.

4. Evaluating flux integral on the faces

It’s worth noting that there are only two kinds of cells in the discretized domain, one is internal cell, the other is boundary cell. All faces of a internal cell are internal faces. The remain cells are boundary cells which contain at least one boundary face.

Internal cell

For a specific internal cell \(C\), e.g. cell 12, shown in Fig. 35, the equation (32) can be expanded as,

(33)\[\color{red}{\iint_{f_{23}}\vec{J}^{T,D}_{f_{23}} \cdot \vec{S}_{f_{23}}} + \iint_{f_{24}}\vec{J}^{T,D}_{f_{24}} \cdot \vec{S}_{f_{24}} + \iint_{f_{21}}\vec{J}^{T,D}_{f_{21}} \cdot \vec{S}_{f_{21}} + \iint_{f_{5}}\vec{J}^{T,D}_{f_{5}} \cdot \vec{S}_{f_{5}} = 0\]

Note that the surface normals \(\vec{S}\) always point out of the cell; for example, \(\vec{S}_{f_{23}}\) would have a positive and \(\vec{S}_{f_{21}}\) a negative sign.

4.1 Considering face \(f_{23}\), calculate the first term on the right hand side of the equation (33) (we have to introduce the first approximation at this step). Using a Gaussian quadrature the integral at the face \(f_{23}\), for example, becomes,

(34)\[\color{red}{\iint_{f_{23}}\vec{J}^{T,D}_{f_{23}} \cdot \vec{S}_{f_{23}}} = \iint_{f_{23}} (\vec{J}^{T,D}_{f_{23}} \cdot \vec{n}_{f_{23}}) dS_{f_{23}} \approx \color{orange}{\sum\limits_{ip\sim ip(f_{23})} (\vec{J}^{T,D} \cdot \vec{n}_{ip})\omega_{ip} S_{f_{23}}}\]

where \(S_f\) is the area of face \(f_{23}\), \(ip\) refers to a integration point and \(ip(f_{23})\) the number of integration points along surface \(f_{23}\), \(\omega_{ip}\) is the integral weights.

Fig. 37 Surface integration of fluxes using (a) one integration point, (b) two integration points, and (c) three integration points [Moukalled et al., 2016]. The total flux \(FluxT_f\) is computed as linear combinations of the two cell values plus a source term.

Important

Only the one integration point scheme, i.e. Fig. 37 (a), is implemented in OpenFOAM (at least before version 8). Therefore the first keyword of Laplacian scheme in system/fvScheme dictionary file of a case has only one option, it is Gauss.

4.2 Choose integral scheme or integral points

To simply explain the calculation process/logic and to also be consistent with OpenFOAM, we here adopt an one integration point scheme with weight \(\omega = 1\), thus equation (34) becomes,

(35)\[\color{orange}{\sum\limits_{ip\sim ip(f_{23})} (\vec{J}^{T,D} \cdot \vec{n}_{ip})\omega S_f } = -\left(D \frac{\partial T}{\partial x} \vec{i} + D\frac{\partial T}{\partial y} \vec{j} \right)_{f_{23}} \cdot \Delta y_{f_{23}} \vec{i} = -\color{blue}{\left( D\frac{\partial T}{\partial x} \right)_{f_{23}} \Delta y_{f_{23}}}\]

where

• \(S_{f_{23}} = \Delta y_{f_{23}}\) is the area of face \(f_{23}\) by assuming \(\Delta z = 1\).

• \(\vec{S}_{f_{23}} = S_{f_{23}} \vec{n}_{f_{23}}\) is the surface vector of face \(f_{23}\).

• \(\vec{n}_{f_{23}} = \vec{i}\) is the normal vector of the face \(f_{23}\) directed out of the cell \(C\) (see Fig. 35).

4.3 Calculate gradient of \(T\) at face centroid, here introduce the second approximation, e.g. assuming linear variation of T and then the gradient term in the equation (35) can be approximated as,

\[-\color{blue}{\left( D\frac{\partial T}{\partial x} \right)_{f_{23}} \Delta y_{f_{23}}} = -D\frac{T_{13} - T_{C}}{x_{13} - x_C}\Delta y_{f_{23}} = -\frac{D \Delta y_{f_{23}}} {\delta x_{f_{23}}} (T_{13} - T_{C})\]

where \(\delta x_{f_{23}}\) represents the distance between center of cells who share face \(f_{23}\), they are cell 12 and cell 13.

Now let’s look at the western face \(f_{21}\) . Using a single integration point, we can again write:

(36)\[\color{orange}{\sum\limits_{ip\sim ip(f_{21})} (\vec{J}^{T,D} \cdot \vec{n}_{ip})\omega S_f } = -\left(D \frac{\partial T}{\partial x} \vec{i} + D\frac{\partial T}{\partial y} \vec{j} \right)_{f_{21}} \cdot \Delta y_{f_{21}} \cdot-\vec{i} = \color{blue}{\left( D\frac{\partial T}{\partial x} \right)_{f_{21}} \Delta y_{f_{21}}}\]

Notice how the surface normal now has the opposite sign!

We continue to spell out the flux across face \(f_{21}\) as:

\[\color{blue}{\left( D\frac{\partial T}{\partial x} \right)_{f_{21}} \Delta y_{f_{21}}} = D\frac{T_{C} - T_{11}}{x_{11} - x_C}\Delta y_{f_{21}} = -\frac{D \Delta y_{f_{21}}} {\delta x_{f_{21}}} (T_{11} - T_{C})\]

Do the same thing for the remain faces (\(f_5, f_{24}\)).

4.3 Construct coefficients of a specific cell \(C_{12}\)

Now we get all the coefficients, thus the equation (33) can be discretized and expressed as a matrix form,

(37)\[a_{C} T_C + a_{F_{23}} T_{F_{23}} + a_{F_{24}} T_{F_{24}} + a_{F_{21}} T_{F_{21}} + a_{F_{5}} T_{F_{5}} = 0\]

with

\[\begin{split}\begin{eqnarray} a_{23} & = & -\frac{D \Delta y_{f_{23}}} {\delta x_{f_{23}}} \\ a_{24} & = & -\frac{D \Delta y_{f_{24}}} {\delta x_{f_{24}}} \\ a_{21} & = & -\frac{D \Delta y_{f_{21}}} {\delta x_{f_{21}}} \\ a_{5} & = & -\frac{D \Delta y_{f_{5}}} {\delta x_{f_{5}}} \\ a_{C} & = & -(a_{23} + a_{24} + a_{21} + a_{5}) \\ \end{eqnarray}\end{split}\]

(38)\[\mathbf{A}_{N\times N}[12,:] \mathbf{T}_{N\times 1} = 0\]

Boundary cell

For a specific boundary cell, e.g. cell 19 shown in Fig. 35, the equation (32) can be rewritten as,

(39)\[\iint_{f_{18}}\vec{J}^{T,D}_{f_{18}} \cdot \vec{S}_{f_{18}} + \iint_{f_{35}}\vec{J}^{T,D}_{f_{35}} \cdot \vec{S}_{f_{35}} +\iint_{f_{37}}\vec{J}^{T,D}_{f_{37}} \cdot \vec{S}_{f_{37}} + \color{red}{\iint_{f_{51}}\vec{J}^{T,D}_{f_{51}} \cdot \vec{S}_{f_{51}}} = 0\]

The first three terms are normal fluxes across internal faces, so there is nothing special and we can just calculate them following steps for internal faces as before. The special thing is the boundary face \(f_{52}\), we can call it \(f_b\) for a general case. Now the problem is how to evaluate flux on the boundary face \(f_b\).

The fluxes on the interior faces are discretized as before, while the boundary flux is discretized with the aim of constructing a linearization with respect to the cell field \(T_C\), e.g. \(T_{C_{19}}\) of cell \(C_{19}\) shown in Fig. 35, thus

(40)\[\color{red}{\iint_{f_{b}}\vec{J}^{T,D}_{f_{b}} \cdot d\vec{S}_{f_{b}}} = a_{F_b} T_C + c_{F_b}\]

All right, now out goal is to determine the coefficients of \(a_{F_b}\) and \(c_{F_b}\) according to the boundary conditions. Generally there are two basic kinds of boundary conditions, they are Dirichlet boundary condition (also called the first type or fixed value boundary condition) and Von Neumann boundary condition (all called the second type or fixed flux boundary condition), respectively.

4.1 Fixed value boundary condition

Fixed value means the value of field \(T\) is given on the boundary face is give, i.e. \(\color{purple}{T_{f_b}}\) is given. There fore we can calculate flux on boundary face, \(f_{51}\) for example,

(41)\[\color{red}{\iint_{f_{51}}\vec{J}^{T,D}_{f_{51}} \cdot \vec{S}_{f_{51}}} = \iint_{f_{51}} (\vec{J}^{T,D}_{f_{51}} \cdot \vec{n}_{f_{51}}) dS_{f_{51}} \approx \color{orange}{\sum\limits_{ip\sim ip(f_{51})} (\vec{J}^{T,D} \cdot \vec{n}_{ip})\omega_{ip} S_{f_{51}}}\]

Here we still use one integration point scheme to explain the calculation process,

(42)\[\color{orange}{\sum\limits_{ip\sim ip(f_{51})} (\vec{J}^{T,D} \cdot \vec{n}_{ip})\omega S_f } = -\left(D \frac{\partial T}{\partial x} \vec{i} + D\frac{\partial T}{\partial y} \vec{j} \right)_{f_{51}} \cdot \Delta y_{f_{51}} \vec{i} = -\color{blue}{\left( D\frac{\partial T}{\partial x} \right)_{f_{51}} \Delta y_{f_{51}}}\]

Now let’s calculate the gradient on the boundary face,

(43)\[-\color{blue}{\left( D\frac{\partial T}{\partial x} \right)_{f_{51}} \Delta y_{f_{51}}} = -D\frac{\color{purple}{T_{f_{51}}} - T_{C}}{x_{f_{51}} - x_C}\Delta y_{f_{51}} = -\frac{D \Delta y_{f_{51}}}{\delta x_{f_{51}}}(\color{purple}{T_{f_{51}}} - T_{C}) = a_{f_{51}}(-T_{C}) + c_{f_{51}}\]

Substituting equation (43) back to equation (42), equation (41) and equation (40), we can get Laplacian discretization coefficients \(a_{f_{51}}\) and \(c_{f_{51}}\),

(44)\[\begin{split}\begin{eqnarray} a_{F_{51}} &=& -D\frac{\Delta y_{f_{51}}}{\delta x_{f_{51}}}\\ c_{F_{51}} &=& a_{F_{51}}\color{purple}{T_{f_{51}}} \end{eqnarray}\end{split}\]

Notice how \(\color{purple}{T_{f_{51}}}\) is known (a boundary condition), so that the term \(-a_{F_{51}}\color{purple}{T_{f_{51}}}\) will end up on the right-hand side and not in the coefficient matrix.

4.2 Fixed flux boundary condition

Fixed flux means normal gradient to the boundary face of field \(T\) is given, i.e. \(\color{purple}{\vec{J}^{T,D}_{f_{b}} \cdot \vec{n}_{f_{b}} = q_{f_b}}\) is given. There fore we can calculate flux on boundary face, \(f_{51}\) for example,

(45)\[\color{red}{\iint_{f_{51}}\vec{J}^{T,D}_{f_{51}} \cdot \vec{S}_{f_{51}}} = \iint_{f_{51}} (\vec{J}^{T,D}_{f_{51}} \cdot \vec{n}_{f_{51}}) dS_{f_{51}} = \iint_{f_{51}} \color{purple}{q_{f_{51}}} dS_{f_{51}} = \color{purple}{q_{f_{51}}} S_{f_{51}} = \color{purple}{q_{f_{51}}} \Delta y_{f_{51}}\]

Substituting equation (45) back to equation (40), we can get Laplacian discretization coefficients \(a_{f_{51}}\) and \(c_{f_{51}}\),

(46)\[\begin{split}\begin{eqnarray} a_{F_{51}} &=& 0 \\ c_{F_{51}} &=& \color{purple}{q_{f_{51}}} \Delta y_{f_{51}} \end{eqnarray}\end{split}\]

4.3 Construct coefficients of the boundary cell \(C_{19}\)

Do the same thing as 4.1 or 4.2 (depends on the type of boundary condition) for all boundary faces of a boundary cell \(C\) and calculate the coefficients of \(a_{F_b}\) and \(c_{F_b}\), then the equation (39) can be expressed as,

(47)\[\begin{split} \begin{eqnarray} b_{C_{19}} & = & a_{F_{18}} (T_{F_{18}} - T_C) + a_{F_{35}} (T_{F_{35}} - T_C) + a_{F_{37}} (T_{F_{37}} - T_C) + a_{F_{51}} (-T_C) +c_{F_{51}} \\ & = & \sum\limits_{f\in[18, 35, 37]} a_{F_f} T_{F_f} + \left(-\sum\limits_{f\in[18, 35, 37, 51]} a_{F_f} \right) T_C + \sum\limits_{f\in[51]} c_{F_f} \\ & = & \sum\limits_{f\sim nb(\text{internal face})} a_{F_f} T_{F_f} + \left(-\sum\limits_{f\sim nb(\text{all face})} a_{F_f} \right) T_C + \sum\limits_{f\sim nb(\text{boundary face})} c_{F_f} \end{eqnarray}\end{split}\]

(48)\[\mathbf{A}_{N\times N}[19,:] \mathbf{T}_{N\times 1} = B_{19}\]

where \(B_{19} = b_{C_{19}} - \sum\limits_{f\sim nb(\text{boundary face})} c_{F_f}\).

5. Construct general matrix form of the equation (32)

(49)\[\begin{split}\begin{eqnarray} -\sum\limits_{f\sim nb(C)} (D\nabla T)_{F_f}\cdot \vec{S}_{F_f} & = & \sum\limits_{f\sim \text{internal faces}(C)} a_{F_f} T_{C} + \left(-\sum\limits_{f\sim nb(C)} a_{F_f} \right)T_{F_f} + \sum\limits_{f\sim \text{boundary faces}(C)} c_{F_{F_f}} \\ & = & a_CT_C - \sum\limits_{f\sim \text{internal faces}(C)} a_{F_f} T_{F_f} + \sum\limits_{f\sim \text{boundary faces}(C)} c_{F_{f}} \end{eqnarray}\end{split}\]

Matrix form of the equation (32) can be written as,

(50)\[\mathbf{A}_{N \times N} \mathbf{T}_{N\times1} = \mathbf{B}_{N\times1}\]

here \(N\) is the number of cell, \(\mathbf{A}\) is a sparse matrix of coefficients, \(\mathbf{T}\) is cell-centered temperature field. The matrix is visualized as Fig. 38. It should be noted that the coefficients matrix of Laplacian term is a symmetric matrix and the boundary conditions only affect diagonal entry of the coefficients matrix. Further, only the fixed value boundary condition affects the matrix and the fixed flux boundary condition affects the RHS(Right hand of side) \(B\).

Regular mesh

Fig. 38 Visualization of \(\mathbf{A}_{N \times N} \mathbf{T}_{N\times1} = \mathbf{b}_{N\times1}\). The coefficients of the selected cell \(C\) (Fig. 35) are marked by green rectangle.

Unstructured mesh

Fig. 39 Visualization of \(\mathbf{A}_{N \times N} \mathbf{T}_{N\times1} = \mathbf{b}_{N\times1}\). The coefficients of the selected cell \(C\) (Fig. 35) are marked by green rectangle.

Step 4, Temporal Discretization: The Transient Term

After spatial discretization, the equation (29) can be expressed as,

(51)\[\frac{\partial T}{\partial t} V_C - L(T_C^t) =0\]

where \(V_C\) is the volume of the discretization cell and \(L(T_C^t)\) is the spatial discretization operator expressed at some reference time \(t\), which can be written as algebraic form (see also equation (49)),

(52)\[L(T_C^t) = a_C T_C^t + \sum\limits_{F\sim NB(C)} a_F T_F^t + \sum\limits_{F\sim NB(C)} c_F\]

where \(a_C\) is the diagonal coefficients of the matrix, \(a_F\) is the off-diagonal coefficients, and the \(c_F\) is the source coefficients as right hand side of matrix of system.

Tip

For specific cell \(C\),

• \(a_F\) is only contributed from internal faces.

• \(a_C\) is the negative summation of \(a_F\), contributed from all faces, but equation of the \(a_F\) for a boundary faces (equation (44) and equation (46)) is a little bit different from internal face.

• \(c_F\) only comes from boundary faces of cell \(C\) if it has boundary face, see also equation (44) and equation (46). \(c_F\) will contribute to RHS (\(\mathbf{B}\)) of the algebraic equation (50).

Let’s do the temporal discretization for the first term of equation (51),

\[\frac{T^{t+\Delta t/2} - T^{t - \Delta t/2}}{\Delta t} V_C + L(T_C^t) = 0\]

To derive the full discretized equation, an interpolation profile expressing the face values at (\(t-\Delta t/2\)) and (\(t+\Delta t\)) in terms of the element values at (\(t\)), (\(t-\Delta t\)), etc., is needed.Independent of the profile used, the flux will be linearized based on old and new values as,

(53)\[b_C = FluxC T_C + FluxC^o T_C^o + FluxV\]

where the superscript \(^{o}\) refers to old values. With the format of the linearization, the coefficient \(FluxC\) will be added to diagonal element and the coefficient \(FluxC^o T_C^o + FluxV\) will be added to the source or RHS (\(B\)).

First order implicit Euler scheme

(54)\[\frac{T^t - T^{t-\Delta t}}{\Delta t} V_C - L(T^{t}) = 0\]

Therefore the coefficients will be

(55)\[FluxC = \frac{V_C}{\Delta t}, ~ FluxC^o = - \frac{V_C}{\Delta t}, ~ FluxV =0\]

First order explicit Euler scheme

(56)\[\frac{T^{t+\Delta t} - T^t}{\Delta t} V_C - L(T^{t}) = 0\]

Note that now the new time is at \(t+\Delta t\) and that the spatial operator (\(L\)) of equation (56) has to be evaluated at time \(t\). Therefore, in order to find the value of \(T\) at time \(t+\Delta t\), we don’t need to solve a set of linear algebraic equations,

(57)\[T^{t+\Delta t} = \frac{\Delta t}{V_C} L(T^{t}) + T^t\]

Step 6: Solution of the discretized equations

The discretization of the differential equation results in a set of discrete algebraic equations, which must be solved to obtain the discrete values of T. The coefficients of these equations may be independent of T (i.e., linear) or dependent on T (i.e. non-linear). The techniques to solve this algebraic system of equations are independent of the discretization method. The solution methods for solving systems of algebraic equations may be broadly classified as direct or iterative.

OpenFOAM implementation

Ok, now let’s do some practical tings, (0) generate mesh; (1) read mesh and do some useful calculation, e.g. cell volume, face area, …, etc; (2) discretize Laplacian term and get coefficients; (3) discretize transient term and get additional coefficients; (4) construct the final coefficient matrix and RHS; (5) solve the system of algebraic equations; (6) write solution to file.

Goal

Deeply look into OpenFOAM and understand how it works!

In order to better understand OpenFOAM’s logic and its work flow, we have to look at the basic structure of the source code of a basic solver, laplacianFoam.

laplacianFoam.C

Listing 20 Source code of laplacianFoam

#include "fvCFD.H"               // Basic head file of OF
#include "simpleControl.H"       // Basic head file of OF
int main(int argc, char *argv[]) // Typical c++ main control function
{
   #include "setRootCaseLists.H" // Do some case/file path-related thing
   #include "createTime.H"       // Create a time object: read controlDict, ...
   #include "createMesh.H"       // (2) Create mesh object: read mesh and do some useful calculation
   simpleControl simple(mesh);   // Time loop control object
   #include "createFields.H"     // (3) Read input data: T and D in the PDE

   Info<< "\nCalculating temperature distribution\n" << endl;
   while (simple.loop(runTime)) // do time loop
   {
      Info<< "Time = " << runTime.timeName() << nl << endl;
      while (simple.correctNonOrthogonal())  // Non-orthogonal correction loop for the unstructured/non-orthogonal mesh
      {
         fvScalarMatrix TEqn                 // (5) construct the final coefficient matrix, include RHS (.source)
            (
               fvm::ddt(T)                   // (4) Discretization of the transient term, return a fvMatrix object
               -
               fvm::laplacian(DT, T)         // (3) Discretization of the Laplacian, return a fvMatrix object
            );
            TEqn.solve();                    // (6) Solve the system of algebraic equations
      }
      #include "write.H"                     // (7) Save solution to file.
      Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
            << nl << endl;
   }
   Info<< "End\n" << endl;
   return 0;
}

Table 1 Comparison between equation and code

Term	Equation	Code	Coefficients(ldu:b)	Reference
Transient	\(\frac{\partial T}{\partial t}\)	fvm::ddt(T)	d(\(FluxC\)), b(\(FluxC^o\))	equation (57)
Laplacian	\(\nabla \cdot D \nabla T\)	fvm::laplacian(DT, T)	d(\(\sum\limits_{F\sim NB(C)}a_F\)), u(\(a_F\) only internal faces)	equation (37), equation (44), equation (46)

createFields.H

Listing 21 createFields.H

volScalarField T //Read input data of T, include ICs and BCs
(
   IOobject
   (
      "T",
      runTime.timeName(),
      mesh,
      IOobject::MUST_READ,
      IOobject::AUTO_WRITE
   ),
   mesh
);

IOdictionary transportProperties // Read dictionary file of transportProperties
(
   IOobject
   (
      "transportProperties",
      runTime.constant(),
      mesh,
      IOobject::MUST_READ_IF_MODIFIED,
      IOobject::NO_WRITE
   )
);
dimensionedScalar DT // Read the diffusivity D (constant) from transportProperties object
(
   transportProperties.lookup("DT")
);

Input data 1

Listing 22 constant/transportProperties

FoamFile
{
   version     2.0;
   format      ascii;
   class       dictionary;
   location    "constant";
   object      transportProperties;
}
DT              DT [0 2 -1 0 0 0 0] 4e-05;

Input data 2

Listing 23 0/T

FoamFile
{
   version     2.0;
   format      ascii;
   class       volScalarField;
   object      T;
}
dimensions      [0 0 0 1 0 0 0];
internalField   uniform 273;
boundaryField
{
   left
   {
      type            fixedValue;
      value           uniform 273;
   }
   right
   {
      type            fixedValue;
      value           uniform 573;
   }
   "(top|bottom)"
   {
      type            zeroGradient;
   }
}

Step0, Generate mesh

Generate mesh

Nothing special, just a meshing process. Here we use the OpenFOAM utility blockMesh to generate a regular mesh (Regular box case) same as Fig. 35.

OpenFOAM script

Listing 24 blockMeshDict of a regular box mesh we shown above

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  5                                     |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
   version     2.0;
   format      ascii;
   class       dictionary;
   object      blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

convertToMeters 1;
xmin -0.05;    //variable definition
xmax 0.05;
ymin -0.015;
ymax 0.015;
Lz 0.01;
vertices    //vertices definition
(
   ($xmin      $ymin   0)  //coordinate of vertex 0
   ($xmax    $ymin   0)  //coordinate of vertex 1
   ($xmax    $ymax   0)  //coordinate of vertex 2
   ($xmin      $ymax   0)  //coordinate of vertex 3
   ($xmin      $ymin   $Lz)//coordinate of vertex 4
   ($xmax    $ymin   $Lz)//coordinate of vertex 5
   ($xmax    $ymax   $Lz)//coordinate of vertex 6
   ($xmin      $ymax   $Lz)//coordinate of vertex 7
);
blocks
(
   hex (0 1 2 3 4 5 6 7) (10 3 1) simpleGrading (1 1 1)
);
boundary
(
   left    //patch name
   {
      type patch ;
      faces   //face list
      (
            (0 4 7 3)
      );
   }
   right
   {
      type patch;
      faces
      (
            (2 6 5 1)
      );
   }
   top
   {
      type patch;
      faces
      (
            (3 7 6 2)
      );
   }
   bottom
   {
      type patch;
      faces
      (
            (1 5 4 0)
      );
   }
   frontAndBack    //patch name
   {
      type empty;
      faces //face list
      (
            (0 3 2 1)   //back face
            (4 5 6 7)   //front face
      );
   }
);
// ************************************************************************* //

Ok, let’s exploring the main steps of the laplacianFoam. Please download (test_laplacianFoam) and do the following practice/debug steps in test_laplacianFoam.C.

Step 1, Read mesh and input field

Read data

The basic properties of cells and faces, e.g. area, volume, face normal vector, will be evaluated after mesh reading, all these processes are happened in the mesh object. It means that after calling creatFields.H all these properties are evaluated and stored in the mesh object. Of course the temperature field object T (volScalarField) with BCs and ICs, and thermal diffusivity DT are also initialized from input data after calling creatFields.H. It should be noted that the part of Laplacian discretization coefficients are also calculated in this step. If the mesh is not changed during simulation time, the mesh related coefficients just need to be calculated one time.

Table 2 Mesh- and Field-related coefficients

Object	Equation	Code	Faces	Reference
mesh	\(1/\delta_{C\leftrightarrow F}\)	mesh.deltaCoeffs()	All faces	equation (37), equation (44), equation (46)
T	Dependents on BCs type	gradientInternalCoeffs (diagonal), gradientBoundaryCoeffs (source)	Boundary patch/faces	equation (44), equation (46)

Access mesh and field properties

Listing 25 Access mesh properties and field boundary properties

// (1). read data
#include "createMesh.H"
#include "createFields.H"
// 1.1 access internal faces
Info<<"\n\nAccess internal mesh"<<endl;
surfaceVectorField Cf = mesh.Cf();
surfaceVectorField Sf = mesh.Sf();
surfaceScalarField S = mesh.magSf();
forAll(Cf, iface)
{
   Info<<iface<<": face center "<<Cf[iface]<<endl;
   Info<<iface<<": face area vector "<<Sf[iface]<<endl;
   Info<<iface<<": face area "<<S[iface]<<endl;
   Info<<iface<<": face delta coeff "<<mesh.deltaCoeffs()[iface]<<endl;
   Info<<iface<<": coeff(D*magSf*deltacoeff) "<<mesh.deltaCoeffs()[iface]*DT*S[iface]<<endl;
}
// 1.2. access boundary mesh
Info<<"\n\nAccess boundary mesh"<<endl;
const fvBoundaryMesh& boundaryMesh = mesh.boundary();
forAll(boundaryMesh, patchI)
{
   const fvPatch& patch = boundaryMesh[patchI];
   forAll(patch, faceI)
   {
      Info<<"Patch "<<patch.name()<<" face "<<faceI<<": face center "<<patch.Cf()[faceI]<<endl;
      Info<<"Patch "<<patch.name()<<" face "<<faceI<<": face area vector "<<patch.Sf()[faceI]<<endl;
      Info<<"Patch "<<patch.name()<<" face "<<faceI<<": face area "<<patch.magSf()[faceI]<<endl;
      Info<<"Patch "<<patch.name()<<" face "<<faceI<<": face delta coeff "<<patch.deltaCoeffs()[faceI]<<endl;
      Info<<"Patch "<<patch.name()<<" face "<<faceI<<": owner cell "<<patch.faceCells()[faceI]<<endl;
   }
}
// 1.3. access boundary field, boundary field coefficients,
forAll(T.boundaryField(), patchI)
{
   Info<<"Boundary patch: "<<mesh.boundary()[patchI].name()<<endl;
   Info<<"Is coupled ? "<<mesh.boundary()[patchI].coupled()<<endl;
   Info<<"gradientInternalCoeffs of field T "<<endl;
   Info<<T.boundaryField()[patchI].gradientInternalCoeffs()<<endl; //Diagonal coeff [A]
   Info<<"gradientBoundaryCoeffs of field T "<<endl;
   Info<<T.boundaryField()[patchI].gradientBoundaryCoeffs()<<"\n"<<endl; //source coeff, [B]
}

Internal mesh

Listing 26 Internal mesh properties of the regular mesh shown in Fig. 35.

46: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07Access internal mesh
0: face center (-0.04 -0.01 0.005)
0: face area vector (0.0001 0 0)
0: face area 0.0001
0: face delta coeff 100
0: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
1: face center (-0.045 -0.005 0.005)
1: face area vector (0 0.0001 0)
1: face area 0.0001
1: face delta coeff 100
1: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
2: face center (-0.03 -0.01 0.005)
2: face area vector (0.0001 0 0)
2: face area 0.0001
2: face delta coeff 100
2: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
3: face center (-0.035 -0.005 0.005)
3: face area vector (0 0.0001 0)
3: face area 0.0001
3: face delta coeff 100
3: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
4: face center (-0.02 -0.01 0.005)
4: face area vector (0.0001 0 0)
4: face area 0.0001
4: face delta coeff 100
4: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
5: face center (-0.025 -0.005 0.005)
5: face area vector (0 0.0001 0)
5: face area 0.0001
5: face delta coeff 100
5: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
6: face center (-0.01 -0.01 0.005)
6: face area vector (0.0001 0 0)
6: face area 0.0001
6: face delta coeff 100
6: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
7: face center (-0.015 -0.005 0.005)
7: face area vector (0 0.0001 0)
7: face area 0.0001
7: face delta coeff 100
7: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
8: face center (0 -0.01 0.005)
8: face area vector (0.0001 0 0)
8: face area 0.0001
8: face delta coeff 100
8: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
9: face center (-0.005 -0.005 0.005)
9: face area vector (0 0.0001 0)
9: face area 0.0001
9: face delta coeff 100
9: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
10: face center (0.01 -0.01 0.005)
10: face area vector (0.0001 0 0)
10: face area 0.0001
10: face delta coeff 100
10: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
11: face center (0.005 -0.005 0.005)
11: face area vector (0 0.0001 0)
11: face area 0.0001
11: face delta coeff 100
11: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
12: face center (0.02 -0.01 0.005)
12: face area vector (0.0001 0 0)
12: face area 0.0001
12: face delta coeff 100
12: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
13: face center (0.015 -0.005 0.005)
13: face area vector (0 0.0001 0)
13: face area 0.0001
13: face delta coeff 100
13: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
14: face center (0.03 -0.01 0.005)
14: face area vector (0.0001 0 0)
14: face area 0.0001
14: face delta coeff 100
14: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
15: face center (0.025 -0.005 0.005)
15: face area vector (0 0.0001 0)
15: face area 0.0001
15: face delta coeff 100
15: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
16: face center (0.04 -0.01 0.005)
16: face area vector (0.0001 0 0)
16: face area 0.0001
16: face delta coeff 100
16: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
17: face center (0.035 -0.005 0.005)
17: face area vector (0 0.0001 0)
17: face area 0.0001
17: face delta coeff 100
17: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
18: face center (0.045 -0.005 0.005)
18: face area vector (0 0.0001 0)
18: face area 0.0001
18: face delta coeff 100
18: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
19: face center (-0.04 0 0.005)
19: face area vector (0.0001 0 0)
19: face area 0.0001
19: face delta coeff 100
19: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
20: face center (-0.045 0.005 0.005)
20: face area vector (0 0.0001 0)
20: face area 0.0001
20: face delta coeff 100
20: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
21: face center (-0.03 0 0.005)
21: face area vector (0.0001 0 0)
21: face area 0.0001
21: face delta coeff 100
21: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
22: face center (-0.035 0.005 0.005)
22: face area vector (0 0.0001 0)
22: face area 0.0001
22: face delta coeff 100
22: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
23: face center (-0.02 0 0.005)
23: face area vector (0.0001 0 0)
23: face area 0.0001
23: face delta coeff 100
23: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
24: face center (-0.025 0.005 0.005)
24: face area vector (0 0.0001 0)
24: face area 0.0001
24: face delta coeff 100
24: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
25: face center (-0.01 0 0.005)
25: face area vector (0.0001 0 0)
25: face area 0.0001
25: face delta coeff 100
25: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
26: face center (-0.015 0.005 0.005)
26: face area vector (0 0.0001 0)
26: face area 0.0001
26: face delta coeff 100
26: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
27: face center (0 0 0.005)
27: face area vector (0.0001 0 0)
27: face area 0.0001
27: face delta coeff 100
27: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
28: face center (-0.005 0.005 0.005)
28: face area vector (0 0.0001 0)
28: face area 0.0001
28: face delta coeff 100
28: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
29: face center (0.01 0 0.005)
29: face area vector (0.0001 0 0)
29: face area 0.0001
29: face delta coeff 100
29: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
30: face center (0.005 0.005 0.005)
30: face area vector (0 0.0001 0)
30: face area 0.0001
30: face delta coeff 100
30: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
31: face center (0.02 0 0.005)
31: face area vector (0.0001 0 0)
31: face area 0.0001
31: face delta coeff 100
31: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
32: face center (0.015 0.005 0.005)
32: face area vector (0 0.0001 0)
32: face area 0.0001
32: face delta coeff 100
32: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
33: face center (0.03 0 0.005)
33: face area vector (0.0001 0 0)
33: face area 0.0001
33: face delta coeff 100
33: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
34: face center (0.025 0.005 0.005)
34: face area vector (0 0.0001 0)
34: face area 0.0001
34: face delta coeff 100
34: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
35: face center (0.04 0 0.005)
35: face area vector (0.0001 0 0)
35: face area 0.0001
35: face delta coeff 100
35: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
36: face center (0.035 0.005 0.005)
36: face area vector (0 0.0001 0)
36: face area 0.0001
36: face delta coeff 100
36: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
37: face center (0.045 0.005 0.005)
37: face area vector (0 0.0001 0)
37: face area 0.0001
37: face delta coeff 100
37: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
38: face center (-0.04 0.01 0.005)
38: face area vector (0.0001 0 0)
38: face area 0.0001
38: face delta coeff 100
38: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
39: face center (-0.03 0.01 0.005)
39: face area vector (0.0001 0 0)
39: face area 0.0001
39: face delta coeff 100
39: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
40: face center (-0.02 0.01 0.005)
40: face area vector (0.0001 0 0)
40: face area 0.0001
40: face delta coeff 100
40: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
41: face center (-0.01 0.01 0.005)
41: face area vector (0.0001 0 0)
41: face area 0.0001
41: face delta coeff 100
41: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
42: face center (0 0.01 0.005)
42: face area vector (0.0001 0 0)
42: face area 0.0001
42: face delta coeff 100
42: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
43: face center (0.01 0.01 0.005)
43: face area vector (0.0001 0 0)
43: face area 0.0001
43: face delta coeff 100
43: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
44: face center (0.02 0.01 0.005)
44: face area vector (0.0001 0 0)
44: face area 0.0001
44: face delta coeff 100
44: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
45: face center (0.03 0.01 0.005)
45: face area vector (0.0001 0 0)
45: face area 0.0001
45: face delta coeff 100
45: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
46: face center (0.04 0.01 0.005)
46: face area vector (0.0001 0 0)
46: face area 0.0001
46: face delta coeff 100
46: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07

Boundary mesh

Listing 27 Boundary mesh properties of the regular mesh shown in Fig. 35.

Patch bottom face 9: owner cell 9Access boundary mesh
Patch left face 0: face center (-0.05 -0.01 0.005)
Patch left face 0: face area vector (-0.0001 0 0)
Patch left face 0: face area 0.0001
Patch left face 0: face delta coeff 200
Patch left face 0: owner cell 0
Patch left face 1: face center (-0.05 0 0.005)
Patch left face 1: face area vector (-0.0001 0 0)
Patch left face 1: face area 0.0001
Patch left face 1: face delta coeff 200
Patch left face 1: owner cell 10
Patch left face 2: face center (-0.05 0.01 0.005)
Patch left face 2: face area vector (-0.0001 0 0)
Patch left face 2: face area 0.0001
Patch left face 2: face delta coeff 200
Patch left face 2: owner cell 20
Patch right face 0: face center (0.05 -0.01 0.005)
Patch right face 0: face area vector (0.0001 0 0)
Patch right face 0: face area 0.0001
Patch right face 0: face delta coeff 200
Patch right face 0: owner cell 9
Patch right face 1: face center (0.05 0 0.005)
Patch right face 1: face area vector (0.0001 0 0)
Patch right face 1: face area 0.0001
Patch right face 1: face delta coeff 200
Patch right face 1: owner cell 19
Patch right face 2: face center (0.05 0.01 0.005)
Patch right face 2: face area vector (0.0001 0 0)
Patch right face 2: face area 0.0001
Patch right face 2: face delta coeff 200
Patch right face 2: owner cell 29
Patch top face 0: face center (-0.045 0.015 0.005)
Patch top face 0: face area vector (0 0.0001 0)
Patch top face 0: face area 0.0001
Patch top face 0: face delta coeff 200
Patch top face 0: owner cell 20
Patch top face 1: face center (-0.035 0.015 0.005)
Patch top face 1: face area vector (0 0.0001 0)
Patch top face 1: face area 0.0001
Patch top face 1: face delta coeff 200
Patch top face 1: owner cell 21
Patch top face 2: face center (-0.025 0.015 0.005)
Patch top face 2: face area vector (0 0.0001 0)
Patch top face 2: face area 0.0001
Patch top face 2: face delta coeff 200
Patch top face 2: owner cell 22
Patch top face 3: face center (-0.015 0.015 0.005)
Patch top face 3: face area vector (0 0.0001 0)
Patch top face 3: face area 0.0001
Patch top face 3: face delta coeff 200
Patch top face 3: owner cell 23
Patch top face 4: face center (-0.005 0.015 0.005)
Patch top face 4: face area vector (0 0.0001 0)
Patch top face 4: face area 0.0001
Patch top face 4: face delta coeff 200
Patch top face 4: owner cell 24
Patch top face 5: face center (0.005 0.015 0.005)
Patch top face 5: face area vector (0 0.0001 0)
Patch top face 5: face area 0.0001
Patch top face 5: face delta coeff 200
Patch top face 5: owner cell 25
Patch top face 6: face center (0.015 0.015 0.005)
Patch top face 6: face area vector (0 0.0001 0)
Patch top face 6: face area 0.0001
Patch top face 6: face delta coeff 200
Patch top face 6: owner cell 26
Patch top face 7: face center (0.025 0.015 0.005)
Patch top face 7: face area vector (0 0.0001 0)
Patch top face 7: face area 0.0001
Patch top face 7: face delta coeff 200
Patch top face 7: owner cell 27
Patch top face 8: face center (0.035 0.015 0.005)
Patch top face 8: face area vector (0 0.0001 0)
Patch top face 8: face area 0.0001
Patch top face 8: face delta coeff 200
Patch top face 8: owner cell 28
Patch top face 9: face center (0.045 0.015 0.005)
Patch top face 9: face area vector (0 0.0001 0)
Patch top face 9: face area 0.0001
Patch top face 9: face delta coeff 200
Patch top face 9: owner cell 29
Patch bottom face 0: face center (-0.045 -0.015 0.005)
Patch bottom face 0: face area vector (0 -0.0001 0)
Patch bottom face 0: face area 0.0001
Patch bottom face 0: face delta coeff 200
Patch bottom face 0: owner cell 0
Patch bottom face 1: face center (-0.035 -0.015 0.005)
Patch bottom face 1: face area vector (0 -0.0001 0)
Patch bottom face 1: face area 0.0001
Patch bottom face 1: face delta coeff 200
Patch bottom face 1: owner cell 1
Patch bottom face 2: face center (-0.025 -0.015 0.005)
Patch bottom face 2: face area vector (0 -0.0001 0)
Patch bottom face 2: face area 0.0001
Patch bottom face 2: face delta coeff 200
Patch bottom face 2: owner cell 2
Patch bottom face 3: face center (-0.015 -0.015 0.005)
Patch bottom face 3: face area vector (0 -0.0001 0)
Patch bottom face 3: face area 0.0001
Patch bottom face 3: face delta coeff 200
Patch bottom face 3: owner cell 3
Patch bottom face 4: face center (-0.005 -0.015 0.005)
Patch bottom face 4: face area vector (0 -0.0001 0)
Patch bottom face 4: face area 0.0001
Patch bottom face 4: face delta coeff 200
Patch bottom face 4: owner cell 4
Patch bottom face 5: face center (0.005 -0.015 0.005)
Patch bottom face 5: face area vector (0 -0.0001 0)
Patch bottom face 5: face area 0.0001
Patch bottom face 5: face delta coeff 200
Patch bottom face 5: owner cell 5
Patch bottom face 6: face center (0.015 -0.015 0.005)
Patch bottom face 6: face area vector (0 -0.0001 0)
Patch bottom face 6: face area 0.0001
Patch bottom face 6: face delta coeff 200
Patch bottom face 6: owner cell 6
Patch bottom face 7: face center (0.025 -0.015 0.005)
Patch bottom face 7: face area vector (0 -0.0001 0)
Patch bottom face 7: face area 0.0001
Patch bottom face 7: face delta coeff 200
Patch bottom face 7: owner cell 7
Patch bottom face 8: face center (0.035 -0.015 0.005)
Patch bottom face 8: face area vector (0 -0.0001 0)
Patch bottom face 8: face area 0.0001
Patch bottom face 8: face delta coeff 200
Patch bottom face 8: owner cell 8
Patch bottom face 9: face center (0.045 -0.015 0.005)
Patch bottom face 9: face area vector (0 -0.0001 0)
Patch bottom face 9: face area 0.0001
Patch bottom face 9: face delta coeff 200
Patch bottom face 9: owner cell 9

Boundary T

Listing 28 Boundary properties of field T of the regular mesh shown in Fig. 35.

0()Boundary patch: left
Is coupled ? 0
gradientInternalCoeffs of field T 
3(-200 -200 -200)
gradientBoundaryCoeffs of field T 
3(54600 54600 54600)

Boundary patch: right
Is coupled ? 0
gradientInternalCoeffs of field T 
3(-200 -200 -200)
gradientBoundaryCoeffs of field T 
3(114600 114600 114600)

Boundary patch: top
Is coupled ? 0
gradientInternalCoeffs of field T 
10{0}
gradientBoundaryCoeffs of field T 
10{0}

Boundary patch: bottom
Is coupled ? 0
gradientInternalCoeffs of field T 
10{0}
gradientBoundaryCoeffs of field T 
10{0}

Boundary patch: frontAndBack
Is coupled ? 0
gradientInternalCoeffs of field T 
0()
gradientBoundaryCoeffs of field T 
0()

Internal cell

Fig. 40 Information of internal cell (\(C_{12}\))

Boundary cell

Fig. 41 Information of boundary cell (\(C_{19}\))

Step 2, discretize Laplacian term

discretize Laplacian term

Because discretization coefficients matrix of Laplacian term is a symmetric matrix, so fvm::Laplacian(DT, T) will return a fvMatrix object only has diagonal and upper. What fvm::Laplacian actually did is evaluate (1) \(a_F\) (see equation (37)) for each internal faces, (2) \(a_C\) for each cells, which is the negative summation of \(a_F\), (3) store the field BCs-related coefficients as internalCoeffs and boundaryCoeffs, respectively. All of these are implemented in gaussLaplacianScheme.C . Note that the Gaussian scheme is the only choice for Laplacian discretization in OF.

Access Laplacian coeffs

Listing 29 Access fvm::Laplacian discretization

fvScalarMatrix Laplacian(fvm::laplacian(DT, T));
Info<<"fvm::laplacian(DT, T): "<<"\n"
   <<"\tLower"<<Laplacian.lower()<<"\n"
   <<"\tDiagonal"<<Laplacian.diag()<<"\n"
   <<"\tUpper"<<Laplacian.upper()<<"\n"
   <<"\tinternalCoeffs"<<Laplacian.internalCoeffs()<<"\n"
   <<"\tboundaryCoeffs"<<Laplacian.boundaryCoeffs()<<"\n"
   <<"\tSource"<<Laplacian.source()<<"\n"
   <<endl;

Tip

The source coefficients come from BCs are not stored in the .source(), but in boundaryCoeffs. So if you print Laplacian.source(), it will display zero. The BCs-related source will be added into the TEqn.source when TEqn.solve() is calling. There is a protected member function named addBoundarySource will be called in solve() function.

Source code of Laplacian

Listing 30 Key source code of fvm::Laplacian (gaussLaplacianScheme.C before nonOrthogonal correcting)

template<class Type, class GType>
tmp<fvMatrix<Type>>
gaussLaplacianScheme<Type, GType>::fvmLaplacianUncorrected
(
   const surfaceScalarField& gammaMagSf,
   const surfaceScalarField& deltaCoeffs,
   const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
   tmp<fvMatrix<Type>> tfvm
   (
      new fvMatrix<Type>
      (
            vf,
            deltaCoeffs.dimensions()*gammaMagSf.dimensions()*vf.dimensions()
      )
   );
   fvMatrix<Type>& fvm = tfvm.ref();

   fvm.upper() = deltaCoeffs.primitiveField()*gammaMagSf.primitiveField();
   fvm.negSumDiag();

   forAll(vf.boundaryField(), patchi)
   {
      const fvPatchField<Type>& pvf = vf.boundaryField()[patchi];
      const fvsPatchScalarField& pGamma = gammaMagSf.boundaryField()[patchi];
      const fvsPatchScalarField& pDeltaCoeffs =
            deltaCoeffs.boundaryField()[patchi];

      if (pvf.coupled())
      {
            fvm.internalCoeffs()[patchi] =
               pGamma*pvf.gradientInternalCoeffs(pDeltaCoeffs);
            fvm.boundaryCoeffs()[patchi] =
               -pGamma*pvf.gradientBoundaryCoeffs(pDeltaCoeffs);
      }
      else
      {
            fvm.internalCoeffs()[patchi] = pGamma*pvf.gradientInternalCoeffs();
            fvm.boundaryCoeffs()[patchi] = -pGamma*pvf.gradientBoundaryCoeffs();
      }
   }

   return tfvm;
}

Laplacian coefficients

Listing 31 fvm::Laplacian Coefficients of the regular mesh shown in Fig. 35.

)
fvm::laplacian(DT, T): 
	Lower
47
(
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
)

	Diagonal
30
(
-8e-07			//0
-1.2e-06
-1.2e-06
-1.2e-06
-1.2e-06
-1.2e-06		//5
-1.2e-06
-1.2e-06
-1.2e-06
-8e-07			//9
-1.2e-06		//10
-1.6e-06
-1.6e-06		//12
-1.6e-06
-1.6e-06
-1.6e-06
-1.6e-06
-1.6e-06
-1.6e-06
-1.2e-06		//19
-8e-07
-1.2e-06
-1.2e-06
-1.2e-06
-1.2e-06
-1.2e-06
-1.2e-06
-1.2e-06
-1.2e-06
-8e-07
)

	Upper
47
(
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
4e-07
)

	internalCoeffs
5
(
3(-8e-07 -8e-07 -8e-07)
3(-8e-07 -8e-07 -8e-07)
10{0}
10{0}
0()
)

	boundaryCoeffs
5
(
3(-0.0002184 -0.0002184 -0.0002184)
3(-0.0004584 -0.0004584 -0.0004584)
10{-0}
10{-0}
0()
)

	Source
30
(
0
1.20371e-36
-1.20371e-36
-1.20371e-36
1.20371e-36
1.20371e-36
1.67048e-52
-1.20371e-36
2.40741e-36
-1.66533e-19
0
0
0
0
0
0
0
0
0
1.66533e-19
0
0
0
0
0
0
0
0
0
0
)

Internal cell

Fig. 42 Information of internal cell (\(C_{12}\))

Boundary cell 19

Fig. 43 Information of boundary cell (\(C_{19}\))

Boundary cell 9

Fig. 44 Information of boundary cell (\(C_{9}\))

Boundary cell 0

Fig. 45 Information of boundary cell (\(C_{0}\))

Boundary cell 10

Fig. 46 Information of boundary cell (\(C_{10}\))

Boundary cell 5

Fig. 47 Information of boundary cell (\(C_{5}\))

Step 3, discretize transient term

discretize Laplacian term

For implicit discretization, fvm::ddt(T) will return a fvMatrix object contains diagonal coefficients and source. The coefficients depend on discretization scheme. For example Euler scheme, the diagonal coefficients are calculated from equation (53) and equation (55). All these are implemented in EulerDdtScheme.C.

Tip

There are 8 schemes for transient discretization in OpenFOAM

1. CoEuler

2. CrankNicolson

3. Euler

4. SLTS

5. backward

6. bounded

7. localEuler

8. steadyState

Access fvm::ddt coeffs

Listing 32 Access fvm::ddt discretization

Info<<"fvm::ddt(T): "<<"\n"
   <<"\tLower"<<ddt.lower()<<"\n"
   <<"\tDiagonal"<<ddt.diag()<<"\n"
   <<"\tUpper"<<ddt.upper()<<"\n"
   <<"\tinternalCoeffs"<<ddt.internalCoeffs()<<"\n" //actually this is not necessary for fvm::ddt, this is definitely equal to zero
   <<"\tboundaryCoeffs"<<ddt.boundaryCoeffs()<<"\n"
   <<"\tSource"<<ddt.source()<<"\n"
   <<endl;

Source code of fvm::ddt

Listing 33 Key source code of fvm::ddt (EulerDdtScheme.C )

template<class Type>
tmp<fvMatrix<Type>>
EulerDdtScheme<Type>::fvmDdt
(
   const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
   tmp<fvMatrix<Type>> tfvm
   (
      new fvMatrix<Type>
      (
            vf,
            vf.dimensions()*dimVol/dimTime
      )
   );

   fvMatrix<Type>& fvm = tfvm.ref();

   scalar rDeltaT = 1.0/mesh().time().deltaTValue(); // 1/dt

   fvm.diag() = rDeltaT*mesh().Vsc(); // Vc/dt (FluxC)

   if (mesh().moving())
   {
      fvm.source() = rDeltaT*vf.oldTime().primitiveField()*mesh().Vsc0();
   }
   else
   {
      fvm.source() = rDeltaT*vf.oldTime().primitiveField()*mesh().Vsc(); // -T_old*Vc/dt
   }

   return tfvm;
}

fvm::ddt coefficients

\(\Delta t = 0.05\ s\)

Listing 34 fvm::ddt Coefficients of the regular mesh shown in Fig. 35.

)fvm::ddt(T): 
	Lower47{0}
	Diagonal
30
(
2e-05
2e-05
2e-05
2e-05
2e-05
2e-05
2e-05
2e-05
2e-05
2e-05
2e-05
2e-05
2e-05
2e-05
2e-05
2e-05
2e-05
2e-05
2e-05
2e-05
2e-05
2e-05
2e-05
2e-05
2e-05
2e-05
2e-05
2e-05
2e-05
2e-05
)

	Upper47{0}
	internalCoeffs
5
(
3{0}
3{0}
10{0}
10{0}
0()
)

	boundaryCoeffs
5
(
3{0}
3{0}
10{0}
10{0}
0()
)

	Source
30
(
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
)

Step 4, construct the final coefficient matrix and RHS

Final matrix

The final coefficient matrix is constructed by simply adding the matrix of Step 2, discretize Laplacian term and Step 3, discretize transient term.

• The diagonal coefficients come from \(a_C\) of Laplacian term and \(FluxC\) of transient term

• The off-diagonal coefficients only come from \(a_F (internal\ face)\) of Laplacian term

• The RHS comes from \(c_F\) of Laplacian term when at boundary faces and \(FluxC^oT_C^o\) of transient term.

Coefficients at the first time step

Listing 35 Final matrix coefficients of the regular mesh shown in Fig. 35 at the first time step.

)TEqn: 
	Lower
47
(
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
)

	Diagonal
30
(
2.08e-05        //0
2.12e-05
2.12e-05
2.12e-05
2.12e-05
2.12e-05        //5
2.12e-05
2.12e-05
2.12e-05
2.08e-05        //9
2.12e-05        //10
2.16e-05
2.16e-05        //12
2.16e-05
2.16e-05
2.16e-05
2.16e-05
2.16e-05
2.16e-05
2.12e-05        //19
2.08e-05
2.12e-05
2.12e-05
2.12e-05
2.12e-05
2.12e-05
2.12e-05
2.12e-05
2.12e-05
2.08e-05
)

	Upper
47
(
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
-4e-07
)

	internalCoeffs
5
(
3(8e-07 8e-07 8e-07)
3(8e-07 8e-07 8e-07)
10{0}
10{0}
0()
)

	boundaryCoeffs
5
(
3(0.0002184 0.0002184 0.0002184)
3(0.0004584 0.0004584 0.0004584)
10{0}
10{0}
0()
)

	Source
30
(
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
0.00546
)

Step 5, solve

For the Regular box case case, we can use PBiCG solver and DILU preconditioner.

Available preconditioner in OpenFOAM

• diagonal : for symmetric & nonsymmetric matrices (not very effective)

• DIC : Diagonal Incomplete Cholesky preconditioner for symmetric matrices

• DILU : Diagonal Incomplete LU preconditioner for nonsymmetric matrices

• FDIC : Fast Diagonal Incomplete Cholesky preconditioner

• GAMG : Geometric Agglomerated algebraic MultiGrid preconditioner

Available solver in OpenFOAM

• BICCG: Diagonal incomplete LU preconditioned BiCG solver

• diagonalSolver: Solver for symmetric and nonsymmetric matrices

• GAMG: Geometric Agglomerated algebraic Multi-Grid solver

• ICC: Incomplete Cholesky Conjugate Gradient solver

• PBiCG: Bi-Conjugate Gradient solver with preconditioner

• PCG: Conjugate Gradient solver with preconditioner

• smoothSolver: Iterative solver with run-time selectable smoother

Krylov-subspace solvers

• CG: The Conjugate Gradient algorithm applies to systems where A is symmetric positive definite (SPD)

• GMRES: The Generalized Minimal RESidual algorithm is the first method to try if A is not SPD.

• BiCG: The BiConjugate Gradient algorithm applies to general linear systems, but the convergence can be quite erratic.

• BiCGstab: The stabilized version of the BiConjugate Gradient algorithm.

Step 6, write

Jupyter notebook

• 1. Read and plot mesh information
• 2. Boundary conditions
• 3. Plot coefficients matrix
  • Implementation in OpenFOAM
• 1. Coordinates and delta coefficients of internal cell
• 2. Coordinates and delta coefficients of boundary cell