Projection Method for Solving Stokes Flow

by   Ryan Hermle, et al.

Various methods for numerically solving Stokes Flow, where a small Reynolds number is assumed to be zero, are investigated. If pressure, horizontal velocity, and vertical velocity can be decoupled into three different equations, the numerical solution can be obtained with significantly less computation cost than when compared to solving a fully coupled system. Two existing methods for numerically solving Stokes Flow are explored: One where the variables can be decoupled and one where they cannot. The existing decoupling method the limitation that the viscosity must be spatially constant. A new method is introduced where the variables are decoupled without the viscosity limitation. This has potential applications in the modeling of red blood cells as vesicles to assist in storage techniques that do not require extreme temperatures, such as those needed for cyropreservation.



There are no comments yet.


page 4

page 5

page 9

page 10

page 11

page 13

page 15

page 32


Successive finite element methods for Stokes equations

This paper will suggest a new finite element method to find a P^4-veloci...

A locally calculable P^3-pressure using a P^4-velocity for incompressible Stokes equations

This paper will suggest a new finite element method to find a P^4-veloci...

A projection method for Navier-Stokes equations with a boundary condition including the total pressure

We consider a projection method for time-dependent incompressible Navier...

Exact pressure elimination for the Crouzeix-Raviart scheme applied to the Stokes and Navier-Stokes problems

We show that, using the Crouzeix-Raviart scheme, a cheap algebraic trans...

Numerical analysis of a projection-based stabilized POD-ROM for incompressible flows

In this paper, we propose a new stabilized projection-based POD-ROM for ...

Assessment of CFD capability for prediction of the Coandă effect

The tendency of a jet to stay attached to a flat or convex surface is ca...

Spontaneous Symmetry Breaking for Extreme Vorticity and Strain in the 3D Navier-Stokes Equations

We investigate the spatio-temporal structure of the most likely configur...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.


Thank you to all of the excellent professors I have had at UW. It has been an amazing experience. A special thank you goes to Chris Vogl*, for believing in me and for giving me the opportunity to do this research with him.


Center for Applied Scientific Computing, LLNL

1 Introduction

Figure 1: Red Blood Cells in Fluid [1]

The fluid dynamics of red blood cells are of great interest to the scientific community. Biological research [2] has discovered organisms that can fill their cells with a specific sugar, increasing the viscosity of the intercellular fluids to such high levels that cell metabolism effectively stops, putting the cell in a state of hibernation. It has been hypothesized that this same method of hibernation could be applied to red blood cells. This would have a tremendous impact on the ability to store blood for critical medical procedures, such as blood transfusions, as the existing method of cryopreservation requires extremely low temperatures to preserve the blood cells.

Efficiently modeling the interaction between fluid velocity, pressure, viscosity, and body forces will catalyze the advancement of this groundbreaking storage technique. In order to begin modeling the fluid dynamics, some preliminary simplifying assumptions are made:

  1. The cell, being roughly the shape of an ellipsoid as seen in figure 2, exhibits radial symmetry, and thus its dynamics are modeled as a two-dimensional cross section.

  2. The dynamics of the inertial forces of the fluid flow, such as the beating of a heart or the stirring of the fluid, are significantly weaker than those of the viscous forces.

  3. The red blood cells are modeled, for the purposes of this project, as a vesicle, which consists only of the outer membrane of the cell.

  4. The fluid surrounding and inside of the vesicle is incompressible, sharing similar properties to water.

Figure 2: Rough Geometric Shape of a Red Blood Cell

The Reynolds Transport Theorem is used to obtain the equations for conservation of mass and conservation of momentum, which are combined to form the incompressible Navier Stokes Equation:



represents the density, p represents the pressure, and the vectors

and are defined as:


Here, and represent the horizontal velocity of the fluid, and represent the horizontal and vertical components of the force. represents the force from the pressure on the fluid, represents the viscous force, represents the inertial force, and represents external body forces, such as gravity or the force from a cell wall.

The incompressibility is imposed through the conservation of mass equation:


Because the density cannot change, . It is assumed that is also constant with respect to space, so . This yields the divergence-free condition:


which will be a very important result throughout this paper.

To properly model the behavior of interest, the full Navier Stokes Equation 1.1 is non-dimensionalized to reflect the idea that the viscous forces are much greater than the inertial forces. Variables with a tilde over them indicate characteristic quantities:


Substituting these terms into equation 1.1 and simplifying yields


Defining the Reynolds Number as


and dropping the star notation converts equation 1.6 to


Physically speaking, the Reynolds number from equation 1.7 represents the ratio of the inertial forces to viscous forces. Thus, in line with the prior assumptions, the dominant viscous forces result in a small Reynolds Number. As the Reynolds Number approaches zero, the full Navier Stokes equation from 1.1 reduces to:


Define the Laplace Operator as


If the viscosity is spatially constant (), the term becomes just . The term goes to zero because of the divergence-free condition:


Thus, when is spatially constant, equation 1.9 becomes:


2 Existing Methods

2.1 Saddle-Point Method

For the Saddle-Point Method, the velocity and pressure are simultaneously solved as the following system of equations:

Figure 3: Pressure Grid,

To perform this task, , , and are discretized and then stacked into a single vector, denoted as . The number of discretization points in a single direction, defined as M, will determine the number of rows in . Since the domain is two-dimensional, each field is represented by discretized values to span the entire domain. The grid for pressure is visualized in Figure 3. The grids for and take on the same form. Thus, has a total of rows, with rows for each of , , and :


where is to the right of , and is above .

Next, a matrix is constructed to approximate the differential operators acting on each discretized variable of , , and using the second-order finite-difference formulas:

Center-Difference Schemes:


Forward- and Backward-Difference Schemes:


The size of is . Each row represents one equation for one of the discretized variables. This creates a large system


which has the form


where G,L, and D are discretized versions of the gradient, Laplacian, and divergence operators. In order to properly index through rows and columns, the intuitive two-dimensional subindices can be mapped to a single index using the MATLAB function sub2ind, and then points can be added to the linear index to get to the next variable. A more robust way to map the indices, which is used in most of the algorithms presented, is to create an mapping matrix for each variable that stores the single index corresponding to each two-dimensional subindex for that variable.

The matrices and are constructed by looping through each of the rows and applying the correct finite differences that apply to the discretized variables that are represented in that row. Once the matrices are constructed, the system is easily solved using the MATLAB "\" operator, also known as mldivide. From there, the solved values of , , and are extracted.

2.2 Fluid in a Pipe

The Saddle-Point Method, as well as methods presented later, is used to model fluid flow in a pipe. The flow is driven by a pressure difference between the left and right edges of the pipe. It is assumed that there are no body forces, so the only force will be in the horizontal direction from the pressure difference. The vertical velocity is assumed to be zero.

A no-slip top and bottom boundary condition is applied, assuming that there is adhesion between the wall of the pipe and the fluid such that the horizontal velocity of the fluid is zero at the top and bottom of the pipe. The pressure is expected to change in the horizontal direction due to the applied gradient, but it is not expected to change in the vertical direction, so a Neumann boundary condition of is applied at the top and bottom of the pipe. Finally, is expected to change in the vertical direction due to the no-slip condition on the walls of the pipe, but it is not expected to vary horizontally, so a Neumann boundary condition of is applied to the left and right edges of the domain. A summary of these conditions is depicted in Figure 4.

Figure 4: Boundary Conditions for Fluid in a Pipe

The values for , , and are chosen arbitrarily to be , , and , respectively. These boundary conditions are then placed into the matrix from equation 2.8, using the indices of the discretized variables that reside near the boundaries. One-sided difference formulas are used for the Neumann conditions.

2.3 Checkerboarding and the Staggered Grid

Figure 5: Surface Plot of Horizontal Velocity
Figure 6: Surface Plot of Pressure, Showing the Checkerboard Pattern

Applying the Saddle-Point method to the setup specified in section 2.2 produces the surface plots for and shown in Figures 5 and 6. The surface plot for looks as expected: zero near the walls of the pipe due to the no-slip condition and a parabolic profile caused by the pressure gradient. However, the pressure graph depicts a clear problem. It seems to be very inconsistent, not smooth, and oscillating wildly. This is what is known as a checkerboard pattern, and it arises from the way the pressure has been discretized and the particular finite difference formulas used.

To explain this phenomenon, consider the following half of equation 2.8 written in component form and with zero body force:


and recall that the second-order finite difference formula used to represent is


The key detail is that the pressure on the left boundary , which helps determine the pressure gradient for the entire system, is only connected to every other pressure node due to how the finite difference formula skips over the node that is being approximated. This creates a chain connecting every other pressure node horizontally across the entire domain. This can be seen in Figure 6, as every other pressure node looks correct, smoothly transitioning from the left boundary condition to the right boundary condition, while the other nodes take on a different value because they are not connected to the boundary conditions.

In order to address this issue, a staggered grid is created as illustrated in Figure 7.

Figure 7: Staggered Grid for p, u, and v

This orientation of cells allows the use of a modified second-order finite-difference formula that uses consecutive grid values of to calculate the first derivative. To derive this formula, Taylor expansions are done about the and nodes, using the p cells just to the sides of these nodes with only a distance of . For the sake of simplicity, it is assumed that the grid spacing is created such that :

Subtracting these two equations and simplifying yields

The same logic is applied for :

In index form, these finite difference formulas can be written as follows:

The same logic is applied to calculate derivatives of at a cell, but the indices are shifted by one because of the way the grid is structured:

This algorithm will connect each cell in a contiguous chain, so the first derivative formula will no longer skip every other cell. Other operators that do not involve cells of a different type are unchanged, using the normal finite-difference formulas. This modification produces the correct pressure surface plot, as shown in Figure 8. A quiver plot, which shows the magnitude and direction of the velocity vector on a two-dimensional plane, is also shown in Figure 9.

Figure 8: Correct Pressure Field Using the Staggered Grid
Figure 9: Vector Field for Fluid in a Pipe

2.4 Decoupling Method

While the Saddle-Point method correctly solved the fluid in a pipe problem, it is computationally expensive. An alternate method utilizes a divide and conquer approach to increase efficiency.

Taking the divergence of both sides of the Stokes Equation (1.9) yields interesting results:


By rearranging the mixed partial derivatives and applying the divergence free condition , one obtains


As a result, the pressure and velocity are decoupled into the following three equations:


This very effective and simple decoupling of the pressure and velocity depends on the fact that no extra terms are generated from applying derivatives to the spatially constant viscosity term . Note, to fully model the behavior of the vesicles in the scenario of interest, the viscosity is not spatially constant, so this method cannot be used in that case.

The solution is obtained using three linear system solves of the form , where each vector now only contains rows, and is now . The matrix in each of the solves is still constructed by looping through the rows and assigning the finite differences that apply. Solving the system created by the boundary conditions from the fluid in a pipe scenario yields the smooth pressure and velocity graphs as seen in Figures 10 and 11. In this case, the values for , , and are chosen as , , and , respectively.

Figure 10: Pressure Field Generated from the Decoupling Method
Figure 11: Horizontal Velocity Field Generated from the Decoupling Method

2.5 Analytic Solution for Fluid in a Pipe

In order to verify that the algorithms are working correctly, the analytic solution for this problem is easily produced. The calculations are especially easy if the domain is defined as and . The decoupled system of equations 2.13 through 2.15 are used to derive the solution.

Since the pressure gradient is horizontal and there are no body forces acting on the fluid, an Ansatz can be made that that . Thus, becomes , which is easily solved by integrating twice:


Enforcing the boundary conditions that on the right side and on the left side yields


A similar Ansatz is made for , the horizontal component of velocity. Since varies only with y, the Ansatz is . Noting that the component of can be easily calculated using the analytic solution. A simplified equation for is found:


This equation can also be easily solved by integrating twice:


Enforcing the no-slip boundary conditions, specifically that at the top and bottom boundaries, yields the analytic solution for :


Additionally, it has already been noted that .

2.6 Convergence Testing

Figure 12: Discrete L2 Error for Pressure. The dashed line is a reference line for second-order.
Figure 13: Discrete L2 Error for Horizontal Velocity
Figure 14: Discrete L2 Error for Vertical Velocity

The numerically solved values for , , and are compared against the values of the analytic solution using the discrete L2 error. Since second-order finite-difference formulas are used in each algorithm, the discrete L2 error should be . The error calculation is defined in the following equations


where , , and represent the values of the analytical solutions calculated at the corresponding spatial locations.

Figures 12 through 14 show the discrete L2 error for , , and , calculated using the Saddle-Point method with the staggered grid and values of , , and for , , and . Since the analytic solutions for the fluid in a pipe model are a first order polynomial for , a second order polynomial for , and identically 0 for , the second-order finite difference scheme approximates the system with machine-order precision. The graphs used to present the L2 error are log scaled in both the x-axis and the y-axis, where the x-axis represents the current discretization distance and the y-axis represents the discrete L2 error for the respective variable.

The reference line on the graphs is a line with a slope of 2, created by graphing vs on the log-log scale. The reasoning behind this is that if the discrete L2 error scales as , graphing on the horizontal axis and on the vertical axis should result in a line with slope 2. The reference line is then used to verify the algorithm is converging with second-order precision, which was of course not necessary in this particular example due to the machine-level precision. Nearly identical values are obtained when using the Decoupling Method for the same scheme.

3 Projection Method for Stokes Flow

3.1 Derivation

The existing methods are established and convergence has been shown, but the problem remains that the Saddle-Point method is relatively slow and the Decoupling method cannot solve for spatially varying viscosity. A method that decouples the pressure solve from the velocity solve but does not assume spatially constant viscosity is needed. This problem motivates the use of the Projection Method, which begins with what is known as the Helmholtz-Hodge Decomposition Theorem:

Theorem 1

A vector field defined on a simply connected domain can be uniquely decomposed into a divergence-free component, , and a curl-free component, :


The full Navier Stokes Equation (1.1) is rearranged to


Taking a derivative with respect to time of the divergence-free condition gives


Thus, equation 3.2 is in the form of the Helmholtz-Hodge Decomposition (3.1). is divergence-free, as previously stated, so it is the divergence-free component . Letting , is the curl-free space. The left-hand side of the equation represents the vector field .

Define an inner product between two vector fields and as


where is the two-dimensional domain. The projection of vector field onto the curl-free is thus defined as




Using this operator, consider . Applying a two-dimensional integration by parts, is given as


Since (equation 3.3), the double integral on the right-hand side vanishes. If boundary conditions are applied such that along the boundary, then




Similar reasoning can be used to show that as well. Note that :


These concepts are well known through the work of Alexandre Chorin [3]. The algorithm, originally known as Chorin’s projection method, is used to numerically solve the incompressible Navier-Stokes Equation (1.1). It takes advantage of the fact that the equation can be written in the form of the Helmholtz-Hodge Decomposition. Its general steps are outlined as follows:

  1. Compute an intermediate velocity, , that ignores the pressure gradient term. This takes the solution off of the true divergence-free solution space.

  2. Use the intermediate velocity to solve for the pressure.

  3. Use the pressure and the intermediate velocity to project the solution back onto the divergence-free solution space, obtaining a solution for that represents its steady state after units of time have passed.

This idea is illustrated in Figure 15.

Figure 15: Illustration of the Projection Method

To solve the scenario of interest in this paper, Chorin’s projection method is modified in order to obtain a solution to the Stokes Flow equation. A number of steps must be taken to obtain from . Starting with a rearranged version of the full Navier Stokes Equation, the operator (3.6) is applied to both sides:


Applying the identities (3.9) and (3.10) reduces equation 3.12 to


Since , can be added within the projection operator and then subtracted outside. defines the amount of time between solutions of :


Both sides of the equation are then integrated in time from to :


The first integral is evaluated using a left-hand rectangular approximation, and the second integral is evaluated using a right-hand approximation. This creates a relationship between and :


Note that this is where this algorithm differs from Chorin’s projection method. Adding the additional terms to the equation allows the user to apply the small Reynolds Number and still obtain a useful equation. The Reynolds Number is now taken to be zero, so the scheme applies to 1.9:


Solving for and distributing the into the operator yields


Denote as


Then can be substituted into equation 3.19 to recover


Now the relationship between and must be established. Subtract equation 3.12 from equation 3.11:


where and are used. Take to obtain


Similar to before, is added to both the P and identity operators. The equation then becomes


Multiply by :


The interior is now replaced with :


Next, apply equation (3.21):


Taking the divergence of both sides causes to vanish due to the divergence-free condition:


The system is complete. The Projection Method is now formally defined as the following three-step algorithm:

  1. Obtain via

  2. Solve for . Note that boundary conditions are required for :

  3. Obtain via


If the viscosity does not vary spatially, then equation (3.29) becomes


3.2 Modeling a Vesicle Membrane

The circular vesicle membrane in fluid applies a force to the fluid around it, causing a jump in pressure across the membrane. This force is modeled by a modified Dirac delta function that is applied in a neighborhood of the membrane and nowhere else. Further details on this approach, known as the Continuous Surface Force approach, can be studied in papers such as [4].

The distance of the vesicle wall from the left edge of the domain is denoted as , and the distance of the vesicle wall from its center is denoted as . Thus, the center of the vesicle is located at coordinates . The signed distance from the vesicle membrane, denoted , is then


The modified Dirac Delta function is defined as


Graphically, this will create a smooth jump upward at the membrane with a height of and a width of , as shown in figure 16 for varying values of .

Figure 16: Modified Dirac Delta Function

The actual value of in the algorithm is chosen to be to avoid numerical issues. If this value is too small, the jump will be under-resolved and the discretized variables will skip right over the force. If it is too big, then the force from the membrane will be spread out too far and will not properly capture the physics. Since represents the distance from the membrane wall, points in the direction of maximum increase, which is always normal to the membrane. Thus, multiplying by will create a force vector that points outward from the membrane. Noting that is the curvature () of the membrane, Young-Laplace’s Law gives that the force should be :


The gradient of is calculated as follows:


where represents the small floating point value returned by the eps Matlab command in order to prevent division by zero when the coordinate lies exactly on the membrane. A graph of the force using , , , and is shown in Figure 17. Dirichlet conditions are used in this scenario, setting , , and to at all of the boundaries. A contour plot shows where to identify the membrane.

Figure 17: Force from the Membrane

3.3 Analytic Solution for the Vesicle Membrane Problem

If a spatially-constant viscosity is assumed, the decoupled system of equations from 2.13 to 2.15 can be used. Equations 2.14 and 2.15 can be simplified to provide the quickest route to the solution for . An ansatz is made such that is identically zero, since the membrane is circular. Under this assumption, vanishes and equations 2.14 and 2.15 become


Because the fluid is stationary, the pressure is the only non-zero variable that needs to be solved. Since the membrane force acts radially outward with the magnitude depending only on , it can be assumed that . This implies converting to polar coordinates to obtain a solution. Center the polar coordinates at the center of the vesicle, and define as the distance from that center. With this definition, the signed distance from the membrane, , becomes


since R is the distance from the center of the vesicle to the membrane. The gradient of a function in polar, assuming that , is . Thus, the gradient of in this coordinate system is trivial: . The modified Dirac Delta Function, as defined in equation 3.34, when converted to polar becomes

Thus, within the region becomes


Applying the polar gradient operator to yields . Using these identities, equation 3.37 becomes


for . Equating the vector components and integrating obtains


for . Note that


for . If , then due to the specified boundary conditions. The pressure must remain continuous at . Thus,


for . If , the pressure must remain continuous at . Thus,


The substitution is used to convert back to cartesian coordinates and produces within the region near the membrane:


Thus, the analytic solution for the circular vesicle in fluid is


3.4 Staggered Grid for Projection Method

Figure 18: Descretized Grid for the Projection Method

A different staggered grid is used for this algorithm than what was used for the Saddle-Point method, but the motivation is the same: to prevent a checkerboard pattern caused by skipping information from adjacent cells. The and cells are arranged as shown in figure 18. The distance is defined as the distance between neighboring cells of the same type.

Discretization is done at each cell, so cells using values of cells of a different type must use the average values of those cells around them. Neighboring cells of a different type are a distance of away in the horizontal direction and away in the vertical direction. To derive a finite-difference approximation to at a node, Taylor expand and about :


Subtracting these two equations and simplifying yields: