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
- 2 Existing Methods
- 3 Projection Method for Stokes Flow
- 4 Comparisons
- 5 Conclusions
- 6 Github Repositories
The fluid dynamics of red blood cells are of great interest to the scientific community. Biological research  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:
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.
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.
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.
The fluid surrounding and inside of the vesicle is incompressible, sharing similar properties to water.
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 vectorsand 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:
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:
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.
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
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.
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.
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.
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
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
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:
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 . 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:
Compute an intermediate velocity, , that ignores the pressure gradient term. This takes the solution off of the true divergence-free solution space.
Use the intermediate velocity to solve for the pressure.
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.
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:
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
Then can be substituted into equation 3.19 to recover
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:
Solve for . Note that boundary conditions are required for :
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 .
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 .
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.
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
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: