3.4.10. Interpolation Functions and Negative Coefficients

A sufficient condition for a monotonic differencing scheme is that all the off-diagonal terms in the stencil be of opposite sign from the diagonal term [232]. Coefficient sets with mixed signs in the off-diagonal entries can potentially admit oscillatory solutions. In this section, the sign convention is that diagonal elements are negative and off-diagonal elements should be greater than or equal to zero. The term “negative coefficients” refers to one or more negative off-diagonal coefficients. Schemes with positive coefficients are usually considered important only when designing upwind convection operators, but they may be just as important for diffusion operators. Monotonic diffusion operators are most useful for artificial viscosity schemes and projection methods in application to the low Mach number Navier-Stokes equations. Positive coefficients are particularly important for the Poisson equation that arises when calculating a velocity correction to the continuity equation. The computed field for the velocity potential should be smooth so that no oscillations are introduced into the pressure field.

Mixed-sign off-diagonal coefficients commonly arise in finite-element-like methods for describing the diffusion operator. Christie and Hall [233] note that applying the Galerkin finite-element method (GFEM) with bilinear quadrilateral elements to harmonic functions sometimes results in negative coefficients. It was later discovered that there is a threshold element aspect ratio for positivity, and negative coefficients are produced on meshes of rectangular elements above that threshold value. Several authors note that the threshold aspect ratio for the quadrilateral element is $\sqrt{2}$ with GFEM and the value is $\sqrt{3}$ with the control volume finite-element method (CVFEM) [141, 234, 235]. The values for the aspect ratio limits only strictly apply to orthogonal structured meshes. Notably, the five-point difference stencil for the 2D finite-difference method never generates negative coefficients. By deduction, the integral formulas that use extra stencil points, introduced by the element-based methods, generate negative coefficients.

A word on oscillations is required before continuing. Smooth solutions are possible with negative coefficients. Finite-element and finite-volume analysis codes for diffusion processes, such as conduction heat transfer, may never experience oscillations. A forcing function is required to induce the oscillations, like a boundary layer with the convection-diffusion equation [236] or an ill-behaved source term in the continuity equation. The mass balance for a control volume is the source term in the projection method. If the mass balance from a time integration step is particularly bad, the projection scheme must smooth large errors. If there are negative coefficients for the velocity potential, then resulting velocity potential field may be non-smooth, which causes the pressure field to be non-smooth. The result is oscillations which grow into the solution and make for a non-robust solution process.

Causes of negative-coefficients are being studied in order to design robust solution algorithms for the Navier-Stokes equations. The numerical method of primary interest is the CVFEM, though results for the GFEM are included for comparison. The GFEM community describes negative-coefficient effects as “hour-glassing”. The “hour-glass” oscillations are most common to reduced-integration formulations for the diffusion equation, and stabilization methods [130, 237] have been developed to damp the oscillations. In the CVFEM [158, 159], negative coefficients are prevented by shifting the integration points for the diffusion flux formulation out towards the edges of the control volumes and elements. The method is termed “integration point shifting” in this section. There is no general way to control the coefficient signs when skewed quadrilateral elements are used with arbitrary connectivities. Coefficient control is not a panacea for negative coefficients since integration point shifting generally reduces the accuracy. Ultimately, the proper mesh will have no negative coefficients at all.

In this section, the numerics behind negative coefficients are discussed for the diffusion operator given by

(3.1145) $\int { {{\partial \phi} \over {\partial x_i}} n_i {\rm d}S },$

where the surface differential is defined by (3.1135). Integration point shift functions are derived for the CVFEM diffusion operator. Shift functions are presented for both two and three dimensions which guarantee coefficient positivity for a particular element aspect ratio. Also, the integration point shifting for CVFEM is shown to be similar to hour-glass stabilization for GFEM.

3.4.10.1. Positive Coefficients for Orthogonal Meshes

Negative coefficients arise in the off-diagonal coefficients when the aspect ratio of an element becomes large. Consider the elemental flux contributions to the control volume centered about Node 3, shown in Fig. Control volume faces in a single element. Contributions to the the control volume centered about node 3.. The first off-diagonal node to have a negative coefficient is the side node farthest from the control volume center, Node 4. The negative coefficient is associated with the vertical flux over the long horizontal face. At the integration point on the long horizontal face, the flux is approximated by an average of the difference between Nodes 3 and 2 and the difference between Nodes 4 and 1. The weighting between the two differences is determined by the location of the integration point. The negative coefficient is removed by removing the influence from the Node 4–1 difference. The integration point is shifted farther from Nodes 4 and 1, towards Nodes 2 and 3. The integration points are shifted such that the element-level coefficients are positive, a sufficient condition for global positivity.

Fig. 3.42 Control volume faces in a single element. Contributions to the the control volume centered about node 3.

Integration point shift functions and the critical aspect ratio are derived for isoparametric quadrilateral elements with bilinear shape functions, and hexahedral elements with trilinear shape functions. Only the orthogonal form is considered. Linear triangles are discussed since they can also produce negative coefficients. For linear elements, only the element geometry (mesh quality) can be modified to control negative coefficients. With isoparametric bilinear and trilinear elements, both the geometry and the location of the integration point control negativity. In addition, integration point shifting is compared to finite-element hour-glass control. Positive coefficients are achieved by either shifting the element integration points or applying the hour-glass stabilization matrix, and in some cases the two are identical.

3.4.10.1.1. Aspect Ratio Definition

In this section, the isoparametric coordinates for an element are oriented such that the aspect ratio is greater than or equal to one. In two dimensions, the aspect ratio for an orthogonal element is the ratio between edge lengths. In three dimensions, each element has two aspect ratios since there are potentially three different edge lengths. The aspect ratios are taken relative to the shortest of the edge lengths which has a reference length of one.

3.4.10.1.2. Quadrilateral Elements

The coefficients for the diffusion operator result from the combination of two basic second-order accurate diffusion operators: the edge operator and the centroid operator, shown in Figure 3.43. The two operators represent the extremes in evaluating derivatives using the bilinear shape function within the element. The edge scheme always gives positive coefficients while the centroid scheme gives rise to negative coefficients above a certain aspect ratio. The centroid scheme results from evaluating all the derivatives for the four control volume sub-faces at the element centroid, equivalent to reduced integration [237] in GFEM. The edge scheme results from evaluating the derivatives out at the ends of the sub-faces in CVFEM or out at the nodes in GFEM. The edge scheme returns the standard five-point finite-difference operator. The traditional CVFEM [159] uses an equal weighting of the edge and the centroid scheme. The GFEM uses one part edge to two parts centroid.

The GFEM is more prone to oscillations with high aspect ratio elements than the CVFEM because it contains a larger weighting of the centroid scheme. The single-point-integrated GFEM element will be the most unstable since it is a pure centroid scheme.

Fig. 3.43 Flux integration points (X) determine nodal ( $\bullet$ ) contributions to the coefficient stencil: a) mid-face rule of CVFEM, b) *edge*–operator, c) *centroid*–operator (one-point integration).

Fig. 3.43 Flux integration points (X) determine nodal ( $\bullet$ ) contributions to the coefficient stencil: a) mid-face rule of CVFEM, b) *edge*–operator, c) *centroid*–operator (one-point integration).

The coefficient signs for a rectangular element are controlled by moving the integration points away from the centroid of the element. The smallest value of the integration point shift that satisfies coefficient positivity is found by symbolically integrating the diffusion flux over a control volume. Consider the diffusion operator evaluated over a collection of equal-size rectangular elements. Each element is longer by a factor of ${\cal AR}$ in the x-direction than the y-direction, where ${\cal AR}$ is the aspect ratio of the elements. The integration points can be shifted in the $\xi$ -direction by $s$ and in the $\eta$ -direction by $t$ . The element coefficients that contribute to the equation centered at Node 3 in Fig. Control volume faces in a single element. Contributions to the the control volume centered about node 3. are:

(3.1146) $\phi_1: \phantom{-} {{\phantom{+} 1 + \phantom{3} {\cal AR}^2 - 2 s {\cal AR}^2 - 2 t} \over {8 {\cal AR}}}$

(3.1147) $\phi_2: \phantom{-} {{- 1 + 3 {\cal AR}^2 + 2 s {\cal AR}^2 + 2 t } \over { 8 {\cal AR}}}$

(3.1148) $\phi_3: - {{ \phantom{-} 3 + 3 {\cal AR}^2 + 2 s {\cal AR}^2 + 2 t } \over { 8 {\cal AR}}}$

(3.1149) $\phi_4: \phantom{-} {{ \phantom{+} 3 - \phantom{3} {\cal AR}^2 + 2 s {\cal AR}^2 + 2 t } \over { 8 {\cal AR}}}$

The positivity constraint for Node 3 comes from coefficient 4 in the element matrix, where the coefficient becomes negative for large values of aspect ratio. The $s$ -shift removes the effect of aspect ratio, while the $t$ -shift has no effect on negativity. For CVFEM, the integration points on vertical faces, in the longer x-direction, should be shifted from the mid-faces out towards the element edges by $s$ . The y-direction flux is the only flux effected by the shift so the y-direction flux is the flux associated with negative coefficients. The minimal amount of $s$ -shift required to maintain positivity depends upon the aspect ratio,

(3.1150) $s > {1 \over 2} {{{\cal AR}^2 - 3} \over {{\cal AR}^2}}, \quad \quad {\cal AR} > \sqrt{3} .$

The maximum aspect ratio for which the unshifted CVFEM remains monotone is $\sqrt{3}$ . A similar formula exists for GFEM, where $s$ and $t$ are shifted from the Gauss points at $\pm 1/\sqrt{3}$ . The element coefficients that contribute to the equation centered at Node 3 in Fig. Control volume faces in a single element. Contributions to the the control volume centered about node 3. are:

(3.1151) $\phi_1: {{\phantom{+} 3 + 3 {\cal AR}^2 - \left( 1 + \sqrt{3} s \right)^2 {\cal AR}^2 - \left( 1 + \sqrt{3} t \right)^2 } \over {12 {\cal AR}}}$

(3.1152) $\phi_2: {{ - 3 + 3 {\cal AR}^2 + \left( 1 + \sqrt{3} s \right)^2 {\cal AR}^2 + \left( 1 + \sqrt{3} t \right)^2 } \over {12 {\cal AR}}}$

(3.1153) $\phi_3: {{ - 3 - 3 {\cal AR}^2 - \left( 1 + \sqrt{3} s \right)^2 {\cal AR}^2 - \left( 1 + \sqrt{3} t \right)^2 } \over { 12 {\cal AR}}}$

(3.1154) $\phi_4: {{ \phantom{+} 3 - 3 {\cal AR}^2 + \left( 1 + \sqrt{3} s \right)^2 {\cal AR}^2 + \left( 1 + \sqrt{3} t \right)^2 } \over { 12 {\cal AR}}}$

Similar to the CVFEM, the $s$ -shift is the only shift that affects the aspect ratio term. The positivity constraint for Node 3 comes from coefficient 4 in the element matrix,

(3.1155) $s > \sqrt{{{3 {\cal AR}^2 - 4} \over {3 {\cal AR}^2 } }} - {1 \over \sqrt{3}}, \quad \quad {\cal AR} > \sqrt{2} .$

The maximum aspect ratio for which the unshifted GFEM remains monotone is $\sqrt{2}$ .

The shift values are very sensitive to the aspect ratio, out to an aspect ratio of about four. The values of the integration point shift function are plotted as a function of the aspect ratio in Figure 3.44. At that aspect ratio, the shifted integration points are near the edge of the element. Since the requisite integration point shift rapidly reaches the element edge for increasing aspect ratio, it can be argued that the maximum shift should always be taken. It will be shown in the section on accuracy that integration point shifting leads to a general loss in accuracy on non-orthogonal meshes. Therefore, it may be desirable to use (3.1150) to compute the minimal shift for each element. The accuracy consideration must be traded against the algorithmic complexity of computing geometry-dependent shape functions for each element.

Fig. 3.44 Integration point locations must be shifted out towards the element-edge with increasing element aspect ratio.

3.4.10.1.3. Reduced Integration

One-point integration methods for quadrilateral and hexahedral elements are popular because they are computationally efficient. Sometimes, oscillations occur and are called hour-glass modes after the displaced element shapes. Hour-glass stabilization methods prevent the hour-glass oscillations from occurring. The hour-glass terms are derived by examining the eigenmodes of the finite-elements and noting that there are missing mode shapes when the elements are integrated at the center [237]. The stabilization term adds the effects of the missing mode shapes back into the element formulation.

It is shown here that the hour-glass stabilization terms have the same effect on the element coefficients as shifting the integration points in two-dimensional elements. Both the hour-glass stabilization and the integration-point shifting modify the discretization to look more like a five-point scheme; or, more like the edge scheme of the previous section. The H-stabilization method is commonly used [130] in the GFEM with reduced integration. For quadrilateral elements, the element matrix for the hour-glass stabilization term is

(3.1156) $C_{\rm hg} \left[ \begin{array}{rrrr} -1 1 -1 1 \\ 1 -1 1 -1 \\ -1 1 -1 1 \\ 1 -1 1 -1 \end{array} \right],$

where the constant $C_{\rm hg}$ contains scaling information. The first row of the hour-glass stabilization matrix, (3.1156), can also be derived by subtracting the matrix coefficients for the one-point integrated GFEM, (3.1157),

(3.1157) $\begin{array}{lc} \phi_1: & \phantom{-}{{1 + {\cal AR}^2}\over { 4 {\cal AR}}} \\ \phi_3: & -{{ 1 + {\cal AR}^2} \over { 4 {\cal AR}}} \\ \phi_2: & \phantom{-}{{{\cal AR}^2 - 1} \over { 4 {\cal AR}}} \\ \phi_4: & \phantom{-}{{ 1 - {\cal AR}^2} \over { 4 {\cal AR}}} \end{array}$

from the coefficients for the five-point difference scheme, (3.1158),

(3.1158) $\begin{array}{lc} \phi_1: & \phantom{-} 0 \\ \phi_3: & - {{ 1 + {\cal AR}^2} \over { 2 {\cal AR}}} \\ \phi_2: & \phantom{-} {{ {\cal AR}^2} \over { 2 {\cal AR}}} \\ \phi_4: & \phantom{-} {{ 1 } \over { 2 {\cal AR}}} \end{array}$

The resulting coefficient set is identical to the hour-glass stabilization matrix, (3.1156), if the multiplier $C_{\rm hg} = (1+{\cal AR}^2)/4{\cal AR}$ . The hour-glass stabilization matrix is added to a diffusion operator to make it look more like a five-point finite difference scheme. Note that the multiplier used by GFEM practitioners [130] is $C_{\rm hg} = 1$ . The hour-glass stabilization matrix can also be used with other schemes. The multiplier for standard GFEM is $C_{\rm hg} = (1+{\cal AR}^2)/6{\cal AR}$ . The multiplier for standard CVFEM is $C_{\rm hg} = (1+{\cal AR}^2)/8{\cal AR}$ , though conservation is only guaranteed on rectangular meshes.

3.4.10.1.4. Triangular Elements

Linear triangular elements can also produce negative off-diagonal coefficients.

There are no shift-functions for triangles since the gradients are constant over the element. Negative coefficients result from the geometry of the element.

Nodal coefficients for a triangular element are computed for the diffusion flux contributions, shown in Figure 3.45. The nodal coefficients for the diffusion flux contribution to the control volume centered about Node 1 are

(3.1159) $\phi_1: \phantom{x} - {1 \over 2} {{\tan \alpha} \over {\tan \beta}} {{\left(1 + \tan^2\beta\right)} \over{\left(\tan \alpha + \tan \beta \right)}}$

(3.1160) $\phi_2: \phantom{x} \phantom{-} {1 \over 2} {{\left( \tan \alpha \tan \beta - 1 \right)} \over {\left(\tan \alpha + \tan \beta \right)}}$

(3.1161) $\phi_3: \phantom{x} \phantom{-} {1 \over 2} {1 \over {\tan \beta}}$

where the base edge length is $r$ and the two adjacent vertex angles are $\alpha$ and $\beta$ . The conditions required to ensure that all the off-diagonal terms remain positive are combined from constraints on all three control volume contributions,

(3.1162) $\tan \alpha > 0$

(3.1163) $\tan \beta > 0$

(3.1164) $\tan \alpha \tan \beta > 1.$

The triangular element yields positive off-diagonal coefficients if $0 < \alpha < \pi/2$ , $0 < \beta < \pi/2$ , and $\pi/2 < \alpha + \beta < \pi$ . The triangle must be acute.

Fig. 3.45 Triangular element geometry in defined by edge length and two vertex angles.

In a previous work [204], quadrilateral elements were subdivided into triangular elements using edge-swapping, along with a Delaunay algorithm, to minimize the effect of negative coefficients. Given a collection of nodes, a Delaunay triangulation of the nodes will generate triangles where the minimum angle between vertices is maximized, leading to near-equilateral triangles that satisfy (3.1162).

3.4.10.1.5. Hex Elements

The aspect ratio limit for the CVFEM diffusion operator with orthogonal, hexahedral elements can be as large as $\sqrt{2}$ if the base is square. The three-dimensional element is more difficult to control because there are two aspect ratios. The local element node numbering for a hexahedron is defined in Fig. Hexahedral element topology and numbering. Let the edge length between Nodes 1 and 2 be of value $A$ , the edge length between Nodes 1 and 4 be of value $B$ , and the edge length between Nodes 1 and 5 be of value $C$ . These are the $\xi$ , $\eta$ , and $\zeta$ directions. The parametric representation of the nodal coefficients for the element contribution to the control volume centered about Node 1 is

(3.1165) $\phi_1: \phantom{x} - {{9 \left( \phantom{-} A^2 B^2 + \phantom{3} B^2 C^2 + \phantom{3} C^2 A^2 \right) } \over {64 A B C}}$

(3.1166) $\phi_2: \phantom{x} \phantom{-} {{3 \left( - A^2 B^2 + 3 B^2 C^2 - \phantom{3} C^2 A^2 \right) } \over {64 A B C}}$

(3.1167) $\phi_3: \phantom{x} \phantom{-} {{3 \left( - A^2 B^2 + 3 B^2 C^2 + 3 C^2 A^2 \right) } \over {64 A B C}}$

(3.1168) $\phi_4: \phantom{x} \phantom{-} {{3 \left( - A^2 B^2 - \phantom{3} B^2 C^2 + 3 C^2 A^2 \right) } \over {64 A B C}}$

(3.1169) $\phi_5: \phantom{x} \phantom{-} {{\phantom{3} \left( 3 A^2 B^2 - \phantom{3} B^2 C^2 - \phantom{3} C^2 A^2 \right) } \over {64 A B C}}$

(3.1170) $\phi_6: \phantom{x} \phantom{-} {{\phantom{3} \left( 3 A^2 B^2 + 3 B^2 C^2 - \phantom{3} C^2 A^2 \right) } \over {64 A B C}}$

(3.1171) $\phi_7: \phantom{x} \phantom{-} {{\phantom{3} \left( \phantom{-} A^2 B^2 + \phantom{3} B^2 C^2 + \phantom{3} C^2 A^2 \right) } \over {64 A B C}}$

(3.1172) $\phi_8: \phantom{x} \phantom{-} {{\phantom{3} \left( 3 A^2 B^2 - \phantom{3} B^2 C^2 + 3 C^2 A^2 \right) } \over {64 A B C}} .$

The region for positive coefficients is plotted in Fig. Limits of Edge-Length Ratio for Positive Coefficients in 3D CVFEM. which comes from examining coefficients for Nodes 2, 4, and 5. The maximum allowable aspect ratios occur for the case of a square base. If the base edges are longer than the vertical edge, then the maximum aspect ratio is $\sqrt{2}$ . If the vertical edge is longer than a base edge, the maximum aspect ratio is $\sqrt{3/2}$ .

The integration points in the CVFEM scheme can be shifted by $s$ , $t$ , and $u$ in the $\xi$ , $\eta$ , and $\zeta$ directions. The conditions for positive coefficients are

(3.1173) ${1 \over {A^2}} {1 \over {1 - 2 s}} - {1 \over {B^2}} {1 \over {3 + 2 t}} > {1 \over {C^2}} {1 \over {3 + 2 u}}$

(3.1174) ${1 \over {B^2}} {1 \over {1 - 2 t}} - {1 \over {C^2}} {1 \over {3 + 2 u}} > {1 \over {A^2}} {1 \over {3 + 2 s}}$

(3.1175) ${1 \over {C^2}} {1 \over {1 - 2 u}} - {1 \over {A^2}} {1 \over {3 + 2 s}} > {1 \over {B^2}} {1 \over {3 + 2 t}}$

The relations are nonlinear and require iteration to extract the limiting values of $s$ , $t$ , and $u$ .

The coefficients that are generated by the standard GFEM operator in three dimensions have no allowable maximum aspect ratio. The only element shape that does not have negative coefficients is a cube, and even then the coefficients of the six nearest nodes are zero. The parametric representation of the nodal coefficients for the element contribution to the equation associated with Node 1 are

(3.1176) $\phi_1: \phantom{x} - {{4 \left( \phantom{-} A^2 B^2 + \phantom{2} B^2 C^2 + \phantom{2} C^2 A^2 \right) } \over {36 A B C}}$

(3.1177) $\phi_2: \phantom{x} \phantom{-} {{2 \left( - A^2 B^2 + 2 B^2 C^2 - \phantom{2} C^2 A^2 \right) } \over {36 A B C}}$

(3.1178) $\phi_3: \phantom{x} \phantom{-} {{\phantom{2} \left( - A^2 B^2 + 2 B^2 C^2 + 2 C^2 A^2 \right) } \over {36 A B C}}$

(3.1179) $\phi_4: \phantom{x} \phantom{-} {{2 \left( - A^2 B^2 - \phantom{2} B^2 C^2 + 2 C^2 A^2 \right) } \over {36 A B C}}$

(3.1180) $\phi_5: \phantom{x} \phantom{-} {{2 \left( 2 A^2 B^2 - \phantom{2} B^2 C^2 - \phantom{2} C^2 A^2 \right) } \over {36 A B C}}$

(3.1181) $\phi_6: \phantom{x} \phantom{-} {{\phantom{2} \left( 2 A^2 B^2 + 2 B^2 C^2 - \phantom{2} C^2 A^2 \right) } \over {36 A B C}}$

(3.1182) $\phi_7: \phantom{x} \phantom{-} {{\phantom{2} \left( \phantom{2} A^2 B^2 + \phantom{2} B^2 C^2 + \phantom{2} C^2 A^2 \right) } \over {36 A B C}}$

(3.1183) $\phi_8: \phantom{x} \phantom{-} {{\phantom{2} \left( 2 A^2 B^2 - \phantom{2} B^2 C^2 + 2 C^2 A^2 \right) } \over {36 A B C}} .$

The contributions from Nodes 2, 4, and 5 are such that there will be negative coefficients for any element shape other than a cube.

Single-point integration is also used for the GFEM hexahedral element. Unfortunately, there is not a clear analogy between the integration point shift and hour-glass stabilization in three dimensions. The element matrix for the hour-glass stabilization term is

(3.1184) $\!\!\!C_{\rm hg} \!\left[\! \begin{array}{rrrrrrrr} -4\! 2\! 0\! 2\! 2\! 0\! -2\! 0\! \\ 2\! -4\! 2\! 0\! 0\! 2\! 0\! -2\! \\ 0\! 2\! -4\! 2\! -2\! 0\! 2\! 0\! \\ 2\! 0\! 2\! -4\! 0\! -2\! 0\! 2\! \\ 2\! 0\! -2\! 0\! -4\! 2\! 0\! 2\! \\ 0\! 2\! 0\! -2\! 2\! -4\! 2\! 0\! \\ -2\! 0\! 2\! 0\! 0\! 2\! -4\! 2\! \\ 0\! -2\! 0\! 2\! 2\! 0\! 2\! -4\! \end{array} \right],$

where the constant $C_{\rm hg}$ contains scaling information. The hour-glass stabilization matrix does act to make the diagonal terms more negative, and the problematic off-diagonal terms more positive.