over any control volume bounded by the surface on which is the normal directed toward the outside of . The integral conservation equation (1) is valid for any control volume . It is applied to every cell of the mesh. These cells usually consist of rectangles or triangles in 2D, and hexaedra, prisms, or tetrahedra in 3D. More generally, any type of volume, bounded by planar surfaces, is possible.
Hexahedra, Prisms, and Tetrahedra The control volumes (CV) are noted (considering Triangles as generic examples), the union of all CV forming a partition of the entire flow domain. The surface integrals appearing in (1) when applied to the CV is simply the sum of fluxes over the number of planar surfaces separating from the neighbor cells:
The contribution in the summation is the flux of the quantity carried by the flow velocity across the surface between cells and . For cell in figure 1, this writes:
Figure 1. Combination of rectangular and triangular CV The calculation for one interfacial surface between cell and cell needs to be described only once, even though it will appear again when computing the fluxes into cell . Indeed only the direction of the external normal is reversed, thus only the sign of the flux needs to be exchanged. For example on fig 1, the absolute value of flux across the boundary between the quadrangle centered on P and the triangle centered on E must be calculated only once, then it contributes to cell with a + sign if if is pointing out of P ( is pointing to the exterior of cell P). When considering cell E, the same expression will appear except that now is pointing to the exterior of cell E toward cell P hence a - sign is affected to the absolute flux if . This is what makes the method absolutely conservative, for mass, momentum and scalars, however coarse the mesh. It is of course very appealing for industrial calculations where mesh refinement is not always possible. However it is quite difficult to evaluate at the interface, since one knows only the average values and as is explained further. This step is often (hardly!) first order accurate on strongly distorted meshes. (It is usually quite accurate on regular meshes, i.e. second order, but the main reason for using unstructured FV methods is precisely that they can be applied to very complex geometries, where the mesh generators usually produce quite distorted cells) .
Figure 2. 1D convection problem The issues of FV against FD can be first outlined on the following 1D example of convection in a pipe with variable section as shown in figure 2. (3) In 1D the velocity vector is simply where is assumed constant, then if is pointing to the right and otherwise. Let the control volume centered on node and bounded by surfaces , . The FV formulation is then: (4) (5)
This formula would be valid for a pipe with variable section. For simplicity let us assume the sections constant and equal to unity. Then:
note that , while and . Furthermore we have defined the average value of inside the CV as: (6)
One could also introduce the average value of the derivative: (7)
Then the FV formulation can be seen as: (8) This is not an approximation , but an exact representation of the physical process occurring within the segment centered around node . However several approximations will need to the introduced:
and: (11)
i.e. equation (8) is then equivalent to the centered FD: (12)
Although FV methods are currently very appealing for their simplicity and possible use on any type of unstructured grid, these methods still contain discretisation errors, which are less obviously exhibited than the truncation errors of finite differences.
Question: What is the order of the approximation ?
Figure 3. Two possibilities for storage of variables
with for instance at the lowest order: (13)
This however introduces some error since ordinarily (except for regular grids) the surface center is not exactly the middle of cell centers and Note that is a constant geometric factor, computed and stored once for all at the beginning of the calculation.
which is obviously much simpler.
Figure 4. Dual mesh for the cell vertex storage
Actually, in the cell-vertex case, new control volumes (called dual mesh, see figure 4) are often constructed around the cell vertex, onto which the cell-centered approach is then applied. To compute the fluxes along the new CV boundaries, variables need to be interpolated from the cell vertices (where they are stored). There seems to be no difference with the cell centered approach except that we have created a new ”dual mesh”, but since this dual mesh is not given by the mesh generator, but by the FV code, it can be defined in an way that will reduce the error in the approximation of the interpolation of the value on the surface (13). In the figure above for example, the dual mesh is constructed from the centres of the triangles (*) and the middles of segments (+). Gradients and interpolated values at these dual nodes are then easily obtained, but the number of sub elements across which the fluxes must be computed is very large (10 on the above example, compared to the 3 sides of the original triangular mesh).
For triangles and quadrangles (hexa and tetrahedra in 3D) it is convenient to use the well-known shape functions of the finite element method. (14)
where is a polynomial whose value is 1 when and 0 on the other vertices In opposition to traditional finite elements, there is no variational formulation, i.e. no multiplication by shape functions, equation (1) is simply used as it stands, and the shape functions are only used for interpolations. This is called the Finite-Volume/Finite-Element (FV/FE) method.
To sum up, the cell centered FV approach uses the original mesh (provided by the mesh generator) as control volumes and stores variables at their center. The cell vertex FV-FE approach stores the variables on the original mesh nodes, and reconstructs new control volumes around them, in a way that will facilitate (or make more accurate) the calculation of the interface fluxes. A disadvantage of the cell-vertex storage is that it involves a significantly larger number of interfacial fluxes to compute. Also, the stored variables lie exactly on the boundary, which does not simplify the use of e.g. ”wall functions” in turbulent flows where Neuman conditions are imposed rather than Dirichlet conditions. In the following, we consider only cell-centered discretisation.
The gradient at the cell centre, is defined by the Gauss theorem: (17)
where the summation is performed on all interface centers, , (for a rectangle ). The geometric coefficient, , is the area of the interface surface projected on the base vector of the Cartesian frame ( is the unit vector along the axis, the procedure is repeated for then ). Next the volume integral on the left of (17) is evaluated using the value at the CV center: (18)
Thus (19)
If the grid is orthogonal, then for instance is obtained by linear interpolation from the CV centers (20)
using the ratio of lengths Once the gradients at the centers are known using (?), then in turn they can be interpolated to give the gradients at the faces: (21)
The problem, as shown in figure 22, is that for a non orthogonal grid, the interpolation (20) is only correct for point , the intersection of line PE with the east face. Thus is must be corrected by: (22) (23)
Figure 5. center of face e, and intersection e’ Unfortunately, the gradients are not known at this stage. The problem, in order to compute the gradients at the cell CV centers, is that the values of on the cell faces are needed, and in order to obtain these, when the mesh is not orthogonal, the interpolations from the CV centers to the centers of the faces involve the gradients. This ”gradient reconstruction” is an implicit relation that can be solved interatively. After the initial sequence (20), (21), (22), one can iterate over (21) and (22). This can be a very time consuming process. For stability reasons, the diffusion step of the Navier Stokes equation solver is usually implicit, and solved iteratively by a Gauss-Seidel or conjugate gradient method. Within each of these iterations, new values of are obtained, and so new ”gradient reconstruction” steps must be performed. Overall, this is an iterative process, within an iterative process. This is clearly a disadvantage compared to finite elements where the matrices are clearly defined at the beginning, separate from the linear system resolution. Alternatively, on can proceed ”as if” the grid was orthogonal, for the implicit part of the code, and use values of the gradient from the previous time step to correct the non-orthogonality in the interpolation, as explained further.
Indeed, there is yet an additional problem which is that the process of first calculating the gradients at the cell centers, then interpolating them on the cell faces leads to a very large discretisation molecule (a Laplacian calculated on cell P involves all neighbours of cell P, but also all neighbours of neighbours). Using values in all the neighbouring cells of , is computed, then values of all neighbours of cell E are needed to obtain . Next, is given by interpolation from the two previous. This involves a very large molecule as seen on the figure below, and moreover the same procedure must be extended to all faces of P. Also, the iteration between (21) and (22) is not local, restricted to cell P, but must be carried over all cells of the mesh, i.e. compute (21) for all cells, then (22) for all cells, and repeat.
CV involved in the calculation of Instabilities resulting from a too large dicretisation molecule can be illustrated by comparison to the more familiar Finite Differences on a rectangular Cartesian grid. The projection of the north and south interfaces in the direction is zero, while that of the east and west interfaces is simply:
Thus
when the interface values are replaced by the interpolation of the nodal values, this results in:
which is the standard second order central difference. Fortunately no gradient reconstruction iterations are needed here since the grid-lines are orthogonal. In the momentum equation of cell will appear the difference between east and west fluxes (resulting from the volume integrated laplacian):
On the last line we have switched to finite difference notation to show that even nodes are linked only to even nodes, (and odd nodes to only odd). Thus the solution on the odd nodes can become decoupled from that on even nodes (excepting the effect of boundary conditions). This is similar to the checkerboard oscillations that can appear for the discretisation of pressure on a collocated grid. In practice however this mainly appears on structured grids. As soon as a few triangles (or non cubic elements in 3D) are introduced, the solutions tend to be recoupled (imagine a checkerboard where two cells are lumped together). Nevertheless, the above construction results in a broad stencil, that, when reaching the stage linear solving linear system leads to a matrix with many non zero terms. Since we will be interested in solving the diffusion implicitly, for stability reasons, this matrix will be difficult to inverse. For a non regular grids, the above ”pseudo laplacian” results in:
When 2 or 3D problems and non orthogonality are considered, many more points are involved. It is thus desirable to seek a more compact approximation to the fluxes. Note that when the control volumes are either rectangles or equilateral triangles, the centre of mass lies also along the normal directions to the interface centres, so that it is possible to directly compute the gradient along the normal: (24)
where is the distance between CV centres. This can then directly be cast in the flux integral: For the case of regular grids we recover the usual discretisation of the (volume integrated) Laplacian.
This obviously leads to a much more compact discretisation, which moreover is free of odd-even oscillations. The problem is that we are in fact calculating the gradient along the direction connecting the CV centres, instead of . However we can use the compact gradient in the implicit formulation, then correct it explicitly using the more exact gradient obtained by interpolations from the CV centres.
The first term is treated as implicit, while the second one is the so called deferred correction and needs to be explicit since it involves a large ”molecule”. On a Cartesian mesh the correction is nill and the discretisation of the diffusion operator is the same as presented in the finite difference section. For very distorted meshes however the correction becomes larger than the implicit part and this may be a source of instabilities.
These relations can be used to compute and from and . We write the previous equations as a linear system: , with:
We have 3 equations for only 2 unknowns, and . In general has lines and columns, where is the number of neighbouring cells sharing an interface with cell , and is = 2 or 3 in 2D and 3D space respectively. In the 2D triangle example considered here and . For the most common case of hexa cells in 3D (bricks), and . There is thus no unique solution to but a ”least-squares method” can be used. For this we mutiply both sides of the linear system by the transpose :
The gradient components in cell P for any variable are thus obtained by:
is a symmetric matrix. is a matrix that depends only on the mesh geometry. only needs to be computed once and stored at the begining of the simulation. is simply the differences .. . This is a simple direct method, relatively low cost compared to the interative ”Gauss formula” method presented earlier. It is however less accurate on stretched or skewed meshes.
Consider the case of a triangular mesh. The three medians of a triangle intersect at the centre of mass or centroid. _ The three _perpendicular bisectors of the sides of a triangle intersect in a common point, called the circumcentre of the triangle. For the storage of the variables the ideal choice for the calculation of volume integrals is the centroid, while the best choice for the diffusive fluxes is obviously the circumcenter. The deferred correction is the price we need to pay in order to allow arbitrary shapes of control volumes. For equilateral triangles, (or other equilateral polygons) this correction is not needed since centroid and circumcenter are identical. Since distorted cells are ordinarily only necessary near the boundaries, it would be more efficient to mesh most of the domain with regular polygons, and to mark the irregular polygons in order to apply the deferred correction only on those cells. Unfortunately this feature is not yet present in standard meshing packages.