3.4.2. Flow Solver

The core flow solver is based on a segregated, projection method approach. The projection method is used to compute the pressure field which is consistent with a velocity field that satisfies continuity. A pressure-smoothing approach similar to the Rhie/Chow scheme [171] is used to prevent pressure decoupling on the collocated mesh. An upwind method is used to interpolate convected values to control volume faces. Detailed descriptions of these methods are discussed in the following sections.

Another prevalent CVFEM method in the literature is the FIELDS method [150, 151]. The continuity and momentum equations are fully coupled in this approach. We experimented with this approach and found that the three-dimensional discrete equations were difficult to solve and open boundary conditions difficult to implement.

3.4.2.1. Projection Method

The role of pressure smoothing, or explicit stabilization, was first developed in the context of collocated finite volume schemes by [171]. Although this original paper did not explore the formal error introduced by this explicit stabilization, [214] later displayed the sensitivity of steady results on relaxation parameters and provided a methodology to circumvent this issue. In general, such early papers (cf. [215]) as well as other more recent papers, (cf. [216]) introduced the role of stabilization almost by happenstance as it entered only through the specific choice of the convecting velocity formula, i.e., the integration point velocity that forms the mass flow rate.

Studies of [217] and [218], each in the context of a finite element algorithm, have commented on the role of stabilization that is provided by the approximation of the derived pressure correction system, namely that {\bf L} \neq {\bf D} {\bf G}, where {\bf L} is the given discrete Laplacian operator and {\bf D} and {\bf G} are the chosen discrete divergence and gradient operators, respectively. Numerical algorithms for which the Laplacian operator does not equal the discrete divergence of gradient operator have been termed “approximate projection” algorithms (cf. [219] and [220]) in the context of solenoidal flow; in general for non-solenoidal flow the formalism of the projection derivation results in an affine projection.

Recent work by Sandia National Laboratories has cast the general approximate projection algorithm within a family of smoothing and time scaling choices. The analysis of choice that has been followed is to cast the algorithm in terms of an approximate factorization (cf. [221]), and note the added stabilization (herein also known as pressure smoothing terms), and splitting errors. This analysis has been extremely useful in understanding the formal accuracy, and even consistency, of a given numerical scheme.

The analysis of a given computational fluids algorithm in the context of an approximate factorization begins with the discrete momentum and continuity equations written in matrix form. The matrix {\bf A} contains discrete, linearized contributions to the momentum equations from the time derivative, convection, and diffusion terms,

(3.738)\begin{bmatrix} {\bf A} & {\bf G} \\ {\bf D} & {\bf 0} \end{bmatrix} \begin{bmatrix} {\bf u}^{n+1} \\ p^{n+1/2} \end{bmatrix} = \begin{bmatrix} {\bf f} \\ {\bf b} \end{bmatrix} .

The discrete nodal gradient and nodal divergence are {\bf G} and {\bf D} respectively (note that the operator {\bf D} may include aspects of the algorithm due to a variable density field). The function {\bf f} contains the additional terms for the momentum equations, e.g., body force terms, lagged stress tensor terms, etc., while the function {\bf b} contains the appropriate terms for a non-solenoidal velocity field, i.e., - \int{ {\partial \rho } \over {\partial t}}dV. The pressure is appropriately interpreted as the pressure at the n+{1\over2} step, (cf. [222]). The form of the matrix operators can be found in the body of literature for control-volume finite element methods (cf. [157]).

Note that (3.738) is not really complete as the boundary condition values are omitted, however, they are not essential in describing the bulk of the splitting and stabilization analysis as noted by [223]. The boundary conditions would simply enter through an additional vector on the right-hand side and modified entries in the matrix operators.

The approximate factorization of (3.738) takes the general form of

(3.739)\begin{bmatrix} {\bf A} & {\bf 0} \\ {\bf D} & {\bf B_1} \end{bmatrix} \begin{bmatrix} {\bf I} & {\bf B_2 G} \\ {\bf 0} & {\bf I} \end{bmatrix} = \begin{bmatrix} {\bf A} & {\bf A B_2 G} \\ {\bf D} & {\bf B_1 + D B_2 G} \end{bmatrix} .

The factor {\bf B_2} determines the projection time scale. The factor {\bf B_1} defines the linear system for pressure. Ideally, {\bf B_1} could be selected to cancel splitting errors in the continuity equation. Practically, the form of {\bf B_1} is governed by implementation and linear solver efficiency.

A completely generalized set of incremental pressure projection methods with potential stabilization can be written by formally defining the operators {\bf B_1} and {\bf B_2} above, here shown as part of the sequence of equations solved,

(3.740){\bf A} \Delta \hat {\bf u} = {\bf f} - {\bf G} p^{n-{1\over2}} -{\bf A}{\bf u}^n,

(3.741)-{\bf L_1} \Delta p^{n+{1\over2}} = -{\bf D}\left(\hat{\bf u}+ \tilde \tau_2 {\bf G}p^{n-{1\over 2}} \right) + {\bf L_2}p^{n-{1\over2}}+ b,

(3.742){\bf u}^{n+1} = \hat{\bf u} - \tilde \tau_3{\bf G}\Delta p^{n+{1\over2}}.

Laplacian operators acting on a general scalar \phi, which define the approximate nature of the projection method, are given by,

(3.743){\bf L}_1 \phi = \tau_1 \nabla \phi \cdot d{\bf A},

(3.744){\bf L_2} \phi = \tau_2 \nabla \phi \cdot d{\bf A}.

For an approximate projection method,

(3.745){\bf L_2} \neq {\bf D} \tilde \tau_2 \bf G,

while for an exact projection,

(3.746){\bf L_2} ={\bf D} \tilde \tau_2 \bf G.

Exact projections can be easily constructed on unstructured collocated meshes (cf. [224]), although classically this results in a wide Laplacian stencil that admits pressure oscillations yet does not add discrete errors in the continuity solve. We assume that \tau_i factors defined above are represented by a diagonal matrix that corresponds to a particular time scale of choice. The relationship between \tilde \tau_i and \tau_i is normalization by a density and volume,

(3.747)\tilde \tau_i = {{\tau_i} \over {\rho V}}.

The choice of these scaling factors defines the scheme in terms of both stabilization and projection scaling. For example, the ideal form for \tilde \tau_3 is the inverse of {\bf A}. The exact choice of \tilde \tau_3 in a practical sense affects the stability of the scheme. The stabilization terms are represented by operators including both \tilde \tau_1 and \tilde \tau_2 that are required to prevent velocity and pressure decoupling in schemes for which {\bf L} \neq {\bf D} {\bf G}.

Rearrangement of (3.742), in terms of \hat {\bf u}, and substitution of this modified equation into (3.740) and (3.741) provides the full set of splitting and stabilization errors:

(3.748)\left[ \begin{array}{lr} {\bf A} {\bf G} \\ {\bf D} {\bf 0} \end{array} \right] \left[ \begin{array}{l} {\bf u}^{n+1} \\ p^{n+1/2} \end{array} \right] = \left[ \begin{array}{l} {\bf f} \\ {\bf b} \end{array} \right] + \left[ \begin{array}{l} \left({\bf I} -{\bf A } \tilde \tau_3\right){\bf G} \Delta p^{n+{1\over2}} \\ ({\bf L_1}-{\bf D}\tilde \tau_3{\bf G})\Delta p^{n+{1\over 2}}+ ({\bf L_2}-{\bf D}\tilde \tau_2{\bf G})p^{n-{1\over 2}} \end{array} \right] .

The error appearing in the momentum equation is due to splitting and generally can be repaired by non-linear iteration, although ideally single iteration methods are desired (as shown).

Again it is emphasized that for approximate projection methods, {\bf L_2} \neq {\bf D}\tilde \tau_2{\bf G}, whereas for exact projection methods, which are usually based on staggering velocity and pressure, {\bf L_2} = {\bf D}\tilde \tau_2{\bf G} and there is no stabilization error (as there is no need to provide stabilization). Frequently, the stabilization terms within (3.741) are included in a modified provisional velocity (cf. [225]), i.e., \tilde {\bf u} =\hat{\bf u}+ \tilde \tau_2 {\bf G}p^{n-{1\over 2}}, that can often hide the true role of stabilization.

A similar analysis for pressure free projection methods (cf. [226]) can be carried out, in which case the equations solved are given by,

(3.749){\bf A} \hat {\bf u} = {\bf f} -{\bf A}{\bf u}^n,

(3.750)-{\bf L_1} \Delta \phi^{n+1} = -{\bf D} \hat{\bf u} + {\bf L_1}\phi^n+ b ,

(3.751){\bf u}^{n+1} = \hat{\bf u} - \tilde \tau_1{\bf G}\phi^{n+1},

with errors,

(3.752)\left[ \begin{array}{lr} {\bf A} {\bf G} \\ {\bf D} {\bf 0} \end{array} \right] \left[ \begin{array}{l} {\bf u}^{n+1} \\ p^{n+1/2} \end{array} \right] = \left[ \begin{array}{l} {\bf f} \\ {\bf b} \end{array} \right] + \left[ \begin{array}{l} -{\bf A } \tilde \tau_3{\bf G}\phi^{n+1} + {\bf G}p^{n+{1\over2}} \\ \left({\bf L_1} - {\bf D}\tilde \tau_1{\bf G}\right)\phi^{n+1} \end{array} \right] .

The error term in the continuity equation is retained to emphasize that this algorithm can be considered in the context of an approximate projection method. Assuming that the Laplacian and gradient operators commute, it is necessary to compute p^{n+1/2} = {\bf A}\tau_3 \phi^{n+1} to obtain the second-order pressure field, while the relationship p^{n+{1\over2}} = \phi^{n+1} will result in a first-order pressure field with splitting error ({\bf I} - {\bf A}\tau_3){\bf G}p^{n+{1\over2}} ( [227]).

Although the above set of algorithms have been written in terms of a two step scheme, i.e., predict \bf {\hat u} and correct \bf {\hat u} by the appropriately scaled scalar gradient, non-linear iterations can also be taken. In this case, the \phi^{n+1} and {\bf u}^{n+1} state are replaced with the k+1 state, whereas the n+{1 \over 2} pressure state is replaced by the k+{1\over2} state. For the residual form, the n^{th} state is replaced with the current iterate, k^{th} state. At convergence within the time step, \phi^{n+1} = \phi^{k+1}, {\bf u}^{n+1} = {\bf u}^{k+1}, and p^{n+{1\over2}} = p^{k+{1\over2}}.

3.4.2.1.1. CVFEM operators

SIERRA/Fuego uses the finite volume technique known as the control volume finite element method of [152]. Control volumes (the mesh dual) are constructed about the nodes, as shown in Fig. Numerics. Each element contains a set of subfaces that define control-volume surfaces. The subfaces consist of line segments (2-D) or surfaces (3-D). The 2-D segments are connected between the element centroid and the edge centroids. The 3-D surfaces are connected between the element centroid, the element face centroids, and the edge centroids. Integration points also exist within the subcontrol volume centroids. Such integration points are used for volume integrals such as source terms, the mass matrix, and, if chosen, gradients.

Defining \phi_K to be the value of \phi at node K, then the variation of \phi within an element that contains the point location {bf x} is given by

(3.753)\phi({\bf x}) = \sum_{K \in {\mathcal N}} N_{K}({\bf x}) \phi_K ,

where N_K({\bf x}) is the shape function associated with node K at position {\bf x}, and {\mathcal N} is the set of all nodes that defines the element. For the CVFEM, either trilinear (3-D) or bilinear (2-D) shape functions are used. Currently, Fuego supports heterogeneous element topologies consisting of hex, tet, pyramid, and wedges.

The discrete nodal gradient operator for direction i can be written as a surface integral on control volume L,

(3.754){\bf G}\phi = (G \phi)_{Li} = \int_{\Gamma_L} \phi ({\bf x}) dn_i \approx \sum_{\alpha \in {\mathcal B}_L} \left(\sum_{K \in {\mathcal N}} N_{K}({\bf x}_{\alpha}) \phi_K\right) n_i ({\bf x}_{\alpha})\Delta A_{\alpha},

where {\mathcal B}_L is the set of surface integration points for control volume L. Similarly, the discrete divergence operator at node L acting on vector u_i is

(3.755){\bf D} {\bf u} = (D u_i)_{L} = \int_{\Gamma_L} \rho({\bf x}) u_i ({\bf x}) dn_i \approx \sum_{\alpha \in {\mathcal B}_L}\rho({\bf x}_\alpha)\left(\sum_{K \in {\mathcal N}} N_{K}({\bf x}_{\alpha}) u_{Ki} \right) n_i ({\bf x}_{\alpha}) \Delta A_{\alpha} ,

and the Laplacian operator that includes spatially varying timescale, \tau, is

(3.756){\bf L}_\tau \phi = (L_\tau \phi)_{L} = \int_{\Gamma_L} \tau({\bf x}) { {\partial \phi} \over {\partial x_j} } dn_j \approx \sum_{\alpha \in {\mathcal B}_L} \tau({\bf x}_\alpha)\left(\sum_{K \in {\mathcal N}} {{\partial N_K ({\bf x}_\alpha)} \over {\partial x_j}} \phi_ K\right) n_j ({\bf x}_{\alpha}) \Delta A_{\alpha} .

Note that an alternative to the gradient operator given in (3.754), which is provided via the CVFEM is

(3.757){\bf G} \phi = ( G \phi)_{Li} = \int_{\Gamma_L} { {\partial \phi} \over {\partial x_i} } dV \approx \sum_{\alpha' \in {\mathcal B}_L} \left(\sum_{K \in {\mathcal N}} {{\partial N_K ({\bf x}_{\alpha'})} \over {\partial x_i}} \phi_K\right) dV_{\alpha'} ,

where {\mathcal B}_L is now the set of all subcontrol volume integration points for control volume L (for clarity, \alpha' denotes the subcontrol volume integration point location).

The general term {\bf D} \tilde \tau {\bf G} \phi deserves a special note in the case of variable density flows. Specifically, the interpolation is currently provided by the following equation:

(3.758){\bf D} \tilde \tau_i {\bf G} \phi = \sum_{\alpha \in {\mathcal B}_L}\rho({\bf x}_\alpha) { {\tilde \tau_i({\bf x}_\alpha)} \over {\rho({\bf x}_\alpha)}} \left(\sum_{K \in {\mathcal N}} N_{K}({\bf x}_{\alpha}) { {G_{Ki}} \over {V_K}} \right) n_i ({\bf x}_{\alpha}) \Delta A_{\alpha},

(3.759)= \sum_{\alpha \in {\mathcal B}_L} \tilde \tau_i({\bf x}_\alpha) \left(\sum_{K \in {\mathcal N}} N_{K}({\bf x}_{\alpha}) { {G_{Ki}} \over {V_K}} \right) n_i ({\bf x}_{\alpha}) \Delta A_{\alpha}.