The Three-Dimensional Jump Conditions for the Stokes Equations with Discontinuous Viscosity, Singular Forces, and an Incompressible Interface ††thanks: This work was supported by the University at Buffalo SUNY and NSF Award #1253739
The three-dimensional jump conditions for the pressure and velocity fields, up to the second normal derivative, across an incompressible/inextensible interface in the Stokes regime are derived herein. The fluid viscosity is only piecewise continuous in the domain while the embedded interface exerts singular forces on the surround fluids. This gives rise to discontinuous solutions in the pressure and velocity field. These jump conditions are required to develop accurate numerical methods, such as the Immersed Interface Method, for the solutions of the Stokes equations in such situations.
Key words. Stokes equations, discontinuous viscosity, singular force, immersed interface method, jump conditions, incompressible interface, inextensible interface
AMS subject classifications. 65N06, 65N12, 76Z05, 76D07, 35R05
Biological system containing cells can be considered as a multiphase fluid problem, with the cell membranes forming a boundary between fluids of differing material properties. Unlike standard multiphase fluid problems, such as droplets or bubbles, biological multiphase systems, such as red blood cells surrounded by blood plasma, are characterized by inextensible membranes and singular forces due to both a surface-tension like contribution and a resistance to bending . Various techniques have been developed to model such systems. These techniques can be split into two categories, based on how they track the location of the interface. Lagrangian methods, such as the boundary integral method [24, 25] explicitly track the location of the interface. While these method can be highly accurate, they are difficult to extend to three-dimensions and are not well suited when the membrane undergoes large deformations. Eulerian method, on the other hand, implicitly track the location of the interface. Examples of Eulerian methods include the phase-field method [2, 13] and the level set method [18, 19]. Advantages of the Eulerian methods include easy extensibility to three dimensions and the ability to handle large membrane deformations naturally. The major disadvantage is the added computational cost of implicit membrane tracking.
Eulerian methods can be further divided into two categories depending on how they treat the singular forces arising at the embedded membrane. The first type, characterized by the immersed boundary method  or the continuum surface force method , distribute (smooth) the singular forces over a small region near the interface. This essentially turns the membrane forces into localized body-force terms, which are then included in the fluid governing equations. For most situations this smoothing results in a first-order accurate method for the fluid equations .
The second type of Eulerian method avoids the smearing of the interface by explicitly incorporating the singular forces and embedded boundary conditions in the field equations. These forces and immersed boundary conditions are included by explicitly considering the jump in the solution and derivatives of the solution across the interface. Mayo  used this idea to solve the Poisson and Biharmonic equation on irregular domains. Leveque and Li called this technique the Immersed Interface Method (IIM) and used it to solve Elliptic equations . They later extended this to Stokes flow with elastic boundaries or surface tension . Later the immersed interface method was further developed to model the Navier-Stokes equations [7, 12, 27]. Interested readers are referred to the book by Li and Ito for more information about the immersed interface method .
The application of the immersed interface method to a multiphase fluid system requires that the jumps in the velocity and pressure field be explicitly calculated. For standard multiphase systems these jump conditions for a discontinuous viscosity across the interface have been determined both in two dimensions [22, 23] and three dimensions [5, 28]. Only recently has work been done on extending the immersed interface method to multiphase flows with inextensible membranes . This recent work, though, is limited to the two-dimensional, constant-viscosity case and is not applicable to the more general discontinuous viscosity situation.
In this work the three-dimensional jump conditions for a multiphase Stokes flow system with an inextensible interface, singular forces, and a discontinuous viscosity are for the first time derived. The jump up to the second normal derivative of both the pressure and velocity fields will be developed. These jump conditions will be used to construct an Immersed Interface solver. Several analytic test cases are also developed to verify the accuracy of the jump conditions.
This paper is organized as follows. First, the basic governing equations for the three-dimensional Stokes system with an inextensible membrane are presented. Next, in Sec. LABEL:sec:IIM the immersed interface method is briefly outlined. The jump conditions for the velocity and pressure in the Stokes equations are developed in Sec. LABEL:sec:JumpConditions. The derivation of analytic test cases is shown in Sec. LABEL:sec:TestCases while results and convergence analysis follow in Sec. LABEL:sec:Results. A brief conclusion and future work follows in Sec. LABEL:sec:Conclusion.
2 Governing Equations
Let a three-dimensional bounded domain, , contain a closed and incompressible material interface . For simplicity assume that the interface is wholly contained in the domain. The region enclosed by the interface is denoted as while the region outside the interface is denoted by , Fig. LABEL:fig:CompDomain. The fluid in each domain is modeled using the Stokes equations:
with boundary conditions on the boundary of . The pressure, , body force term, , and viscosity, , are all piecewise constant in the domain, with a finite jump occurring across the interface. The velocity field, , is assumed to be continuous across the interface. It is also assumed that the velocity obey a no-slip condition on the interface. It is thus possible to state that the velocity of the interface is given by .
Application of the volume incompressibility constraint, Eq. (LABEL:eq:divfree), to the momentum equations, Eq. (LABEL:eq:momentum), leads to a pressure-Poisson equation in each domain,
with a boundary condition on the outer domain of
where is the outward facing normal to the domain .
The incompressible interface requires that the velocity on the interface be divergence-free:
where is the surface gradient and is the surface divergence operator.
The presence of the interface results in a singular force, , exerted on the surrounding fluids. This singular force is application specific. In vesicle/cell simulations this force would be composed of a variable surface tension, , and other singular forces, :
where is the outward facing normal to the interface (from into ) and is the total curvature. An example of other singular forces which might be applied are the bending forces in vesicle and red-blood cell simulations [19, 20] while the surface tension will be chosen to enforce the surface incompressibility constraint, Eq. (LABEL:eq:surface_constraint). Future work will consider singular forces of this form.
The singular force is balanced by a jump in the normal component of the stress tensor:
where the stress tensor is given by
3 Immersed Interface Method
The governing equations as described in Sec. LABEL:sec:StokesEquations result in two coupled, but discontinuous, fluid fields. In this section the scheme developed in Ref.  and used to solve the two-dimensional Stokes equations in Ref.  is reviewed. This is also the impetus for deriving the normal jump conditions across the interface.
Begin by discritizing the domain using a uniform Cartesian mesh with a grid spacing of in all directions. While this section will focus on two-dimensional domains the extension to three-dimensions is straight-forward and used in the subsequent sections. Now consider the solution to a generic Poisson problem of the form
where both and are (possibly) discontinuous across the interface. Upon discritizing Eq. (LABEL:eq:generic_poisson) it is possible to classify all grid points as either regular or irregular, Fig. LABEL:fig:IIMStencils. Regular nodes are those nodes where the interface does not cross the discritized stencil. As the solution is continuous away from the interface the discretization at regular nodes does not need to be modified. Irregular nodes are defined as those where the interface does cross the discritized stencil. As the solution may be discontinuous across the interface it is not possible to use the standard discretization. Instead the system is modified at irregular points to take into account the discontinuity.
Consider the irregular point shown in Fig. LABEL:fig:IIMStencils. Let this point exist in the domain. A standard second-order finite difference discretization results in
The issue is the values at and exist in the domain, not the domain. This fact must be taken into account in the discretization.
Let a jump across the interface be defined as
where is a point on the interface and is the outward normal vector. Assume that the jumps in the solution, , the first normal derivative, , and the second normal derivative are known on the interface for the problem given in Eq. (LABEL:eq:generic_poisson). Using a Taylor Series expansion in the normal direction about an interface point the jumps may be extended to a grid point by
where is the signed-distance from the grid point to the interface, Fig. LABEL:fig:IIMStencils. Note that by using the signed distance functions it is possible to account for extension of the jumps into either domain.
By extending the solution jumps from the interface to the grid points it can be written that . Use this expression in the Eq. (LABEL:eq:poisson_discritization) to get
As the jumps and are known they can be moved to the right-hand side as explicit corrections on the linear system,
where is the total correction needed at an irregular node at location . In this case . This method is easily extended to irregular nodes on either side of the interface and to three-dimensional systems.
It should be noted that the extension of the jumps must be calculated up to third-order accuracy to ensure that irregular nodes have a local truncation error of . Despite this higher local truncation error the overall method will retain the second-order accuracy of the underlying discretization [1, 8, 10]. If the second-normal derivative is not taken into account the local truncation error for irregular nodes will be reduced to and the overall scheme would only be first-order.
4 Derivation of Jump Conditions
The governing equations in Sec. LABEL:sec:StokesEquations result in two coupled Poisson problems. Using the Immersed Interface Method of Sec. LABEL:sec:IIM the two systems can be written over the entire domain by including the appropriate corrections:
where the corrections and are non-zero only at irregular grid points. To calculate these corrections at irregular grid points the jumps up to the second-normal derivatives need to be derived for both the pressure and the augmented velocity, . This section will outline the derivation of these jump conditions. It should be noted that a -continuous velocity indicates that . It should also be noted that the interface can be approached from either the or domains. It is thus defined that
where is a quantity, such as velocity, of interest.
Let be the outward (pointing into ) unit normal while and are two unit tangent vectors chosen such that form a Darboux Frame. Also chose and to lie along the principle directions of the interface, Fig. LABEL:fig:LocalVectors. Denote the directions in the normal and two tangent vectors as , , and , respectively. In the following sections all derivative in a particular direction are denoted as , , and . Note that a lack of bold-face indicates a direction and not a vector.
Theorems LABEL:thm:normal_jump to LABEL:thm:slapdn are non-trivial and non-obvious relations that will be used to develop the jump conditions for the pressure and velocity.
The jump in the viscosity times the normal derivative of the velocity in the normal direction is zero:
Proof. For the fluid in each domain the incompressibility conditions holds, . Thus, the jump in the incompressibility condition is zero across the interface. In terms of the local vectors this is
The velocity as you approach the interface must also be surface-divergent free, . Therefore the the jump in the surface divergence of the velocity must also be zero,
Combining Eqs. (LABEL:eq:volume_incompressible) and (LABEL:eq:surface_incompressible) results in
As the normal vector is continuous across the interface, , the jump condition, Eq. (LABEL:eq:normal_derivative), is obtained.
The following expression holds on the interface,
Proof. The incompressibility constraint, , must hold identically in both domains. Thus, the normal derivative of this constraint must also be zero, . It is thus possible to write
where the fact that , , and are held to be true.
The following expression holds on the interface,
Proof. The full Laplacian in the local frame can be written as
while the surface Laplacian can be written as 
Thus the following holds,
Using the result of Theorem LABEL:thm:normal_jump the expression in Eq. (LABEL:eq:slapdn) is shown to hold true.
The next three theorems outline the derivation of the pressure jump conditions.
The jump in the pressure is given by
Proof. The balance of forces on the interface, Eq. (LABEL:eq:stress_balance), is
Taking the normal component of this balance results in
Using Eq. (LABEL:eq:normal_derivative) results in the jump of the pressure field, Eq. (LABEL:eq:pressure_jump).
The jump in the normal derivative of the pressure may be written as
or alternatively as
where and are the principle curvatures along the and directions and recalling that is the velocity on the interface.
Proof. As Eq. (LABEL:eq:momentum) holds in each domain the jump in the system must be zero, . Dotting this by results in Eq. (LABEL:eq:dpdn1). As has been noted elsewhere  this form is not useful for numerical simulations due to the second-order partial derivatives of the velocity.
To obtain the alternative form begin by taking the surface divergence of the force balance on the interface, , or in expanded form:
Next, write the surface divergence of a vector on the interface as
When considering the derivatives of along the tangent directions it is necessary to also take derivatives of the normal vector along the same directions. This results in the following expression along the -direction
where the derivative along the -direction is similar. In the Darboux Frame it can be written that , where is the geodesic torsion along the direction. As the -direction is a principle direction then the geodesic torsion is zero, . Using this results in the following expression,
Combining this with a similar result in the -direction gives
In a similar fashion the final portion of Eq. (LABEL:eq:dpdn_der1) can be written as
Combining Eqs. (LABEL:eq:dpdn_der1), (LABEL:eq:dpdn_der2), (LABEL:eq:dpdn_der3), and (LABEL:eq:dpdn_der4) and using Theorem LABEL:thm:dudn2 and Theorem LABEL:thm:slapdn gives the expression
Next, use Eq. (LABEL:eq:split_lap) to split the Laplacian term in the base jump condition, Eq. (LABEL:eq:dpdn1), and combine with Eq. LABEL:eq:dpdn_der5 to obtain
Letting be the geodesic curvature along the direction allows for the following to hold,
Due to the continuity of tangential derivatives, holds on the entire interface. It is thus possible to state that . Therefore
When Eq. (LABEL:eq:dpdn_der_8) is dotted by the unit normal vector and noting that , see Ref. , the expression holds. A similar expression holds in the direction.
The exact form given in Eq. (LABEL:eq:dpdn2) can be obtained by writing jumps of the form as . This allows for the following to hold,
Replacing the appropriate terms in Eq. (LABEL:eq:dpdn_der_6) results in the alternative jump in the normal pressure derivative.
The jump in the second-normal derivative of the pressure is given by
where the jump in the first-normal derivative of the pressure is given in Eq. (LABEL:eq:dpdn2).
Proof. As the pressure-Poisson equation, Eq. (LABEL:eq:pressue-poisson), holds in each domain the jump is given by . Splitting the Laplacian of the pressure into the local frame results in
Solving for the jump in the second derivative and using the result of Theorem LABEL:thm:pressure_jump results in Eq. (LABEL:eq:pnn).
The final three theorems outline the derivation of the velocity jump conditions.
The jump in the augmented velocity is given by
Proof. Expanding the jump results in
due to and on the interface.
The jump in the normal derivative of the augmented velocity is given by
where is the projection operator given by .
Proof. Dot the force balance on the interface, Eq. (LABEL:eq:stress_balance), by the tangent directions to obtain
In conjunction with Eq. (LABEL:eq:normal_derivative) a linear system can be written as
where , , and are the components of the local vectors.
As the inverse of the matrix in Eq. (LABEL:eq:jun) is the transpose of the matrix the jump can be obtained as
Realizing that the term is simply the tangential projection of the force onto the interface the jump given in Eq. (LABEL:eq:dudn) is obtained.
The jump in the second normal derivative of the augmented velocity is given by
Proof. Consider the jump of the augmented momentum balance equation, Eq. (LABEL:eq:augmented_momentum). Expand the Laplacian of the velocity and the gradient of the pressure into normal and tangential components:
Due to continuity of tangential derivatives, the unit normal, and the total curvature, it can be stated that
The jump in the pressure is simply while
Finally, noticing that is the surface gradient the jump condition given in Eq. (LABEL:eq:dudn2) can be obtained.
5 Sample Test Case
To verify the accuracy of the derived jump conditions a sample analytic test case has been created. To create this test case valid inner and outer velocity fields, and , are determined for a pre-determined interface. These velocity fields are chosen to enforce the conditions and on the interface and in the domain . Next, arbitrary inner and outer pressure fields, and are chosen. Using the velocity and pressure field the inner and outer body forces, and can be calculated using Eq. (LABEL:eq:momentum). The singular force, is obtained through the force balance on the interface, Eq. (LABEL:eq:stress_balance).
The analytic body and singular forces are then used in conjunction with the known interface velocity and appropriate boundary conditions to solve the discontinuous Stokes equations in the domain using the Immersed Interface Method. The numerical pressure and velocity fields are compared to their analytic versions to determine the overall spatial accuracy of the method.
In this example the interface is a unit sphere centered at the origin. The outer viscosity is taken to be 1 while the inner viscosity is given by . The exact inner velocity and pressure for the region are given by
This results in an inner body force of the form
The exact outer velocity and pressure for the region is
which results in a body force of
A force balance on the interface results in a singular force of
Clearly on the interface holds on and can be applied in the derived jump conditions.
In the following results the domain is a taken to be a cube. The pressure and velocity fields are obtained using a split calculation where the pressure is calculated first by solving Eq. (LABEL:eq:pressue-poisson) and then obtaining the velocity field through the solution of Eq. (LABEL:eq:momentum). It should be noted that the jump conditions allow for any Stokes field solution technique, not only the one used here.
Sample pressure results for the plane on a grid are shown in Fig. (LABEL:fig:SamplePressure). The external viscosities shown are and . It is clear that the sharp discontinuity of the pressure is captured by the method and the result is not greatly affected by the viscosity ratio.
To demonstrate that the derived jump conditions, when used with a second-order Immersed Interface Method, achieve the desired accuracy the errors for both the pressure and velocity fields are now presented. In calculating the error external viscosities ranging from to have been considered. Grid sizes ranging from to have been used, corresponding to grid spacings of to .
The error for the pressure field are presented in Fig. LABEL:fig:PressureConvergence while the corresponding velocity errors are shown in Fig. LABEL:fig:VelocityConvergence. For both the pressure and velocity fields the second-order accuracy of the underlying discretization is clearly preserved despite the presence of the discontinuities in the solution field. It should be noted that the overall accuracy of the method depends on the viscosity difference between the inner and outer fluid, with matched viscosity producing the lowest error. This is to be expected as the jump conditions are heavily dependent on this viscosity difference. Discretization errors in the surface derivative quantities, such as the surface Laplacian of Eq. (LABEL:eq:dpdn2), are amplified as the viscosity difference increases. While more accurate surface derivative calculations would reduce this dependence, it will always remain and must be taken into account.
7 Concluding Remarks
In this paper the three-dimensional, normal jump conditions for the pressure and velocity fields have been derived for an embedded incompressible/inextensible interface with a singular force and in the presence of a discontinuous viscosity field. The jump conditions take into account the additional constraint of interface incompressibility given by a surface-divergent free velocity field. The jump conditions are applicable to any Stokes solver using the Immersed Interface Method. As a demonstration of the accuracy of the jump conditions a sample analytic test case has been created. Error analysis using this test case demonstrates that a second-order accurate discretization maintains the expected accuracy despite the discontinuity of the solution fields. Future work will use these jump conditions to construct generalized Immersed Interface Methods which are capable of enforcing surface incompressibility while solve the Stokes equations for multiphase Stokes flow with discontinuous viscosity.
-  Loyce Adams and Zhilin Li, The Immersed Interface/Multigrid Methods for Interface Problems, SIAM J. Sci. Comput., 24 (2002), pp. 463–479.
-  T Biben, K Kassner, and C Misbah, Phase-field approach to three-dimensional vesicle dynamics, Phys. Rev. E, 72 (2005), p. 041921.
-  MP Carmo, Differential geometry of curves and surfaces, Prentice-Hall, 1976.
-  YC Chang, TY Hou, B Merriman, and S Osher, A level set formulation of eulerian interface capturing methods for incompressible fluid flows, J. Comput. Phys., 124 (1996), pp. 449–464.
-  K Ito, Z Li, and X Wan, Pressure Jump Conditions for Stokes Equations with Discontinuous Viscosity in 2D and 3D, Methods Appl. Anal., 13 (2006), pp. 199–214.
-  M.C. Lai and H.C. Tseng, A simple implementation of the immersed interface methods for Stokes flows with singular forces, Comput. & Fluids, 37 (2008), pp. 99–106.
-  D.V. Le, B.C. Khoo, and J. Peraire, An immersed interface method for viscous incompressible flows involving rigid and flexible boundaries, J. Comput. Phys., 220 (2006), pp. 109–138.
-  RJ Leveque and ZL Li, Immersed interface methods for Stokes flow with elastic boundaries or surface tension, SIAM J. Sci. Comput., 18 (1997), pp. 709–735.
-  Randall J Leveque and Zhilin Li, The immersed interface method for elliptic equations with discontinuous coefficients and singular sources, SIAM J. Numer. Anal., 31 (1994), pp. 1019–1044.
-  Z. Li and K. Ito, Maximum Principle Preserving Schemes for Interface Problems with Discontinuous Coefficients, SIAM J. Sci. Comput., 23 (2001), pp. 339–361.
-  , The Immersed Interface Method: Numerical Solutions of PDEs Involving Interfaces and Irregular Domains, Frontiers in Applied Mathematics, Society for Industrial and Applied Mathematics, 2006.
-  Zhilin Li and Ming-Chih Lai, The Immersed Interface Method for the Navier–Stokes Equations with Singular Forces, J. Comput. Phys., 171 (2001), pp. 822–842.
-  JS Lowengrub, A Ratz, and A Voigt, Phase-field modeling of the dynamics of multicomponent vesicles: Spinodal decomposition, coarsening, budding, and fission, Phys. Rev. E, 79 (2009), p. 031926.
-  Anita Mayo, The Fast Solution of Poisson’s and the Biharmonic Equations on Irregular Regions, SIAM J. Numer. Anal., 21 (1984), pp. 285–299.
-  Charles S Peskin, Numerical analysis of blood flow in the heart, J. Comput. Phys., 25 (1977), pp. 220–252.
-  C Pozrikidis, Modeling and simulation of capsules and biological cells, CRC Press, 2003, ch. Shell theory for capsules and shells., pp. 35–102.
-  D Russell and ZJ Wang, A cartesian grid method for modeling multiple moving objects in 2D incompressible viscous flow, J. Comput. Phys., 191 (2003), pp. 177–205.
-  D. Salac and M. Miksis, A level set projection model of lipid vesicles in general flows, J. Comput. Phys., 230 (2011), pp. 8192–8215.
-  David Salac and Michael J Miksis, Reynolds number effects on lipid vesicles, J. Fluid Mech., 711 (2012), pp. 122–146.
-  U Seifert, Configurations of fluid membranes and vesicles, Adv. Phys., 46 (1997), pp. 13–137.
-  Z Tan, D.V. Le, K.M Kim, and B.C. Khoo, An Immersed Interface Method for the Simulation of Inextensible Interfaces in Viscous Fluids , Commun. Comput. Phys., 11 (2012), pp. 925–950.
-  Z. Tan, D.V. Le, Z. Li, K.M. Lim, and B.C. Khoo, An immersed interface method for solving incompressible viscous flows with piecewise constant viscosity across a moving elastic membrane, J. Comput. Phys., 227 (2008), pp. 9955–9983.
-  Z Tan, D Le, K Lim, and B Khoo, An Immersed Interface Method for the Incompressible Navier–Stokes Equations with Discontinuous Viscosity Across the Interface, SIAM J. Sci. Comput., 31 (2009), pp. 1798–1819.
-  SK Veerapaneni, D Gueyffier, D Zorin, and G Biros, A boundary integral method for simulating the dynamics of inextensible vesicles suspended in a viscous fluid in 2D, J. Comput. Phys., 228 (2009), pp. 2334–2353.
-  SK Veerapaneni, A Rahimian, G Biros, and D Zorin, A fast algorithm for simulating vesicle flows in three dimensions, J. Comput. Phys., 230 (2011), pp. 5610–5634.
-  Jian-Jun Xu and Hong-Kai Zhao, An Eulerian Formulation for Solving Partial Differential Equations Along a Moving Interface, J. Sci. Comput., 19 (2003), pp. 573–594.
-  Sheng Xu and Z. Jane Wang, An immersed interface method for simulating the interaction of a fluid with moving boundaries, J. Comput. Phys., 216 (2006), pp. 454–493.
-  Sheng Xu and Z Jane Wang, Systematic derivation of jump conditions for the immersed interface method in three-dimensional flow simulation, SIAM J. Sci. Comput., 27 (2006), pp. 1948–1980.