Discretizing the 2D Stokes-Continuity Equations

Coefficients and right-hand-side terms of the system of equations for the Stokes and continuity equations are generated by looping over basic grid cells (fig:staggered-grid) and discretizing the governing equations associated with unknowns $v_x$, $v_y$ and $P'$ using finite difference stencils where $v_x$ is the x-velocity unknown defined on the vx-staggered grid, $v_y$ is the y-velocity unknown defined on the vy-staggered grid and $P'$ is the scaled pressure unknown defined on the staggered pressure grid. The following equation shows how the scaled pressure unknown $P'$ is related to the Pressure $P$:

eq:scaled-pressure

\[P' = \frac{P}{f_{P}}\]

The scaling factor $f_{P}$ is defined by the following equation:

eq:pressure-scaling-factor

\[f_{P} = \left(\frac{ 2\eta_{vp(0,0)b} }{ \Delta x_{avg} + \Delta y_{avg} }\right)\]

where $\eta_{vp(0,0)b}$ is the visco-plastic viscosity at the upper-left node of the basic grid, $\Delta x_{avg}$ and $\Delta y_{avg}$ are the average grid spacings in the $x$- and $y$-directions, respectively, and $P$ is pressure.

x-Stokes Discretization

The finite difference representation of the x-Stokes equation eq:x-stokes is given by:

eq:x-stokes-fd-representation

\[\frac{\Delta \sigma_{xx}}{\Delta x} + \frac{\Delta \sigma_{xy}}{\Delta y} - \frac{\Delta P}{\Delta x} = -g_x \rho \text{.}\]

As described by (Gerya, 2010), a conservative finite difference formulation on the stencil shown in Fig. allows Eq. to be re-written as:

eq:x-stokes-discretized

\[\begin{split} \left[ \frac{ 2 \eta_{vp,R} \frac{(v_{xR} - v_{xC})}{\Delta x_R} -2 \eta_{vp,L} \frac{(v_{xC} - v_{xL})}{\Delta x_L} } {\Delta x_C} \right] \\ + \left[ \frac{ \eta_{vp,D} \left( \frac{(v_{xD} - v_{xC})}{\Delta y_D} + \frac{(v_{yLR} - v_{yLL})}{\Delta x_C} \right) - \eta_{vp,U} \left( \frac{(v_{xC} - v_{xU})}{\Delta y_U} + \frac{(v_{yUR} - v_{yUL})}{\Delta x_C} \right) } {\Delta y_C} \right] \\ - f_{P} \frac{(P_R' - P_L')}{\Delta x_C} = -g_x \frac{(\rho_U + \rho_D)}{2} \end{split}\]

where the subscripts $L$, $R$, $U$, $D$, $UL$, $LL$, $UR$ and $LR$ denote nodal position to the left, right, up, down, upper-left, lower-left, upper-right and lower-right, respectively, with respect to the central velocity node $v_{xC}$ located along the right right boundary of the basic grid cell where the stencil is applied. Subscripts $b$, $p$, $x$ and $y$ refer to the basic, pressure, x-velocity and y-velocity grids, respectively. The grid spacings $\Delta x_C$, $\Delta x_L$, $\Delta x_R$ denote the x-direction grid spacings centered on the central velocity node $v_{xC}$, to the left of this node and to the right of this node, respectively. The grid spacings $\Delta y_C$, $\Delta y_U$, $\Delta y_D$, denote the y-direction grid spacings centered on the central velocity node $v_{xC}$, above this node and below this node, respectively.

fig:stencil-x-stokes

Stencil used to discretize Stokes equations in the x-direction to produce the finite difference representation shown in equation Eq.. This stencil is used to generate an equation associated with the central velocity unknown $v_{xC}$ for the large system of equations with coefficients for unknowns $v_{xC}$, $v_{xL}$, $v_{xR}$, $v_{xU}$, $v_{xD}$, $v_{yUL}$, $v_{yLL}$, $v_{yUR}$, $v_{yLR}$, $P'_R$, $P'_L$ and the right-hand-side term as shown in equation Eq.. Red circles denote nodes on the x-velocity grid, blue circles denote nodes on the y-velocity grid, light blue circles denote nodes on the pressure grid, and black circles denote nodes on the basic grid. This stencil is applied to the basic grid cell colored light gray with upper-left node indices $(i,j)_b$. Staggered grid indices are shown with small font for nodes used in the stencil. Unknowns and staggered grid spacings use the following subscript notation to denote position relative to the central velocity node $v_{xC}$: C, L, R, U, D, UL, LL, UR, and LR denote the central, left, right, up, down, upper-left, lower-left, upper-right, and lower-right positions, respectively.

Rearranging Eq. leads to the following equation with coefficients for unknowns $v_{xC}$, $v_{xL}$, $v_{xR}$, $v_{xU}$, $v_{xD}$, $v_{yUL}$, $v_{yLL}$, $v_{yUR}$, $v_{yLR}$, $P'_R$, $P'_L$ and the right-hand-side term $R_x$:

eq:x-stokes-coeff

\[\begin{split} & -\left( \frac{2 \eta_{vp,L}}{\Delta x_L \Delta x_C} + \frac{2 \eta_{vp,R}}{\Delta x_R \Delta x_C} + \frac{\eta_{vp,U}}{\Delta y_U \Delta y_C} + \frac{\eta_{vp,D}}{\Delta y_D \Delta y_C} \right) v_{xC} \\ & + \frac{2 \eta_{vp,L}}{\Delta x_L \Delta x_C} v_{xL} + \frac{2 \eta_{vp,R}}{\Delta x_R \Delta x_C} v_{xR} + \frac{\eta_{vp,U}}{\Delta y_U \Delta y_C} v_{xU} + \frac{\eta_{vp,D}}{\Delta y_D \Delta y_C} v_{xD} \\ & + \frac{\eta_{vp,U}}{\Delta x_C \Delta y_C} v_{yUL} - \frac{\eta_{vp,D}}{\Delta x_C \Delta y_C} v_{yLL} -\frac{\eta_{vp,U} }{\Delta x_C \Delta y_C} v_{yUR} - \frac{\eta_{vp,D}}{\Delta x_C \Delta y_C} v_{yLR} \\ & - \frac{f_{P}}{\Delta x_C} P_R' + \frac{f_{P}}{\Delta x_C} P_L' = -g_x \frac{\rho_U + \rho_D}{2}. \end{split}\]

The large matrix coordinates $(i',j')$ of coefficients in equation Eq. are related to the global index $I_{v_x}$ from equation Eq. as follows:

eq:stokes-coeff-global-indices

\[\begin{split} v_{xC} \text{ : } & (I_{v_x} \text{, } I_{v_x}) \\ v_{xL} \text{ : } & (I_{v_x} \text{, } I_{v_x} - \Delta I_{v_{xR}}) \\ v_{xR} \text{ : } & (I_{v_x} \text{, } I_{v_x} + \Delta I_{v_{xR}}) \\ v_{xU} \text{ : } & (I_{v_x} \text{, } I_{v_x} - 3) \\ v_{xD} \text{ : } & (I_{v_x} \text{, } I_{v_x} + 3) \\ v_{yUL} \text{ : } & (I_{v_x} \text{, } I_{v_x} - 2) \\ v_{yUR} \text{ : } & (I_{v_x} \text{, } I_{v_x} - 2 + \Delta I_{v_{xR}}) \\ v_{yLL} \text{ : } & (I_{v_x} \text{, } I_{v_x} + 1) \\ v_{yLR} \text{ : } & (I_{v_x} \text{, } I_{v_x} + 1 + \Delta I_{v_{xR}}) \\ P'_L \text{ : } & (I_{v_x} \text{, } I_{v_x} + 2) \\ P'_R \text{ : } & (I_{v_x} \text{, } I_{v_x} + 2 + \Delta I_{v_{xR}}) \\ \end{split}\]

where $\Delta I_{v_{xR}}$ is the cell index shift required to get to the right velocity node and is defined as follows:

\[\Delta I_{v_{xR}} = 3(N_{y(b)} - 1)\]

The x-Stokes stencil is not applied along basic cells located along the right boundary since the $v_{xC}$ is prescribed. For these cases the matrix element $(3I_{cell} \text{, } 3I_{cell})$ is equal to 1 and the right-hand-side is equal to the prescribed $v_x$ velocity ($v_x$ = zero for free slip and no slip). These nodes are referred to as ghost unknowns by (Gerya, 2010). When ghost boundary nodes are included in the stencil formulation some unknowns become constants that are moved to the right-hand-side of the system of equations.

y-Stokes Discretization

The finite difference representation of the y-Stokes equation eq:y-stokes-stabilization is given by:

eq:y-stokes-fd-representation

\[\frac{\Delta \sigma_{yy}}{\Delta y} + \frac{\Delta \sigma_{yx}}{\Delta x} - \frac{\Delta P}{\Delta y} - g_y \Delta t \left( v_x\frac{\Delta \rho}{\Delta x} + v_y\frac{\Delta \rho}{\Delta y}\right) = -g_y \rho \text{,}\]

which can be rewritten by applying conservative finite difference formulation of (Gerya, 2010) as:

eq:y-stokes-discretized

\[\begin{split} \left[\frac{ 2\eta_{vp,D}\frac{(v_{yD} - v_{yC})}{\Delta y_D} - 2\eta_{vp,U}\frac{(v_{yC} - v_{yU})}{\Delta y_U} }{\Delta y_C}\right] \\ + \left[\frac{ \eta_{vp,R} \left( \frac{(v_{yR} - v_{yC})}{\Delta x_R} + \frac{(v_{xLR} - v_{xUR})}{\Delta y_C} \right) - \eta_{vp,L} \left( \frac{(v_{yC} - v_{yL})}{\Delta x_L} + \frac{(v_{xLL} - v_{xUL})}{\Delta y_C} \right) }{\Delta x_C} \right]\\ - f_P\frac{(P'_D - P'_U)}{\Delta y_C}\\ - g_y \Delta t \frac{1}{4} \left(v_{xUL} + v_{xUR} + v_{xLL} + v_{xLR}\right) \frac{1}{2}\left(\frac{\rho_{RR} - \rho_R}{\Delta x_R} + \frac{\rho_R - \rho_{L}}{\Delta x_L}\right)\\ - g_y \Delta t v_{yC} \frac{1}{2} \left(\frac{\rho_{LR} - \rho_R}{\Delta y_D} + \frac{\rho_R - \rho_{UR}}{\Delta y_U}\right)\\ = -g_y\frac{(\rho_{bR}+\rho_{bL})}{2} \end{split}\]

where the subscripts $L$, $R$, $RR$, $U$, $D$, $UL$, $LL$, $UR$ and $LR$ denote nodal position to the left, right, right-right, up, down, upper-left, lower-left, upper-right and lower-right, respectively, with respect to the central velocity node $v_{yC}$ located along the upper boundary of the basic grid cell where the stencil is applied. Subscripts $b$, $p$, $x$ and $y$ refer to the basic, pressure, x-velocity and y-velocity grids, respectively. The grid spacings $\Delta y_C$, $\Delta y_U$, $\Delta y_D$ denote the y-direction grid spacings centered on the central velocity node $v_{yC}$, above this node and below this node, respectively. The grid spacings $\Delta x_C$, $\Delta x_L$, $\Delta x_R$ denote the x-direction grid spacings centered on the central velocity node $v_{yC}$, to the left of this node and to the right of this node, respectively.

fig:stencil-y-stokes

Stencil used to discretize Stokes equations in the y-direction to produce the finite difference representation shown in equation Eq.. This stencil is used to generate an equation associated with the central velocity unknown $v_{yC}$ for the large system of equations with coefficients for unknowns $v_{yC}$, $v_{yL}$+, $v_{yR}$, $v_{yU}$, $v_{yD}$, $v_{xUL}$, $v_{xLL}$, $v_{xUR}$, $v_{xLR}$, $P'_U$, $P'_D$ and the right-hand-side term as shown in equation Eq.. See figure Fig. for additional details.

Rearranging Eq. leads to the following equation with coefficients for unknowns $v_{yC}$, $v_{yU}$, $v_{yD}$, $v_{yL}$, $v_{yR}$, $v_{xUL}$, $v_{xLL}$, $v_{xUR}$, $v_{xLR}$, $P'_U$, $P'_D$ and the right-hand-side term $R_y$:

eq:y-stokes-coeff

\[\begin{split} - \left( \frac{2\eta_{vp,U}}{\Delta y_C\Delta y_U} + \frac{2\eta_{vp,D}}{\Delta y_C\Delta y_D} + \frac{\eta_{vp,L}}{\Delta x_C\Delta x_L} + \frac{\eta_{vp,R}}{\Delta x_C\Delta x_R} + c_{y,stab} \right) v_{yC} \\ + \frac{2\eta_{vp,L}}{\Delta x_C\Delta x_L}v_{yL} + \frac{2\eta_{vp,R}}{\Delta x_C\Delta x_R}v_{yR} + \frac{2\eta_{vp,U}}{\Delta y_C\Delta y_U}v_{yU} + \frac{2\eta_{vp,D}}{\Delta y_C\Delta y_D}v_{yD} \\ + \left(\frac{\eta_{vp,L}}{\Delta x_C\Delta y_C}-c_{x,stab}\right)v_{xUL}\\ - \left(\frac{\eta_{vp,L}}{\Delta x_C\Delta y_C}+c_{x,stab}\right)v_{xLL}\\ - \left(\frac{\eta_{vp,R}}{\Delta x_C\Delta y_C}+c_{x,stab}\right)v_{xUR}\\ - \left(\frac{\eta_{vp,R}}{\Delta x_C\Delta y_C}+c_{x,stab}\right)v_{xLR} \\ - \frac{f_P}{\Delta y_C}P'_D + \frac{f_P}{\Delta y_C}P'_U = -g_y\frac{(\rho_{bR}+\rho_{bL})}{2} \end{split}\]

where $c_{y,stab}$ is given by:

eq:y-stokes-stab

\[c_{y,stab} = g_y \Delta t \frac{1}{2} \left(\frac{\rho_{LR} - \rho_R}{\Delta y_D} + \frac{\rho_R - \rho_{UR}}{\Delta y_U}\right) \text{,}\]

and $c_{x,stab}$ is given by:

eq:x-stokes-stab

\[c_{x,stab} = g_y \Delta t \frac{1}{8} \left(\frac{\rho_{RR} - \rho_R}{\Delta x_R} + \frac{\rho_R - \rho_{L}}{\Delta x_L}\right) \text{.}\]

The large matrix coordinates $(i',j')$ of coefficients in equation Eq. are related to the global index $I_{v_y}$ from equation eq:global-indices-of-stokes-unknowns as follows:

eq:y-stokes-matrix-indices

\[\begin{split} v_{yC} \text{ : } & (I_{v_y} \text{, } I_{v_y}) \\ v_{yU} \text{ : } & (I_{v_y} \text{, } I_{v_y} - 3) \\ v_{yD} \text{ : } & (I_{v_y} \text{, } I_{v_y} + 3) \\ v_{yL} \text{ : } & (I_{v_y} \text{, } I_{v_y} - \Delta I_{v_{xR}}) \\ v_{yR} \text{ : } & (I_{v_y} \text{, } I_{v_y} + \Delta I_{v_{xR}}) \\ v_{xUL} \text{ : } & (I_{v_y} \text{, } I_{v_y} - 1 - \Delta I_{v_{xR}}) \\ v_{xLL} \text{ : } & (I_{v_y} \text{, } I_{v_y} + 2 - \Delta I_{v_{xR}}) \\ v_{xUR} \text{ : } & (I_{v_y} \text{, } I_{v_y} - 1) \\ v_{xLR} \text{ : } & (I_{v_y} \text{, } I_{v_y} + 2) \\ P'_U \text{ : } & (I_{v_y} \text{, } I_{v_y} + 1) \\ P'_D \text{ : } & (I_{v_y} \text{, } I_{v_y} + 4) \\ \end{split}\]

Continuity Discretization

The finite difference representation of the continuity equation Eq is given by:

eq:continuity-fd-representation

\[\frac{\Delta v_x}{\Delta x} + \frac{\Delta v_y}{\Delta y} = R_c\]

which can be rewritten by applying conservative finite difference formulation of \cite{gerya2010} as:

eq:continuity-discretized

\[\begin{split} \left[ \frac{v_{xR} - v_{xL}}{\Delta x_b} + \frac{v_{yD} - v_{yU}}{\Delta y_b} \right] = 0 \end{split}\]

where the subscripts $L$, $R$, $U$, $D$ denote nodal position to the left, right, up and down, respectively, with respect to the central pressure node $P'_{pC}$ located at the center of the basic grid cell where the stencil is applied. Subscripts $b$, $p$, $x$ and $y$ refer to the basic, pressure, x-velocity and y-velocity grids, respectively. The grid spacings $\Delta x_b$ and $\Delta y_b$ denote the x and y-direction grid spacings of the basic grid centered on the central pressure node $P'_{pC}$. Rearranging Eq. leads to the following equation with coefficients for unknowns $v_{xL}$, $v_{xR}$, $v_{yU}$, $v_{yD}$ and the right-hand-side term $R_c$:

eq:continuity-coeff

\[\begin{split} \frac{f_P}{\Delta x_b}v_{xR} - \frac{f_P}{\Delta x_b}v_{xL} + \frac{f_P}{\Delta y_b}v_{yD} - \frac{f_P}{\Delta y_b}v_{yU} = 0 \end{split}\]

fig:stencil-continuity

Stencil used to discretize the continuity equation to produce the finite difference representation shown in equation Eq.. This stencil is used to generate an equation associated with the central pressure unknown $P'$ for the large system of equations with coefficients for unknowns $v_{xL}$, $v_{xR}$, $v_{yU}$, $v_{yD}$, $P'$ and the right-hand-side term as shown in equation Eq.. See figure Fig. for additional details.

The large matrix coordinates $(i',j')$ of coefficients in equation Eq. are related to the global index $I_{P}$ from equation Eq. as follows:

eq:continuity-matrix-indices

\[\begin{split} v_{xL} \text{ : } & (I_p \text{, } I_p - 2 - \Delta I_{v_{xR}})\\ v_{xR} \text{ : } & (I_p \text{, } I_p - 2)\\ v_{xU} \text{ : } & (I_p \text{, } I_p - 4)\\ v_{yD} \text{ : } & (I_p \text{, } I_p - 1) \end{split}\]

Sparse Stokes Matrix Example

Appling the stencils from x-Stokes Discretization, y-Stokes Discretization and Continuity Discretization for a small system leads to the system of equations shown in Fig..

fig:sparse-matrix-stokes

Sparse matrix for discretized Stokes-continuity equations where $N_{x,b} = N_{y,b} = 4$. For this case $\Delta I_{v_{xR}}$ is equal to 9. Unknowns for $v_x$, $v_y$ and $P'$ are denoted along the horizontal and vertical matrix axes as x, y, and p, respectively. The associated stencil is applied to each row of the matrix where the x-Stokes stencil is applied to rows labeled as x, the y-Stokes stencil is applied to rows labeled as y, and the continuity stencil is applied to rows labeled as p. Coefficients for unknowns associated with the equation generated from a given stencil are applied to the appropriate column of the matrix. Locations of unknowns relative to a main central node denoted by C are denoted as U, D, L, R, UL, UR, LL, LR which represent up, down, left, right, upper-left, upper-right, lower-left, and lower-right, respectively. Coefficients with italic font are equal to zero since the product of the prescribed unknown and associated coefficient from the x-stokes-coeff equation is moved to the right-hand-side term. Central main nodes associated with a given stencil with coefficients not equal to 0 or 1 are denoted as C. Matrix values equal one in the central node location are associated with prescribed boundary unknowns. Zeros along the diagonal represent pressure coefficients that are equal to zero since pressure does not appear in the incompressible continuity equation. Blank spaces are zero-value coefficients for non-diagonal unknowns. Note that the central pressure node $P'$ has one non-zero coefficient equal to 1 at matrix location $(2,2)$ where the pressure boundary condition is applied.