3-Numerical Methods

3 Numerical Methods

Overview

CMS-Flow has both implicit and explicit solution schemes. The explicit solver is designed for dynamic problems with extensive wetting and drying that require small computational time steps, while the implicit solver is intended for simulating tidal- and wave-induced circulation at tidal inlets, navigation channels, and adjacent beaches. A detailed description of the numerical formulation of the explicit solver of CMS-Flow can be found in Buttolph et al. (2006) and is not repeated here. The sections below specifically refer to the implicit solver of CMS-Flow.

The implicit solver uses the SIMPLEC algorithm (Van Doormal and Raithby 1984) on a non-staggered grid to handle the coupling of water level and velocity. Primary variables u-, v-velocity, and water level are stored on the same set of grid points, and fluxes at cell faces are determined using a Rhie and Chow (1983) type momentum interpolation method (Wu et al. 2011). The explicit solver uses a staggered grid with velocities at the cell faces and the water levels and water depths at the cell centers (Buttolph et al. 2006a). CMS-Flow also calculates salinity, sediment transport, and morphology change.

CMS-Flow Computational Grid

The implicit version of CMS-Flow uses a generic Cartesian grid which can be regular, nonuniform, or locally refined by splitting a cell into four subcells. Three requirements are necessary for the input grid:

1. Cells must have a rectangular shape. Irregularly shaped cells are not allowed.

2. Cells may have a total of four to six neighboring cells (faces).

3. Only two neighboring cells are allowed in the same direction (i.e. North, South, East, West).

Mesh refinement can be achieved either by locally decreasing the grid spacing (nonuniform Cartesian grid) as shown in the left panel of Figure 3.1 or by subdividing or splitting a cell into multiple cells as shown in the right panel of Figure 3.1. The refined mesh can be further refined and split in multiple levels, as needed. This mesh is referred to here as a telescoping mesh.

Figure 3.1. Examples of a Cartesian grids allowed in CMS: Stretched Cartesian (left) and telescoping grid (right)

One important aspect of incompressible flow models is the location of primary variables: velocity and pressure (water level). On a staggered grid, the pressure (water level) is located at the center of cells and the u- and v-velocities are located along the faces of cells (Harlow and Welsh 1965; Patankar 1980). On a non-staggered grid, all of the primary variables are located at the center of cells. A staggered grid can more easily eliminate the checkerboard oscillations when compared to a non-staggered grid; however, a non-staggered grid involves a simpler source code and can minimize the number of coefficients that must be computed and stored during a simulation because many of the terms in the equations are approximately equal. In particular, a staggered grid is more complicated when handling the interface between coarse and fine cells where five- or six-face control volumes are used. Therefore, a non-staggered (collocated) grid approach is adopted for CMS-Flow, with a Rhie and Chow (1983) momentum interpolation technique used to eliminate the checkerboard oscillations. Figure 3.1 shows the location of primary variables and the 5- and 7-point stencils (computational molecule) used in the calculations for CMS-Flow.

The data structure for a grid can be approached in three ways: 1) block-structured, 2) hierarchical tree, and 3) unstructured. The block-structured approach divides the domain into multiple blocks and each block is treated as structured. A special treatment is applied between blocks to ensure mass and momentum balance using this approach. The hierarchical tree approach is memory intensive and requires parent-child relationships and a tree traverse to determine mesh connectivity. For the unstructured approach, all cells are numbered in a one-dimensional sequence and tables are used to determine the connectivity of neighboring cells. Among these three approaches, the unstructured approach is the most simple and is therefore applied in CMS-Flow.

Computational cells are numbered in an unstructured manner via a 1-D index array. Inactive cells (permanently dry) are not included in the 1-D index array to save memory and computational time. All active computa-tional nodes are numbered sequentially. For convenience with handling boundary conditions, each boundary cell has a neighboring ghost cell outside of the computational domain. Each ghost cell corresponds to a boundary face of the boundary cell. The ghost cells are stored at the end of the 1-D index array.

General Transport Equation

In order to avoid redundant derivations of discretized equations, discretization of a general transport equation is outlined and described below. All of the governing equations are some form of this transport equation; therefore, the same discretization may be applied to all of the governing equations. The general transport equation is given by

$\underbrace {\frac {\partial (h\phi )}{\partial t}} _{\text{Temporal Term}}+\underbrace {\frac {\partial (hV_{j}\phi )}{\partial x_{j}}} _{\text{Advection Term}}=\underbrace {{\frac {\partial }{\partial x_{j}}}\left(\Gamma h{\frac {\partial \phi }{\partial x_{j}}}\right)} _{\text{Diffusion Term}}+\underbrace {S_{\phi }} _{\text{Source Term}}$(3-1)

where $\phi$ is a general scalar, t is time, h is the total water depth, $V_{j}$ is the transport velocity, $\Gamma$ is the diffusion coefficient for $\phi$ , and $S_{\phi }$ includes all remaining terms. Note that in the case of the continuity and momentum equations, $\phi$ is equal to 1 and $U_{j}$ , respectively.

Spatial Discretization

A control-volume technique is used in which the governing equations are integrated over a control volume to obtain an algebraic equation that can be solved numerically. Integration of Equation (3-1) over a control volume yields:

$\int _{A}{\frac {\partial (h\phi )}{\partial t}}dA+\int _{A}{\frac {\partial }{\partial x_{j}}}\left(hV_{j}\phi -\Gamma h{\frac {\partial \phi }{\partial x_{j}}}\right)dA=\int _{A}S_{\phi }dA$(3-2)

${\frac {\partial (h_{P}\phi _{P})}{\partial t}}\Delta A_{P}+\int _{L}h[({\hat {n_{i}}}V_{i})\phi -\Gamma ({\hat {n_{i}}}\bigtriangledown _{i}\phi )]dL=S_{\phi }\Delta A_{p}$(3-3)

${\frac {\partial (h_{P}\phi _{P})}{\partial t}}\Delta A_{P}+\sum _{f}h_{f}\left[U_{f}\phi _{f}-\Gamma _{f}(\bigtriangledown _{\perp }\phi )_{f}\right]\Delta l_{f}=S_{\phi }\Delta A_{p}$(3-4)

where

${\hat {n_{i}}}=({\hat {n_{1}}},{\hat {n_{2}}})$ = outward unit vector normal to cell face f

$V_{f}=({\hat {n_{i}}}V_{i})_{f}$ = outward cell face velocity

h_f = linearly interpolated total water depth at the cell face f

$(\bigtriangledown _{\perp }\phi )_{f}=({\hat {n_{i}}}\bigtriangledown _{i}\phi )_{f}$ is the outward normal gradient of $\phi$ at cell face f

The above equations, the Green-Gauss Theorem has been used to convert the area integral to a boundary integral. The symbol () indicates the cell face linear interpolation operator. The cell face velocity,$V_{f}$ , is calculated using a momentum interpolation method similar to that of Rhie and Chow (1983) and is described in a subsequent section.

Temporal Discretization

The general transport equation is rewritten as

$\int {\frac {\partial (h\phi )}{\partial t}}dt=\int Fdt$(3-5)

where includes all the remaining terms, such as …. For stability and efficiency, a fully implicit time-stepping scheme is used in the form

${\frac {(1+0.50\theta )h^{n+1}\phi ^{n+1}-(1+\theta )h^{n}\phi ^{n}+0.5\theta h^{n-1}\phi ^{n-1}}{\Delta t}}=F^{n+1}$(3-6)

where $\theta$ is a weighting factor between 0 and 1. For $\theta =0$, the scheme reduces to the first-order backward Euler scheme, and with $\theta =1$, the scheme reduces to the second-order backward scheme (Ferziger and Peric 1997). The superscripts indicate the time step levels, with n+1 being the current timestep and n being the previous timestep.

Cell-face interpolation operator

The general formula for estimating the cell-face interpolation operator (?) value of $\phi _{f}$ is given by

$\phi _{f}=f_{\perp }\phi _{N}+(1-f_{\perp })\phi _{P}+f_{\perp }(r_{\Box }\bigtriangledown _{\Box }\phi )_{N}+(1-f_{\perp })(r_{\Box }\bigtriangledown _{\Box }\phi )_{P}$(3-7)

where $f_{\perp }$ is a linear interpolation factor, $\bigtriangledown _{\Box }$ is the gradient operator in the direction parallel to face f, and $r_{\Box }$ is the distance from the cell center to the ghost point O parallel to the cell face f (see Figure 3.2).

Figure 3.2. Schematic showing two types of refined cells.

By definition,$\Box =2\mid {\hat {n_{1}}}\mid +1\mid {\hat {n_{2}}}\mid$ and $\perp =1\mid {\hat {n_{1}}}\mid +2\mid {\hat {n_{2}}}\mid$ . Note that for neighboring cells without any refinement, $r_{\Box }$ is equal to zero and the above equation is consistent with non-refined cell faces. The linear interpolation factor is defined as

$f_{\perp }={\frac {x_{\perp ,f}-x_{\perp ,P}}{x_{\perp N}-x_{\perp ,P}}}={\frac {\Delta x_{\perp ,P}}{\Delta x_{\perp ,P}+\Delta x_{\perp ,N}}}$(3-8)

where $x_{\perp ,f}$ is the coordinate of f perpendicular to the face and $\Delta x_{\perp }$ is the cell dimension perpendicular to the face f.

Advection Schemes

Hybrid Scheme

The hybrid scheme is a composed of a first-order upwind scheme and a second-order central difference scheme. When the Peclet number, $P_{f}$, is larger than 2, the first-order upwind scheme is used; otherwise, the second-order central difference scheme is used:

$\phi _{f}=\left\{{\begin{aligned}&(\phi _{D}+\phi _{c})/2&for\mid P_{f}\mid <2\\&\phi _{c}&for\mid P_{f}\mid >2\end{aligned}}\right.$(3-9)

where the subscripts D and C indicate the first downstream and first upstream nodes and $P_{f}=U_{f}\mid \delta _{\perp }\mid /{\bar {\Gamma _{f}}}$ is the Peclet number at the cell face in which $\mid \delta _{\perp }\mid =\mid x_{\perp ,N}-x_{\perp ,P}$ .

Exponential Scheme

The exponential scheme interpolates the face value using an exact solution to the 1-D steady advection-diffusion equation:

${\frac {\phi _{f}-\phi _{0}}{\phi _{L}-\phi _{0}}}={\frac {exp(P_{f}x_{f}/\mid \delta _{\perp }\mid )-1}{exp(P_{f})-1}}$(3-10)

where $\phi _{0}=\phi \mid _{x=0},\phi _{L}=\phi \mid _{x=L}$ . The exponential scheme has automatic upwinding and is stable, but is usually less than second order.

Hybrid Linear/Parabolic Scheme

The Hybrid Linear/Parabolic Approximation scheme of Zhu (1991) may be written as

$\phi _{f}=\left\{{\begin{aligned}&\phi _{c}+(\phi _{d}-\phi _{c}){\hat {\phi _{c}}}&for0\leq {\hat {\phi _{c}}}\leq 1\\&\phi _{c}&otherwise\end{aligned}}\right.$(3-11)

where the subscripts D, C and indicate the first downstream and first and second upstream cells, respectively. The normalized variable, ${\hat {\phi _{c}}}$ , is determined based on the formulation of Jasak et al. (1999)

${\hat {\phi _{c}}}={\frac {\phi _{C}-\phi _{U}}{\phi _{D}-\phi _{U}}}=1-{\frac {\phi _{D}-\phi _{C}}{2(\bigtriangledown _{\perp }\phi )_{C}\delta _{\perp ,C}}}$ (3-12)

where $\sigma _{\perp ,C}=x_{\perp ,D}-x_{\perp ,C}$ . The Hybrid Linear/Parabolic Approximation scheme is second order.

Cell-face gradient operator

A linearly exact second-order approximation for the normal gradient at cell face f is calculated using the auxiliary node concept of Ferziger and Peric (1997)

$(\bigtriangledown _{\perp }\phi )_{f}={\frac {\phi _{N}-\phi _{P}+(r_{\Box }\bigtriangledown _{B}ox\phi )_{N}-(r_{\Box }\bigtriangledown _{\Box }\phi )_{p}}{\mid \sigma _{\perp }\mid }}$(3-13)

where the subscripts P and N refer to two neighboring cells, $\mid \sigma _{\perp }\mid =\mid x_{\perp ,N}-x_{\perp ,P}\mid$ is the distance between cells P and N, normal to the cell face (see Figure 3.2), and $\bigtriangledown _{\Box }$ is the gradient operator in the direction par-allel to face f. Ham et al. (2002) compared the auxiliary node formulation to the fully-unstructured discretization proposed by Zwart et al. (1998) for the viscous terms and found that the auxiliary node formulation is significantly more stable.

Cell-centered gradient operator

The cell-centered gradient operator is calculated using the Green-Gauss Theorem as

$\int _{A}\bigtriangledown _{i}\phi dA=\sum _{f}{\hat {n_{i}}}\phi _{f}\Delta l_{f}$ (3-14)

Equation 3-14 above is second order and conservative for regular and nonuniform grids.