A Simple Finite Difference Method for Time-Dependent, Variable Coefficient Stokes Flow on Irregular Domains
We present a simple and efficient variational finite difference method for simulating time-dependent Stokes flow in the presence of irregular free surfaces and moving solid boundaries. The method uses an embedded boundary approach on staggered Cartesian grids, avoiding the need for expensive remeshing operations, and can be applied to flows in both two and three dimensions. It uses fully implicit backwards Euler integration to provide stability and supports spatially varying density and viscosity, while requiring the solution of just a single sparse, symmetric positive-definite linear system per time step. By expressing the problem in a variational form, challenging irregular domains are supported implicitly through the use of natural boundary conditions. In practice, the discretization requires only centred finite difference stencils and per-cell volume fractions, and is straightforward to implement. The variational form further permits generalizations to coupling other mechanics, all the while reducing to a sparse symmetric positive definite matrix. We demonstrate consistent first order convergence of velocity in and norms on a range of analytical test cases in two dimensions. Furthermore, we apply our method as part of a simple Navier-Stokes solver to illustrate that it can reproduce the characteristic jet buckling phenomenon of Newtonian liquids at moderate viscosities, in both two and three dimensions.
The equations of Stokes flow describe the motion of fluids where the nonlinear advection terms present in the Navier-Stokes equations have been eliminated or ignored. This approximation is important in many physical scenarios where inertia is essentially negligible (), however it typically leads to the steady-state Stokes equations. In this paper we are instead interested in the time-dependent or unsteady case, motivated by its use as a building block for solving the Navier-Stokes equations. That is, a number of methods apply operator splitting to the Navier-Stokes equations in order to treat the advective terms separately from the time-dependent Stokes equations.
The goal of this paper is to derive and validate a simple, efficient, and stable method for solving the time-dependent Stokes flow equations in the presence of irregular free surfaces and solid boundaries, while supporting spatially varying viscosity and density coefficients. We solve these equations on the classic staggered Cartesian grid using a finite difference approach. However, in order to support irregular domains which arise frequently in the presence of evolving liquid interfaces and moving objects, we must determine a proper discretization for non-grid-aligned (ie., embedded) boundary conditions. The free surface condition poses a particular challenge because it involves a delicate coupling between the boundary normal, the pressure, and the components of the deviatoric stress tensor. Our approach will be to take inspiration from the finite element method, and express the Stokes flow problem in a variational form relating velocity, pressure, and deviatoric stress. We will show that this yields a hybrid approach that is discretized with simple finite differences, but which implicitly enforces complex embedded boundaries through the natural boundary conditions of the problem. We believe this to be the first method that can provide a fully implicit discretization of this problem on Cartesian grids, while simultaneously yielding a symmetric positive-definite linear system and correctly handling both variable coefficients and the difficult free surface case.
While there exist finite element and finite volume methods that can be applied to the Stokes equations on irregular domains, these methods require meshes comprised of well-shaped elements that align with the physical boundaries. Frequent remeshing becomes necessary in settings with rapidly evolving boundaries, and as others have noted, such remeshing is a challenging and expensive problem in its own right . Unstructured meshes often also incur a performance penalty as compared to regular grids due to the overhead costs of manipulating more complex data structures. We therefore pursue the simpler and more efficient embedded finite difference approach, allowing boundaries to cut arbitrarily through an underlying regular Cartesian grid.
Viscous free surface flows have long posed challenges for finite difference methods. The core problem is ensuring that the region external to the liquid applies zero force on the surface itself through the application of appropriate boundary conditions. The seminal marker-and-cell (MAC) paper of Harlow and Welch  discussed these difficulties, and dispensed with the free surface viscous stress conditions for simplicity. However, they noted that at lower Reynolds numbers the more accurate condition becomes necessary. Hirt and Shannon proposed a simple correction to the normal stress , which Nichols and Hirt subsequently improved to address the tangential stress . Because an explicit discretization of viscous terms suffers from a stringent time step restriction of that becomes problematic at high viscosities, Pracht incorporated the preceding ideas into a fully implicit integration scheme to improve stability . All of these schemes assume that the surface normal is either aligned with coordinate axes, or at a angle, in order to explicitly design conditions for these special cases. This effectively rasterizes or voxelizes the domain into a “stairstep” approximation that may not realistically reflect the physical geometry. This can limit the generality and accuracy of the approach, and the case-by-case analysis somewhat complicates the implementation. Nevertheless, Tomé et al. developed a three-dimensional extension that convincingly simulates a range of free surface phenomena , while using explicit time integration and a pressure projection method that splits the Stokes equations into separate pressure and viscous components (see  for an exploration of the issues raised by this type of splitting). Oishi et al. later implemented a more stable implicit time integration scheme  within the same framework to enable larger time steps and improve efficiency. Aside from our own work, this is the only other implicit finite difference method that has addressed the phenomenon of viscous jet buckling with proper free surface conditions. However, it requires the same case-by-case analysis as its predecessors, as well as the use of the less robust bi-conjugate gradient method to solve a large, sparse, non-symmetric linear system. A few authors have also advocated a semi-implicit treatment of spatially varying viscosity, in which terms that couple different velocity components are treated explicitly, and non-coupled terms are treated implicitly .
In computer animation, staggered grid schemes have also been used extensively for free surface flows, again using a pressure projection approach . Recently, Batty and Bridson showed that improper viscous free surface conditions in these methods tend to destroy angular momentum . To address this, they proposed a fully implicit scheme for variable viscosity that yields symmetric positive-definite linear systems, and presented animations of rotating and buckling viscous liquids. While that method shares much with the present work, it differs in several key respects. First, it uses a pressure projection approach to split up the Stokes equations which entails weaker coupling between pressure and viscous terms. Secondly, it enforces a simplified free surface boundary condition that incorrectly neglects the interaction between pressure and viscous stresses at the surface. Thirdly, the present approach yields a linear system in terms of stresses, rather than velocities. Finally, the current work provides numerical experiments indicating that our scheme is indeed convergent.
There exists a broad family of methods that seek to accurately enforce boundary conditions of various kinds on irregular domains, while solving equations on Cartesian finite difference grids. We will collectively refer to this family, which includes the current work, as embedded boundary methods (although other authors have sometimes used the term immersed boundary methods ). As noted above, the motivation for embedded approaches is to avoid the substantial computational cost of generating and manipulating high quality, fully unstructured, boundary-conforming meshes as required by finite element or finite volume methods. Common examples of embedded boundary methods include the traditional immersed boundary methods (IBM) , ghost fluid methods (GFM) , immersed interface methods (IIM) , matched interface and boundary methods (MIB) , and cut-cell methods . These have been quite effective for a number of problems and some can achieve higher order accuracy. Nevertheless, we are not aware of any symmetric fully implicit embedded boundary method with a comparably simple implementation that can address the variable coefficient unsteady Stokes equations with free surfaces. As an example, recent work that primarily emphasizes the use of the ghost-fluid approach for simulating turbulent atomization , Desjardins et al. revert to a diffuse interface method for viscous terms, stating that this choice was motivated by the complexity of viscous GFM discretizations  and the difficulty of achieving an implicit formulation.
Among existing embedded boundary methods, the most closely related to our approach are staggered grid methods that use standard centred-difference stencils, scaled by simple diagonal weighting matrices to support irregular boundaries. These methods yield symmetric positive-definite linear systems, but have primarily been applied to Poisson and diffusion problems to date. First, a simple ghost fluid method has been used to enforce Dirichlet boundary conditions  for pressure projection methods and for implicit integration of spatially constant viscosity with solid boundaries. Secondly, a finite volume-like technique has been proposed for handling Neumann (and Robin) boundary conditions in the context of pressure projection methods for solid-fluid interaction . Work by Batty et al. on solid-fluid coupling and viscous flows  also falls into this category, although they are derived from variational principles rather than ghost fluid or finite volume concepts. The current work directly extends and unifies these last two methods, while providing numerical experiments that validate the approach.
Robinson-Mosher et al. have developed a related symmetric positive-definite formulation of time-dependent Stokes flow, focusing on solid-fluid coupling . However, their viscosity discretization is limited to voxelized (axis-aligned) solid geometry and constant viscosity, and does not address the free surface boundary condition we consider. Moreover, our derivation identifies their algebraic transformation as a change of variables from a velocity-pressure formulation to a pressure-stress formulation, lending physical insight to the method, and recovers a more fundamental variational form of the mechanics that offers improved specialized solvers.
There are of course a host of methods for treating the more common symmetric indefinite form of the Stokes equations for steady and unsteady problems . Among these are several methods carefully designed to be optimal for the Stokes problem, in the sense that their convergence rates are independent of the size of the discrete simulation mesh, with “appropriate” choices of preconditioners or smoothers . Nonetheless, such methods face additional challenges in addressing boundary conditions and variable coefficients, particularly for the case of irregular embedded boundaries. For example, Wan et al. applied substantial modifications to a geometric multigrid method to efficiently handle embedded interfaces for the simpler Poisson equation . Similarly, many fast solvers for the Stokes problem rely on a simplification (discussed in Section 3.1) that decouples the viscous contributions into three independent Poisson-type problems, enabling the use of existing fast Poisson solvers as sub-components. This simplification does not readily apply in the free surface or variable viscosity settings. We therefore prefer to pursue a symmetric positive-definite formulation for which effective black-box iterative solvers and preconditioners are more readily available. The development of special-purpose solvers and preconditioners for the positive-definite formulation is a potentially exciting research direction, which we defer to future work.
3Time-Dependent Stokes Flow
The governing equations for time-dependent Stokes flow are:
where is the fluid density, is the fluid velocity, is the pressure, is the symmetric deviatoric stress tensor, is the dynamic viscosity coefficient, and the subscript indicates a time derivative. We allow both and to vary spatially. At solid boundaries, the no-slip condition applies, which dictates that the fluid velocities match those of the solid boundary, :
At a free surface, the fluid is subject to the constraint that the traction applied at the surface is zero:
In the above, is taken to be the outward normal of the free surface and represents the identity matrix.
3.1A Note on the Simplified Stokes Equations
In situations where the viscosity is spatially constant, some manipulations are often carried out to reduce the contribution of viscosity to a simpler form. Specifically:
where a simple vector calculus identity and the divergence free property of have been used to eliminate the second term on the right hand side. This reduces the contribution of viscosity to a simple Laplace operator applied to each component of velocity independently. However, the natural boundary conditions of this modified problem are quite different from those of the original problem, as discussed by Limache et al. ; if care is not taken this can lead to non-physical solutions, particularly for free surface flows. We therefore prefer the fully general form of the viscous terms, even for constant viscosity .
4Variational Formulations for Time-Dependent Stokes Flow
A backwards Euler time discretization of the governing equations yields:
In these expressions is the final velocity field, is the input (possibly divergent) velocity field, and is the size of the time step.
An equivalent variational formulation is the following:
The domain of integration is the liquid (non-air) region, . Calculus of variations can be used to show that the optimality conditions for this problem yield precisely the equations of the original PDE problem above, with the zero traction free surface boundary condition (Equation 3) enforced as a natural boundary condition.
Similarly, the following variational formulation enforces the Stokes equations inside the domain, but with static solid boundary conditions ().
Here the integration is performed over the fluid (non-solid) region, . Note that these two variational formulations differ essentially by an integration by parts operation on the middle terms.
With a particular variational formulation in hand, we proceed to directly discretize the required integrals in a manner similar to that proposed by Batty et al. . This results in a discrete optimization problem in which the boundary conditions are enforced naturally and implicitly. This is in contrast with the more common approach of first discretizing the PDE form itself, which necessitates the explicit handling of the potentially difficult boundary conditions outlined earlier.
We discretize the derivatives using centred finite differences on the classic staggered (MAC) grid, with pressures at cell centres, and velocity components on cell faces . To support viscosity we must also place the components of the deviatoric stress tensor on our grid. The most natural way to do this is to locate diagonal components of the stress tensor () at cell centres, and off-diagonal components () on cell edges (nodes in two dimensions). This arrangement is illustrated in 2D and 3D in Figures Figure 1 and Figure 2, respectively. Straightforward centred differencing can then be used to compute the required derivatives in the correct locations. This approach was first proposed by Darwish et al.  for two dimensions, and later extended to three dimensions by Mompean and Deville . While it has traditionally been applied to non-Newtonian fluids in which the constitutive equations are more complex, we find its simplicity and elegance useful for the purely Newtonian flows we consider.
Note that because the deviatoric stress does not measure compression, we are assured that which in 2D implies that . We can therefore simplify the necessary computations by solving just for rather than both quantities; in the final linear system, this will yield equations of the form , thereby retaining symmetry. Similar transformations apply in 3D to eliminate , yielding and . Naturally, the stress tensor will also be symmetric, so that , , and , which further reduces the number of variables to be computed.
We approximate the integral comprising the variational forms by scaling each term by the fractional volume of material in a cell-sized control volume surrounding the appropriate sample point, and sum over all cells. For example, terms that lie on faces are scaled by the volume of fluid in a square (or cubic) control volume surrounding the face centre. The necessary two-dimensional control volumes are illustrated in Figure 3. (While we could apply a higher order approximation of the integrals, it is this simple piecewise constant approximation that ensures we retain the same stencils as a standard grid-aligned finite difference discretization.) For free surface boundary conditions, we need to estimate volume fractions interior to the liquid (ie. not air), indicated by weights . Later, we will also need the corresponding air fractions, (where the identity). Likewise, for solid boundary conditions we estimate the volume fractions of a cell that is inside the fluid (ie. not solid), and its complementary solid fraction, . These volume fractions are illustrated in Figure 4, and can, for example, be estimated from a level set representation using the method of Min and Gibou . We will see that this simple piecewise constant approximation of the integrals ensures that we retain simple centred finite difference stencils, while volume fraction weighting handles boundary conditions.
In equation (Equation 6), we are integrating over the liquid region , so we must use volume fractions with subscript . The first term of the equation consists of velocity data that lies on faces, so we estimate the integral using fractional volumes associated to each velocity face sample. We indicate this using superscripts on the weights that indicate the type of control volume being considered; in this case the weight used is . Similarly, the second term consists of pressures and divergences that are conceptually located at cell centres, so we use cell-centred weights, . Finally, the last two terms are associated with stresses, which lie on cell edges and centres, so we use the fractions .
Applying this framework to (Equation 6) leads to the following discrete Stokes problem with free surface boundaries:
In this expression is a diagonal matrix of densities per velocity sample and is a diagonal matrix of viscosity coefficients per stress sample. The derivatives are discretized with staggered grid finite difference operators: is the usual discrete gradient, and is the discrete deformation rate applied to velocity (ie. ). Note that the negative transposes of and are the corresponding discrete vector and tensor divergence operators, respectively. The various terms are diagonal matrices consisting of the volume fractions described above.
Since this discrete optimization problem is a quadratic, the optimality conditions yield the following symmetric indefinite linear system:
The solid boundary Stokes problem (Equation 7) can be discretized in much the same manner, except that integrals are computed over the fluid (non-solid) region using volume fractions . The discrete form is:
The linear system for the solid wall Stokes problem is:
Because these two systems have nearly identical forms, and differ only by the weighting matrices , we can straightforwardly combine them to handle both free surfaces and solid boundaries in the same problem:
By combining the two formulations only at the discrete level, we are able to exploit the natural boundary conditions to handle both boundary types together, despite the fact that in each of the two continuous variational formulations only one of the two boundaries can be considered natural.
For this Stokes system, eliminating stress by applying the Schur complement to the centre block yields the symmetric indefinite velocity-pressure system most commonly associated with the Stokes problem. However, as outlined earlier we prefer to solve symmetric positive-definite systems, as they are generally more amenable to common black box solvers such as preconditioned conjugate gradient methods, domain decomposition, etc. Conveniently, the upper-left block of the full system is a diagonal matrix which can be trivially inverted. We therefore perform a Schur complement on this block to eliminate velocity and arrive at a sparse symmetric positive-definite system for pressure and stress. This is exactly analogous to classic pressure projection methods for incompressibility, in which eliminating velocity yields a symmetric positive-definite Poisson equation for pressure alone. Of course, it should be emphasized that this technique can only be applied in the time-dependent Stokes case; otherwise the upper-left block is simply zero making a Schur complement impossible.
The symmetric positive-definite form is:
where the blocks of the matrix are
Our results for the Stokes problem consistently indicate first order convergence for all variables in , and first order convergence for velocity in . Stress fails to converge in due to noisy errors along the boundaries. Given that stress is computed as the gradient of velocity, it is not entirely surprising that it loses one order of accuracy in . Fortunately, for most practical applications the behaviour of velocities is of greater importance.
6Non-Homogeneous Boundary Conditions
The preceding formulation for the Stokes problem exploits natural homogeneous boundary conditions to simplify handling of irregular domains on Cartesian grids. However, many practical situations will call for non-homogeneous boundary conditions where the boundary values are non-zero. For example, moving solid boundaries will require non-zero boundary velocities to be enforced, as will inflow and outflow boundaries. Similarly, non-zero pressure boundary conditions have been used to support surface tension (eg. ). To incorporate non-homogeneous boundary values into our framework in a consistent manner, we introduce additional terms to account for the work done by the boundary itself.
6.1Prescribed Traction Boundaries
To add a prescribed traction boundary, we need to account for the work done by the traction at the surface. We do this by adding the following boundary term to (Equation 6):
where and are the prescribed pressure and deviatoric stress, respectively, and in this case is the outward normal with respect to the air (non-liquid) region, . (Of course, only the resulting normal component, ie. surface traction, will be enforced.) This surface integral can be converted to a volume integral through integration by parts:
This ensures that all the relevant quantities are at grid locations that are consistent with the functionals considered previously. This also avoids the need to discretize surface integrals which would add further complications and potentially sacrifice the convenience of centered difference stencils.
Using terms to indicate the air fraction of a particular control volume, the discretized form is:
The new terms modify the right hand side of the linear system (Equation 8), to become:
Although the air region may extend far from the actual liquid surface, we only need to apply modifications to the right hand side for rows of the system in which the matrix has non-zero entries, indicating that there is liquid present. This means that in practice only the interface between the air and liquid regions plays a role.
6.2Prescribed Velocity Boundaries
In the common case of moving solid boundaries with prescribed velocities , we account for the work done by the solid on the fluid by adding the following term to (Equation 7):
where is the outward normal to the solid region, . In volume integral form we have:
Labelling solid fractions and discretizing, we arrive at the following term:
This results in a modification to the right hand side of the linear system (Equation 9), to become:
These right hand side modifications are also only applied to rows in which the matrix has valid non-zero entries.
7Generalization and Further Factorization
At this point, we highlight the abstract form of the problem studied above, which helps see the discrete nature of the transformation we use and illustrates how to generalize our technique to other problems such as two-way fluid-solid interaction. This is a special case of a weighted linear least-squares problem with block diagonal regularization subject to linear equality constraints, which we state using bold elements to avoid confusion with prior notation:
For the Stokes discretization above, the velocity vector here corresponds to the vector of all fluid velocities and , the block diagonal regularization matrix which we term the mass matrix for the abstract problem, corresponds to the (simply diagonal) product . The first term in the objective is thus a kinetic energy norm of the difference between the old (or predicted) velocity and the new velocity which we are finding. The constraint matrix and constraint vector correspond to the negative weighted discrete divergence operator and the zero vector (for the divergence of the Stokes velocity in the interior) modified as necessary by non-homogeneous boundary conditions. Finally, in the weighted least-squares term the matrix corresponds to (i.e. the weighted measure of deformation rate), the vector is zero in the interior of the Stokes problem but modified by non-homogenous solid boundaries as needed, and the weighting matrix corresponds to (i.e. the weighted, time-scaled viscous coefficients).
Rearranged into this format, our Stokes discretization looks like:
The solution is a balance between staying close to the old velocity (weighted by density) and minimizing the deformation rate (weighted by viscosity and time step) with cell fractions accounting for irregular geometry, subject to incompressibility. This makes evident the connection to Helmholtz’s minimum dissipation theorem for Stokes flow .
As before, we can write down the KKT optimality conditions, introducing a Lagrange multiplier for the linear constraint:
Up to a scale factor, the new variable corresponds to pressure in the Stokes problem.
This is of course a symmetric indefinite matrix, but our variable change to pressure and viscous stress generalizes here as well to arrive at a sparse SPD matrix. First introduce the weighted least-squares residual
Eliminating velocity with the first equation leaves us with the SPD Schur complement problem for the residual and the Lagrange multiplier:
A key advantage to this system is that, under the assumption that the mass matrix and the weight matrix are block diagonal or more specifically have sparse inverses, and that and don’t have overly dense columns, this is a sparse matrix itself. (If isn’t of this form, using can restore sparsity.)
Being sparse and SPD, many black-box linear solvers are available which may not apply to the original indefinite system. In our studies, we used PCG with a parallelized algebraic multiplicative Schwartz overlapping domain-decomposition preconditioner, for example.
We can go further, however, reducing the requirement on and . In particular, the matrix in equation (Equation 15) can be factored conveniently as:
These factors are obviously as sparse as the original components. Note also that the right-hand side of equation (Equation 15) can be written as
In the case of Stokes with homogeneous boundary conditions, , and it is then clear that equation (Equation 15) gives the normal equations for a sparse, weighted, and unconstrained least-squares problem:
Specialized solvers such as LSQR  may then offer a significant advantage.
The primary advantage of this generalization, and factorization, lies however in the ease of extending the formulation to more general dynamics. For example, Batty et al. coupled rigid body dynamics with inviscid incompressible flow using a pressure-only subset of equation () expressed as an unconstrained least-squares problem. However, for rigid bodies overlapping many fluid grid cells the system suffered a large dense block corresponding to those cells, leading Robinson-Mosher et al. to pursue an indefinite form . In the factored form, however, the rigid body merely corresponds to six fairly dense rows in and , which can be handled much more efficiently in both storage and multiplication with the factored form. Elastic and/or constrained solid dynamics can be phrased in the same minimization form as equation (Equation 12), e.g. English and Bridson’s isometrically deforming membranes , further opening a clear route to efficiently solvable general solid-fluid coupling: the additional degrees of freedom of the solid are appended to , the solid mass matrix is appended to , additional solid constraints are appended to and , and elastic potential energy terms are expressed as a possibly non-linear least-squares terms appended to and (as is quite natural for hyperelastic finite element models, and points out the utility of a Hessian-free Gauss-Newton iteration for solving them). Finally, we note that while Robinson-Mosher et al. discuss some of the above manipulations , they do so for a less general linear algebra problem, and do not discuss the optimization form for either the initial constrained problem nor the final unconstrained least-squares problem.
8Null Space Elimination
An issue that arises with our approach is the presence of null spaces due to overlapping volume weights assigned to different terms of the discrete variational problem. For example, consider the divergence operator for a cell near a free surface boundary; there are volume weights associated with both velocities (face centres) and pressures (cell centres). Cases frequently arise in the discretized system in which a pressure with a non-zero volume weight enforces a divergence constraint on at least one velocity face with zero associated volume. This velocity sample will appear in no other equations because of its zero volume weight, and therefore it may take on an arbitrary value as long as it satisfies the constraint. An example of such a null space scenario is shown in Figure 5.
In the final result, only samples with a positive associated volume weight are considered, so the physical solution is not adversely affected. However, such large spurious null spaces can pose problems when applying standard solvers for sparse linear systems; we would therefore like to eliminate them.
To do so, we identify each variable that enforces a relationship on a sample with zero volume weights. For example, in the free surface pressure case a non-zero weighted pressure is tagged as invalid if it enforces the divergence-free condition on one or more velocity faces with zero weights. Likewise, in the solid boundary case if a non-zero weighted velocity sample borders a zero-weighted pressure sample, that velocity sample is tagged as invalid. This process can likewise be extended to the viscous terms. Once these invalid variables have been identified, they can be straightforwardly eliminated from the linear system, and replaced with the value of the boundary condition, resulting in a modification to the right-hand-side. This reduced set of equations retains symmetry and has the same solution, but no longer suffers from large null spaces.
There can be one additional null space when the fluid domain is completely enclosed by solid walls. Because only the gradient of pressure affects the final velocities, pressure solutions that differ by a constant are effectively equivalent. This rank 1 null space doesn’t pose substantial problems, so we have not eliminated it; if necessary, it could be removed by arbitrarily fixing one pressure value in the domain.
We have verified that the method computes the exact solution for linear problems even for irregular domains. This includes the case of hydrostatic fluid with solid (and possibly free surface) boundaries, as well as rigid translations and rotations of liquid bodies with free surface boundaries. We provide a range of examples below to illustrate the convergence orders achieved by our methods in more difficult scenarios. All of the examples make use of curved boundaries which do not align with the underlying Cartesian grid, and we consider a single time step of the time-dependent problem in question. The examples test our methods in the presence of free surfaces and both static and moving boundaries. We compute the and errors, where our discrete norm is computed as for a uniform grid spacing in spatial dimensions.
For each case, we transformed the linear system to the sparse symmetric positive-definite form as described above, and solved it with the conjugate gradient method preconditioned with overlapping multiplicative Schwarz domain decomposition. The preconditioner was determined purely algebraically, from a simple graph partition of the sparse matrix; we expect the positive-definite form would allow the use of many other black-box solvers as well. Though we do not consider this iterative solver a core contribution of the paper, Table ? presents some representative iteration counts to illustrate its scaling.
|PCG Iterations||PCG Iterations|
|Grid||Free Surface (Section 9.1)||Solid Wall (Section 9.2)|
9.1Stokes Flow with Free Surface Boundaries (2D)
Our free surface Stokes test case is a fluid disk of radius centred at the origin, with density and viscosity , computed over a timestep . For simplicity of presentation, we describe the final velocity field in terms of a streamfunction, , where the velocity field can be derived as . This also guarantees that the velocity field is divergence free. The streamfunction is:
This is a non-trivial velocity field designed to fulfill the free surface zero traction condition at , smoothly blended into a zero velocity at the origin (. The zero traction condition (Equation 3) enforces a relationship between the surface pressure and the viscous stress resulting from this velocity field. To satisfy this condition, we use the following expression for pressure:
The pressure in this expression will be non-zero at the interface; any method to solve this problem will need to correctly handle the coupling between pressure and viscous stresses. From this information, the expressions for the input velocity and final stresses can be derived using equations (Equation 5)-( ?). We used a computer algebra system for this purpose. The convergence results are shown in Table ? and Figure ?.
9.2Stokes Flow with Solid Wall Boundaries (2D)
Our Stokes solid boundary test case is an annulus centred at the origin with inner radius , outer radius , density , and viscosity , computed over a timestep . Inner and outer boundaries are static solids. We will again use a streamfunction to dictate our velocity field and ensure it is divergence free:
For pressure, we use:
Equations (Equation 5)-( ?) can be used to derive the input velocities and final stresses. The convergence results are shown in Table ? and Figure ?.
9.3Stokes Flow with Both Solid Wall and Free Surface Boundaries (2D)
To test a scenario featuring both solid and free surface boundaries, we solve for fluid motion in an annulus centred at the origin with inner radius , outer radius , density , and viscosity , over a timestep . The outer boundary is a free surface, and the inner boundary is a static solid. We will again use a streamfunction to dictate our velocity field and ensure it is divergence free:
For pressure, we use:
Equations (Equation 5)-( ?) can be used to derive the input velocities and final stresses. The convergence results are shown in Table ? and Figure ?.
9.4Stokes Flow with Prescribed Velocity Solid Boundaries (2D)
To test solid boundaries with prescribed (non-zero) boundary velocities, we solve for fluid motion in an annulus centred at the origin, with inner radius , outer radius , density , and viscosity , over a timestep . The outer boundary is static, while the inner boundary rotates rigidly with a clockwise angular velocity . For the final velocity we use the streamfunction :
For pressure, we use:
As in the preceding examples, stresses and input velocities can be computed from equations (Equation 5)-( ?). Convergence results are shown in Table ? and Figure ?.
9.5Three Dimensional Stokes Flow
To test our method in three dimensions where analytical solutions are substantially more difficult to derive, we created a numerical solution at resolution , and tested convergence towards this solution. The test case consists of a sphere of liquid centred at the origin with a free surface at , containing a nested static solid sphere of radius . Density was set to and viscosity to . We compute a time step of length starting from an input velocity field:
Convergence results are shown in Table ?. Interestingly, this numerical experiment suggests first order convergence in even for stress variables, an improvement over the results in 2D. We cannot yet explain why the move to 3D would bring greater accuracy, but several similar tests showed the same pattern.
Table ? summarizes the approximate orders of convergence suggested by our experiments in two dimensions. As noted earlier, we achieve essentially first order convergence, with the exception of pressure and stress in . We now proceed to make some general comments.
A rigourous convergence theory to describe the method is beyond the scope of the current work. However, in the absence of boundaries, our methods are equivalent to a straightforward discretization of the original PDE form with second order centred finite differences on a staggered grid, and can therefore expect to achieve uniform second order convergence. Near boundaries this clearly does not hold, and the effective quadrature is likely only first order due to the use of piecewise constant approximations of integral terms. This is in line with our results in for velocity. Solution gradients can generally be expected to be one order less accurate than the solution itself, and this is also evident in the errors in the results for fluid pressure and stress. Interestingly, Chen et al. recently showed that the immersed boundary method for Stokes flow likewise exhibits errors in pressure in .
The virtual node method of Bedrossian et al.  for the Poisson equation uses a variational formulation similar to ours, but makes use of piecewise bilinear Cartesian elements near the boundary to estimate the relevant integrals, at the cost of denser stencils for boundary cells. Their results indicate second order convergence which is consistent with the fact that our use of piecewise constant estimates yields first order convergence. This also suggests that applying bilinear elements near boundaries may be effective in raising the convergence order of our method for Stokes flow, while maintaining the benefits of sparsity and positive-definiteness.
Ng et al.  pointed out that in the final discretized form, their method for the Poisson equation with Neumann solid boundaries is identical to that of Batty et al.  if the face volume weights suggested by the variational perspective are replaced with face area (finite volume) weights. The latter choice leads to an increase in accuracy for pressure from first to second order and velocity from zeroth to first order. The related ghost fluid method for the Poisson equation with Dirichlet boundaries  likewise exhibits second order in pressure and first order in velocity, with a different choice of weights. While this hints that alternative diagonal weighting matrices might raise the order of accuracy of the current method, we have found that directly introducing ghost fluid or finite volume weights in this setting breaks the symmetry of the linear system. Nonetheless, a deeper exploration of the connections between our method and ghost fluid/immersed interface methods, finite volume methods, and finite element methods might provide the key to an improved weighting scheme.
10Application to Viscous Jet Buckling
One particularly fascinating phenomenon exhibited by highly viscous liquids is jet buckling. When a falling liquid column of sufficient viscosity impacts a solid surface, it will fold or coil over on itself rather than spreading out smoothly. Relatively few researchers have looked at simulating Newtonian viscous buckling, despite its prevalence in many common liquids such as honey. To the best of our knowledge, the GENSMAC code of Tomé, McKee and co-authors is the only prior finite difference scheme to do so in three dimensions . However, as noted earlier this approach requires a case-by-case analysis of discrete surface orientations and its implicit formulation entails solving a large non-symmetric linear system. Jet buckling has also been addressed in a finite element setting  and with an SPH approach .
With this problem in mind, we incorporate our Stokes solver into a simple two-stage fractional step Navier-Stokes routine (similar to that presented in section 2 of the paper by Ng et al. ). First, starting with a velocity field at time , we compute advection and body forces to produce an intermediate velocity field :
We account for advection terms with a first order semi-Lagrangian scheme using bilinear interpolation of velocities. We then simply add any external body forces (gravity in our examples). From this intermediate velocity, we then simultaneously incorporate viscous forces and project the velocity field to be divergence free using our Stokes solver, to arrive at time :
with the appropriate free surface and solid boundary conditions applied. Tracking of the liquid surface position is performed using a basic semi-Lagrangian level set method (eg. ).
10.1Two Dimensional Jet Buckling
Figure ? presents the results of a two-dimensional simulation of planar viscous jet buckling. The simulation domain is a circle of radius centred at . A horizontal ceiling is placed at , featuring a liquid jet inflow centered at with a fixed vertical velocity of and a width of . This configuration yields a drop height of . The density of the liquid is and the dynamic viscosity of the liquid is . Gravity is set at . The simulation grid used a resolution of cells.
Following Tomé and McKee , this yields a Reynolds number of and an aspect ratio for the liquid jet of . This falls within the regime in which planar buckling is expected to occur () according to Cruikshank and Munson; as shown our method reproduces the buckling phenomenon.
10.2Three Dimensional Jet Buckling
Figures ? and ? present the results of a three-dimensional simulation of cylindrical viscous jet buckling (ie. coiling). The simulation domain is a sphere of radius centred at . A circular inlet is centred at , with a fixed vertical velocity of and a diameter . This configuration yields a drop height of . The density of the liquid is and the dynamic viscosity of the liquid is . Gravity is set at . The simulation grid used a resolution of cells.
This yields a Reynolds number of and an aspect ratio for the liquid jet of . This falls within the guidelines for when axisymmetric buckling typically arises () according to Cruikshank and Munson, and the result does indeed exhibit substantial buckling.
11Conclusions and Future Work
We have shown that a Cartesian grid finite difference method derived from a variational principle can correctly capture difficult irregular boundary conditions in Stokes flow problems, while providing stability for large time steps and yielding a sparse, symmetric positive-definite linear system. To do so, we have unified and extended recent work on embedded boundary methods for pressure projection and viscosity. In practice the method’s implementation is remarkably simple, yet it correctly captures the challenging free surface boundary condition that enables simulation of jet buckling phenomena. This work raises a number of questions and directions for future work.
In our numerical study, we observed improved convergence in 3D compared to 2D, and plan to investigate if this truly holds—perhaps beginning by deriving a full-fledged analytic test case as we have done in 2D. The 2D convergence test cases we presented also considered only scenarios where the two different boundary conditions (solid and free surface) do not meet. We suspect that this is an inherently more difficult problem to address, giving rise to issues analogous to those which occur in the presence of sharp boundary features, and our preliminary experiments support this conjecture. Nevertheless, such configurations occur frequently in the buckling examples we have included, illustrating that the method remains stable and provides qualitatively reasonable results.
While the viscous jet buckling example provides a practical validation of our method’s boundary condition enforcement, the current underlying Navier-Stokes simulator is fairly basic. A thorough study of this phenomenon would likely need to consider improved advection and time-splitting methods in place of the first order approaches applied here.
In terms of handling related phenomena, we have noted that our work is closely related to that of Batty et al. who considered considered the simpler Poisson problem for incompressibility with rigid bodies . An obvious extension of the current work would therefore be to consider Stokes flow coupled with fully dynamic deformable structures. Two-phase flow would also be a useful direction to pursue, as would non-Newtonian fluid models. Finally, it would be interesting to consider whether exploiting variational principles in this manner might be useful for handling irregular boundaries in other problems that are commonly discretized on staggered grids. Some potential examples include vorticity-based formulations of fluid flow, diffusion problems, or porous flow.
The following symbols and letters are used throughout the paper, with units given for dimensionful quantities:
- G. K. Batchelor, An Introduction to Fluid Dynamics, 1967.
- Christopher Batty, Florence Bertails, and Robert Bridson, A fast variational framework for accurate solid-fluid coupling, ACM Trans. Graph. (SIGGRAPH), 26 (2007), p. 100.
- Christopher Batty and Robert Bridson, Accurate viscous free surfaces for buckling, coiling, and rotating liquids, in Symposium on Computer Animation, 2008, pp. 219–228.
- Jacob Bedrossian, James H. von Brecht, Siwei Zhu, Eftychios Sifakis, and Joseph Teran, A second order virtual node method for Poisson interface problems on irregular domains, J. Comp. Phys., (in press) (2010).
- Michele Benzi, Gene H. Golub, and Jorg Liesen, Numerical solution of saddle point problems, Acta Numerica2, 14 (2005), pp. 1–37.
- Andrea Bonito, Marco Picasso, and Manuel Laso, Numerical simulation of 3D viscoelastic flows with free surfaces, J. Comp. Phys., 215 (2006), pp. 691–716.
- Mark Carlson, Peter J. Mucha, R. Van Horn, and Greg Turk, Melting and flowing, in Symposium on Computer Animation, 2002, pp. 167–174.
- Robert K.-C. Chan and Robert L Street, A computer study of finite amplitude water waves, J. Comp. Phys., 6 (1970), pp. 68–94.
- Kuan-Yu Chen, Ko-An Feng, Yongsam Kim, and Ming-Chih Lai, A note on pressure accuracy in immersed boundary method for Stokes flow, J. Comp. Phys., In Press (2011).
- J. O. Cruikshank, Low-Reynolds-number instabilities in stagnating jet flows, Journal of Fluid Mechanics, 193 (1988), pp. 111–127.
- J. O. Cruikshank and B. R. Munson, Viscous-fluid buckling of plane and axisymmetric jets, Journal of Fluid Mechanics, 113 (1981), pp. 221–239.
- M. S. Darwish, J. R. Whiteman, and M. J. Bevis, Numerical modelling of viscoelastic liquids using a finite-volume method, Journal of Non-Newtonian Fluid Mechanics, 45 (1992), pp. 311–337.
- Olivier Desjardins, Vincent Moureau, and Heinz Pitsch, An accurate conservative level set/ghost ﬂuid method for simulating turbulent atomization, J. Comp. Phys., 227 (2008), pp. 8395–8416.
- Howard Elman, Multigrid and Krylov subspace methods for the discrete Stokes equations, International Journal for Numerical Methods in Fluids, 22 (1998), pp. 755–770.
- R. Elliot English and Robert Bridson, Animating developable surfaces using nonconforming elements, ACM Trans. Graph. (SIGGRAPH)2, 27 (2008), p. 66.
- Doug Enright, Frank Losasso, and Ron Fedkiw, A fast and accurate semi-Lagrangian particle level set method, Computers and Structures, 83 (2005), pp. 479–490.
- Doug Enright, Duc Nguyen, Frédéric Gibou, and Ron Fedkiw, Using the particle level set method and a second order accurate pressure boundary condition for free surface flows, in Proceedings of the 4th ASME-JSME Joint Fluids Engineering Conference, 2003, pp. 337–342.
- Henrik Fält and Doug Roble, Fluids with extreme viscosity, in SIGGRAPH Sketches, 2003, p. 1.
- Ronald Fedkiw, Tariq Aslam, Barry Merriman, and Stanley Osher, A non-oscillatory Eulerian approach to interfaces in multimaterial flows (the ghost fluid method), J. Comp. Phys., 152 (1999), pp. 457–492.
- Frédéric Gibou, Ron Fedkiw, L.-T. Cheng, and Myungjoo Kang, A second order accurate symmetric discretization of the Poisson equation on irregular domains, J. Comp. Phys., 176 (2002), pp. 205–227.
- J. Guermond, P. Minev, and J. Shen, An overview of projection methods for incompressible flows, Computer Methods in Applied Mechanics and Engineering, 195 (2006), pp. 6011–6045.
- F. H. Harlow and J. E. Welch, Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface, Phys. Fluids, 8 (1965), pp. 2182–2189.
- C. W. Hirt and J. P. Shannon, Free surface stress conditions for incompressible-flow calculations, J. Comp. Phys., 2 (1968), pp. 403–411.
- Hans Johansen and Phillip Colella, A Cartesian grid embedded boundary method for Poisson’s equation on irregular domains, J. Comp. Phys., 147 (1998), pp. 60–85.
- Myungjoo Kang, Ron Fedkiw, and Xu-Dong Liu, A boundary condition capturing method for multiphase incompressible flow, SIAM J. Sci. Comput., 15 (2000), pp. 323–360.
- Randall J. LeVeque and Zhilin Li, Immersed interface methods for Stokes flow with elastic boundaries or surface tension, SIAM J. Sci. Comput., 18 (1997), pp. 709–735.
- Jie Li, Yuriko Y. Renardy, and Michael Renardy, Numerical simulation of breakup of a viscous drop in simple shear flow through a volume-of-fluid method, Phys. Fluids, 12 (2000), pp. 269–282.
- A. Limache, S. R. Idelsohn, R. Rossi, and E. Onate, The violation of objectivity in Laplace formulations of the Navier-Stokes equations, International Journal for Numerical Methods in Fluids, 54 (2007), pp. 639–664.
- A. Limache, P. J. Sanchez, L. D. Dalcin, and S. R. Idelsohn, Objectivity tests for Navier-Stokes simulations: The revealing of non-physical solutions produced by Laplace formulations, Computer Methods in Applied Mechanics and Mechanical Engineering, 197 (2008), pp. 4180–4192.
- Chohong Min and Frédéric Gibou, Geometric integration over irregular domains with application to level-set methods, J. Comp. Phys., 226 (2007), pp. 1432–1443.
- Rajat Mittal and Gianluca Iaccarino, Immersed boundary methods, Annual review of fluid mechanics, 37 (2005), pp. 239–261.
- G. Mompean and M. Deville, Unsteady finite volume simulation of Oldroyd-B fluid through a three-dimensional planar contraction, Journal of Non-Newtonian Fluid Mechanics, 72 (1997), pp. 253–279.
- Yen Ting Ng, Chohong Min, and Frédéric Gibou, An efficient fluid-solid coupling algorithm for single-phase flows, J. Comp. Phys., 228 (2009), pp. 8807–8829.
- B. D. Nichols and C. W. Hirt, Improved free surface boundary conditions for numerical incompressible-flow calculations, J. Comp. Phys., 8 (1971), pp. 434–448.
- Cassio M. Oishi, Murilo F. Tomé, José A. Cuminato, and Sean McKee, An implicit technique for solving 3D low Reynolds number moving free surface flows, J. Comp. Phys., 227 (2008), pp. 7446–7468.
- C. C. Paige and M. A. Saunders, LSQR: An algorithm for sparse linear equations and sparse least squares, TOMS, 8 (1982), pp. 43–71.
- Joseph Papac, Frédéric Gibou, and Christian Ratsch, Efficient symmetric discretization for the Poisson, heat and Stefan-type problems with Robin boundary conditions, J. Comp. Phys., 229 (2010), pp. 875–889.
- Charles S. Peskin, The immersed boundary method, Acta Numerica, 11 (2002), pp. 479–517.
- William E. Pracht, A numerical method for calculating transient creep flows, J. Comp. Phys., 7 (1971), pp. 46–60.
- James W. Purvis and John E. Burkhalter, Prediction of critical Mach number for store configurations, AIAA Journal, 17 (1979), pp. 1170–1177.
- A. Rafiee, M. T. Manzari, and M. Hosseini, An incompressible SPH method for simulation of unsteady viscoelastic free surface flows, International Journal of Non-Linear Mechanics, 42 (2007), pp. 1210–1223.
- Nick Rasmussen, Doug Enright, Duc Nguyen, Sebastian Marino, N. Sumner, Willi Geiger, Samir Hoon, and Ron Fedkiw, Directable photorealistic liquids, in Symposium on Computer Animation, 2004, pp. 193–202.
- Avi Robinson-Mosher, Craig Schroeder, and Ron Fedkiw, A symmetric positive definite formulation for monolithic fluid structure interaction, In submission, (2010).
- Avi Robinson-Mosher, Tamar Shinar, Jon Gretarsson, Jonathan Su, and Ronald Fedkiw, Two-way coupling of fluids to rigid and deformable solids and shells, ACM Trans. Graph. (SIGGRAPH), 27 (2008), p. 46.
- Doug Roble, Nafees bin Zafar, and Henrik Falt, Cartesian grid fluid simulation with irregular boundary voxels, in SIGGRAPH Sketches, 2005, p. 138.
- P. A. Stewart, N. Lay, Mark Sussman, and Mitsuhiro Ohta, An improved sharp interface method for viscoelastic and viscous two-phase flows, J. Sci. Comput., 35 (2008), pp. 43–61.
- Mark Sussman and Mitsuhiro Ohta, Improvements for calculating two-phase bubble and drop motion using an adaptive sharp interface method, Fluid Dynamics and Materials Processing, 3 (2007), pp. 21–36.
- Murilo F. Tomé, L. Grossi, Antonio Castelo, José A. Cuminato, Norberto Mangiavacchi, Valdemir G. Ferreira, F. S. de Sousa, and Sean McKee, A numerical method for solving three-dimensional generalized Newtonian free surface flows, Journal of Non-Newtonian Fluid Mechanics, 123 (2004), pp. 85–103.
- Murilo F. Tomé and Sean McKee, GENSMAC: A computational marker and cell method for free surface flows in general domains, J. Comp. Phys., 110 (1994), pp. 171–186.
- M F Tomé and Sean Mckee, Numerical simulation of viscous flow: Buckling of planar jets, International Journal for Numerical Methods in Fluids, 29 (1999), pp. 705–718.
- Justin Wan and Xu-Dong Liu, A boundary condition-capturing multigrid approach to irregular boundary problems, SIAM J. Sci. Comput., 25 (2004), pp. 1982–2003.
- Y. C. Zhou, Shan Zhao, Michael Feig, and G. W. Wei, High order matched interface and boundary method for elliptic equations with discontinuous coefficients and singular sources, J. Comp. Phys., 213 (2005), pp. 1–30.