- Email: [email protected]
Using integral equations and the immersed interface method to solve immersed boundary problems with stiff forces
Computers & Fluids 38 (2009) 266–272
Contents lists available at ScienceDirect
Computers & Fluids j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m p fl u i d
Using integral equations and the immersed interface method to solve immersed boundary problems with stiff forces q Anita T. Layton Department of Mathematics, Duke University, Durham, NC 27708, USA
a r t i c l e
i n f o
Article history: Received 17 August 2007 Received in revised form 20 February 2008 Accepted 22 February 2008 Available online 2 March 2008
a b s t r a c t We propose a fast, explicit numerical method for computing approximations for the immersed boundary problem in which the boundaries that separate the fluid into two regions are stiff. In the numerical computations of such problems, one frequently has to contend with numerical instability, as the stiff immersed boundaries exert large forces on the local fluid. When the boundary forces are treated explicitly, prohibitively small time-steps may be required to maintain numerical stability. On the other hand, when the boundary forces are treated implicitly, the restriction on the time-step size is reduced, but the solution of a large system of coupled non-linear equations may be required. In this work, we develop an efficient method that combines an integral equation approach with the immersed interface method. The present method treats the boundary forces explicitly. To reduce computational costs, the method uses an operator-splitting approach: large time-steps are used to update the non-stiff advection terms, and smaller substeps are used to advance the stiff boundary. At each substep, an integral equation is computed to yield fluid velocity local to the boundary; those velocity values are then used to update the boundary configuration. Fluid variables are computed over the entire domain, using the immersed interface method, only at the end of the large advection time-steps. Numerical results suggest that the present method compares favorably with an implementation of the immersed interface method that employs an explicit time-stepping and no fractional stepping. Ó 2008 Elsevier Ltd. All rights reserved.
1. Introduction Frequently in problems of fluid flow or wave propagation, a boundary between different regions exerts a force on the material, or an interface separates regions of different material properties. Substantial effort has been directed to developing computational techniques for simulating moving boundaries or interfaces within a viscous fluid domain. The most notable example is perhaps the immersed boundary method, developed by Peskin [26] for solving the full incompressible Navier–Stokes equations with moving boundaries. The immersed boundary method was originally developed for studying blood flow through a beating heart [25]. In the immersed boundary method, moving boundaries are represented by means of Lagrangian markers, at which boundary forces are computed. These forces are then transferred to an underlying Cartesian grid via a regularization using discrete (smoothed) delta functions, and the Navier–Stokes equations are solved on the Cartesian grid. In principle, the Lagrangian representation of the immersed boundary allows complicated boundaries to be incorporated without significant difficulty. The immersed boundary methq
This work was supported in part by the National Science Foundation, Grant DMS0715021. E-mail address: [email protected]
0045-7930/$ - see front matter Ó 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.compfluid.2008.02.003
od has since been applied to a wide variety of problems [1,3,7,8]. Because the approximate Dirac delta functions, which transfers the singular boundary forces onto the underlying fluid, typically have OðhÞ support, where h denotes the spatial mesh spacing, the immersed boundary method does not capture the jump discontinuity in solution (e.g., pressure) but rather approximates the solution as a continuous function with a large gradient. In general, the immersed boundary method computes approximations with firstorder spatial accuracy. An alternative approach that captures jumps in solution and its derivatives sharply and that generates approximations with second-order spatial accuracy is the immersed interface method developed by LeVeque and Li [17,18]. The immersed interface method is similar to a method developed earlier by Mayo for solving elliptic problems on an irregular domain [21,22]. Both methods are Cartesian grid methods with second-order spatial accuracy (though higher-order immersed interface methods have recently been developed [12,20]), with the key idea being the incorporation of known jumps in the solution or its derivatives into the finite difference schemes. One difficulty that researchers of immersed boundary problems have to contend with is that boundary forces arising from stiff immersed boundaries frequently impose a severe restriction on timestep size in order to maintain numerical stability. When the
A.T. Layton / Computers & Fluids 38 (2009) 266–272
267
boundary forces are treated explicitly, stiff forces impose a prohibitive restriction on the time-step, thereby increasing the computational costs. Thus, much effort has been invested in developing implicit and semi-implicit versions of the immersed boundary method and related methods, e.g. [9,16,18,23,28]. Nonetheless, the treatment of stiff forces remains a challenge. Because of the coupling between boundary forces and the underlying fluid, an implicit or semi-implicit treatment of the boundary forces requires the solution of a large system of non-linear equations via an iterative procedure, the convergence of which can be a concern. Perhaps owing to that difficulty, the majority (though by no means all) of the implementations of the immersed boundary and interface methods are explicit ones. As a result, when the computations become unstable, one is left to wonder whether the problem is attributable to a programming error, or whether one should try a yet smaller time-step, or in the case of an implicit method, a better initial guess for the iterations. The principal goal of this work is to develop an efficient timestepping approach for advancing solutions of immersed boundary problems for the case where the boundary forces are stiff. The proposed method offers the following advantages:
boundary is considered here, the generalization to multiple, closed immersed boundaries is straightforward. We seek to compute a solution on a bounded computational domain X ¼ C þ X þ Xþ , where Xþ 2 X1þ . The immersed boundary C is assumed to be elastic, so that when it is distorted from its relaxed state by stretching or compressing, it exerts a force on the local fluid. The model equations are the Navier–Stokes equations with forcings: ou ð1Þ þ u ru ¼ rp þ lr2 u þ f; q ot
Boundary forces are treated explicitly, thus avoiding the solution of a large system of coupled non-linear equations. With the forces treated explicitly, a small time-step may be needed, which can lead to large overall computational costs. The proposed method keeps the overall computational costs low by means of an integral equation, which provides a fast way of updating the boundary forces.
for i ¼ 1; 2. L0 is the length of C at its relaxed state; x 2 X, and Xðs; tÞ gives the location of C at time t, parametrized by the arclength s at the unstretched state for 0 6 s 6 L0 ; F i is the force strength at point s; and d is the two-dimensional delta function. Note that d is not the approximate (smoothed) delta function used in the immersed boundary method, and that f is singularly supported along C. We consider an elastic band, where the force strength F is given by
r u ¼ 0;
ð2Þ
where u ¼ ðu; vÞ denotes the fluid velocity; p is the pressure; q and l are the fluid density and viscosity, respectively, assumed to be constant; and f is a force term singularly support along C, with the xand y-components denoted by f1 and f2 , respectively. The forces f1 and f2 are elastic or tension forces that arise from C being stretched or compressed. These singular forces can be written as Z L0 fi ðx; tÞ ¼ F i ðs; tÞdðx Xðs; tÞÞds; ð3Þ 0
Our approach is based on a combination of the integral solution of the unsteady Stokes equations and the immersed interface method. Boundary integrals have been used to remove the stiffness from interfacial flows with surface tension [10]. Indeed, the idea of combining the boundary integral method and immersed interface method has been previously applied in other fluid problems. For example, using that approach, Beale and Layton solved the Stokes problem with discontinuous viscosity [2,14], and Biros et al. computed solutions for the Stokes equations [4] and for the Navier– Stokes equations [5] on irregular domains. In this work, the integral solution of the unsteady Stokes equations is used to approximate fluid velocity local to the immersed boundary. That boundary velocity is then used to update the boundary configuration and forces. Numerical results suggest that the efficiency of the proposed method compares favorably to an implementation of the immersed interface method that employs explicit time-stepping. 2. Numerical method We consider an evolving, closed boundary C that separates the two-dimensional domain X into two subdomains, X1þ (exterior) and X (interior); see Fig. 1. Note that although only one immersed
Ω+
Γ Ω−
Fig. 1. The immersed boundary model configuration.
Fðs; tÞ ¼
o ðTðs; tÞsðs; tÞÞ; os
with the tension Tðs; tÞ given by oX Tðs; tÞ ¼ T 0 1 : os
ð4Þ
ð5Þ
The tension coefficient T 0 describes the elastic properties of the immersed boundary and is assumed to be a constant in this model. The vector tangential to C is given by sðs; tÞ, where sðs; tÞ ¼
oX=os : joX=osj
ð6Þ
Thus, the force density can be computed directly from the location X of the boundary C. Note that at the relaxed state, oX=os ¼ 1 and tension vanishes. In the spatial discretization, let hx and hy denote the spatial grid interval along the x- and y-axis, respectively, and ðihx ; jhy Þ be the ith and jth spatial grid points along the x- and y-axis, respectively, where i; j ¼ 0; 1; 2; . . . ; N. For notational simplicity, let hx ¼ hy h. The position of the immersed boundary at t is represented by a set of boundary markers XðtÞ ðX k ðtÞ; Y k ðtÞÞ, for k ¼ 0; 1;2; .. .; N k , where ðX 0 ðtÞ;Y 0 ðtÞÞ ¼ ðX Nk ðtÞ; Y Nk ðtÞÞ because the boundary is assumed to be a simple closed curve. The kth boundary marker approximates ðXðsk ;tÞ; Yðsk ;tÞÞ, where sk ¼ kL0 =N k ; that is, the boundary markers are chosen such that in an unstretched or relaxed state, they are equally spaced. The present method updates boundary forces by means of the integral solution of the unsteady Stokes equations (see below). To apply that solution to the Navier–Stokes equation, we first discretize the momentum Eq. (1) in time using a semi-implicit approach, with the inertia and force terms treated explicitly. An operator-splitting approach is also applied, where the stiff forces are updated using a smaller time-step, which we refer to as the forcing substep, than the advection or overall time-step used to update the advection terms. Suppose we take N m (small) forcing sub-
268
A.T. Layton / Computers & Fluids 38 (2009) 266–272
steps Dt m per (large) advection time-step Dtn ; thus, Dt n ¼ N m Dt m . Let tm and t n denote the forcing and advection time-levels, respectively, where t m ¼ tn þ mDt m ¼ tn þ mDtn =N m , for m ¼ 1; 2; . . . ; N m . With a first-order time-differencing, the discretized form of Eq. (1) is given by frac1Dtm umþ1 þ mþ1
ru
1 l 1 m fm u un run þ ; rpmþ1 r2 umþ1 ¼ q q q Dtm
¼ 0:
i;k¼1;2
ð7Þ ð8Þ
n
We then select a point x0 2 Xþ , and integrate Eq. (11) over X . Using the divergence theorem to convert the surface integral over X into a boundary integral over C, we obtain X Z Gij ð^ xÞrij ðxÞ Tijk ð^ xÞui ðxÞ nk dx
n
Note that the advection term u ru is evaluated at t n . A simple first-order time-stepping is used, in part, because the operatorsplitting is first-order. Possible extensions to higher-order time discretization are noted in Section 4. Because the boundary forces are treated explicitly, the size of Dtm may be severely restricted to maintain numerical stability. For the proposed method to be efficient, the computation cost in advancing the solution by one Dtm must be sufficiently cheap. To achieve that, we note that to update the boundary forces, one must update the boundary position; to update the boundary position, one needs to update the fluid velocity local to the boundary—but nowhere else. With that observation, we propose an efficient time-stepping method that updates only the boundary velocity during each Dtm using boundary and volume integrals; fluid velocity is then computed (less frequently) in the entire domain at each advection step Dtn .
¼
C
X Z
i¼1;2
X
Gij ð^ xÞbi ðxÞdx;
ð14Þ
Next, we select a point x0 2 C, and integrate Eq. (11) over X [24]: X1 Z uj ðx0 Þ ¼ Gij ð^ xÞrij ðxÞ Tijk ð^ xÞui ðxÞ nk dx pl C i¼1;2 Z 1 Gij ð^ xÞbi ðxÞdx : ð15Þ þ 2pl X To express the integral solution uj ðx0 Þ in terms of the boundary forces f, we multiply Eq. (14) by 1=ðplÞ and take the difference between Eq. (15) and the resulting equation, which yields Z X 1 Z 1 uj ðx0 Þ ¼ Gij ð^ xÞfi ðxÞdx Gij ð^ xÞbi ðxÞdx : ð16Þ pl C 2pl X i¼1;2 Thus, given boundary forces f, boundary velocity, denoted uC , can be obtained by computing the boundary integrals and the volume integrals on the right side of (16). 2.2. Evaluation of boundary and volume integrals
2.1. Integral solution To derive the method, we drop the superscripts in Eq. (7) and let k2 q=lDtm . Using that notation, we rewrite Eq. (7) as ðr2 k2 Þu rp ¼ b;
ð9Þ
where p is a modified pressure that contains a 1=l factor. The left side contains unknown variables at t mþ1 , and the right side b um =Dt m þ un run contains variables at previous time-levels t n and t m . To develop an integral formulation for Eq. (9), which is also known as Brinkman’s equation, we rewrite Eq. (9) in terms of the stress tensor r: r r k2 u ¼ b:
ð10Þ
Let r and u be a solution of the inhomogeneous Eq. (10), and let r and u be a solution of the homogeneous equation with b ¼ 0. Then by Green’s identity, r, r , u, and u satisfy X o X ui rij ui rij ¼ ui bi : ð11Þ oxj i¼1;2 i¼1;2 We identify r and u with the fundamental solution of Brinkman’s equation, denoted T and G, respectively. In two-dimensional free space, Gij and Tijk are given by ! K 1 ðkrÞ 1 Gij ð^ xÞ ¼ dij K 0 ðkrÞ þ kr ðkrÞ2 ! ^ 2K 1 ðkrÞ 2 xi ^ xj ; ð12Þ þ 2 K 0 ðkrÞ þ r kr ðkrÞ2 ! 4K 1 ðkrÞ 4 xÞ ¼ ^ xj dik 1 þ 2K 0 ðkrÞ þ Tijk ð^ kr ðkrÞ2 ! 4K 1 ðkrÞ 4 xi djk Þ krK 1 ðkrÞ þ þ ð^ xk dij þ ^ kr ðkrÞ2 ! ^ 16K 1 ðkrÞ 16 xi ^ xj ^ xk ; ð13Þ 8K 0 ðkrÞ 2krK 1 ðkrÞ þ 2 þ r kr ðkrÞ2 where ^x x0 x, r j^xj, and K 0 and K 1 are modified Bessel functions of the second kind [24].
Normally, the evaluation of a volume integral can be computaR tionally expensive. However, in the case of X Gij bi dx, which is the second term in the summation on the right side of (16), one might note that Gij decays rapidly away from the poles (i.e., as a function of kr, see Eq. (12)). Thus, in practice, the volume integral can be approximated by a narrow band D neighboring C. In particular, note that k scales as 1=ðlDt m Þ, thus Gij decays faster at lower viscosity l and smaller time-step Dtm . For stiff forces, Dtm must be sufficiently small to maintain numerical stability, rendering the decay of Gij extremely fast. Consequently, the band D along which the volume integral is computed can be made narrow without significantly sacrificing accuracy. A challenge in computing both the boundary and volume integrals in (16) is that when k is large, Gij has large derivatives. (Spex, the distance from the pole, cifically, Gij decreases rapidly as ^ increases.) To reduce quadrature errors, the adaptive trapezoidal rule, which employs an increasingly refined grid near the poles, was used. One technical complication in using an adaptive quadrature grid is that the boundary force f is known only at the boundary markers but not at every quadrature nodes in the refined grid. Although spatial interpolation could be used to estimate f values at the quadrature nodes, one notes that, compared to Gij , f likely does not vary significantly from one quadrature node to its neighbor. Thus, in the evaluation of the boundary integral, we allow Gij to vary at every quadrature node (see below), but assume that f is constant within the interval ½ðXk1 þ Xk Þ=2; ðXk þ Xkþ1 Þ=2. This approach avoids the additional cost of interpolations without significantly sacrificing numerical accuracy. A similar approach is used in the approximation of the volume integral. An adaptively refined grid is used to compute the quadrature that approximates the volume integral over D. The values of b, which varies slowly compared to Gij , are known at grid points but not at each quadrature node. Following the approach used in the evaluation of the boundary integral, we assume that b is constant within the h h square centered at each grid point. Another issue here is that away from C, um (the first term in b) is not known. This is because u is updated at tm only along C, thus elsewhere u is known only at previous t n levels. The values of um away from C can be estimated by extrapolating in time using values at tn and
269
A.T. Layton / Computers & Fluids 38 (2009) 266–272
t n1 . Extrapolation may introduce inaccuracy, but because Gij decays rapidly away from C, many of those terms in the integral that require extrapolation, which are located away from C, likely have a small contribution. The most expensive step in the computation of uC is the evaluation of Green’s function using Eq. (12). Because Eq. (12) involves Bessel functions, the evaluation of which is computationally nontrivial, it is not clear that, in a naïve implementation, computing the boundary velocity using the integral Eq. (16) is less expensive than computing fluid velocity over X. To speed up the evaluation of Gij , we place quadrature nodes at fixed distances from the boundary markers; i.e., the adaptively refined grid used in the trapezoidal rule is the same for every time-step throughout the computation. Specifically, the quadrature nodes used in estimating boundary integral are chosen to be points along C that are located at distances of r q from each boundary marker Xk , for q ¼ 1; 2; . . . ; N q . The locations of the quadrature nodes in the computation of the volume integral are chosen analogously. Because the distances r q between the quadrature nodes and the boundary markers remain fixed, the two terms in Eq. (12) involving the Bessel functions, i.e., ! ! K 1 ðkr q Þ 1 1 2K 1 ðkrq Þ 2 ; ; K ðkr Þ þ K 0 ðkrq Þ þ 0 q kr q r 2q kr q ðkrq Þ2 ðkrq Þ2 can be precomputed. The evaluation of the boundary and volume integrals requires OðN q N k Þ and OðN 2q N k Þ work, respectively. On the other hand, if fluid velocity and pressure were computed everywhere to advance the boundary configuration, as in an implementation of the original immersed interface method (i.e., without the fractional stepping and the integral approach), then the computations would take OðN 2 log NÞ work if fast Fourier transforms are used. If one assumes that N q scales as N, then a comparison of the relative efficiency of those two approaches depends, in large part, on the values of N q and N. The efficiency of the two approaches will be compared using a numerical test below. 2.3. The algorithm Once the boundary velocity um C is obtained by evaluating Eq. (16), the immersed boundary position can be updated. As previously noted, the position X of C is determined by the location of a set of boundary markers along C; and at a given time t, C is represented by a cubic spline through those markers [18]. To advance the boundary by one time-step, we move the boundary markers at the same velocity as the local fluid: o Xðs; tÞ ¼ uðXðs; tÞ; tÞ: ot
method because the method is well-known among the immersed boundary community.) The boundary and volume integrals are not used to compute fluid variables everywhere at the advection time-level t n because such computations cost OðN 2q N2 Þ time and are likely more expensive compared to the immersed interface method. The algorithm for advancing the solution from tn to t nþ1 is summarized below: (1) For m ¼ 1; 2; . . . ; N m 1. (a) Compute boundary velocity uðXk ; t m Þ along C using the integral Eq. (16). using (18). (b) Update boundary position Xmþ1 k mþ1 using (3) and (4). (c) Update boundary forces f (2) Express jumps in solution and derivatives across C in terms N nþ1 of boundary forces f m ¼ f . Compute fluid solution unþ1 and pnþ1 over X using the immersed interface method [19].
3. Numerical examples In this section we present numerical results for two examples. Using these examples, we aim to demonstrate the second-order accuracy of the method in space, and to assess its efficiency. All calculations reported below were performed using Fortran programs, which were executed in double precision on a computer with an Intel Core 2 Duo 6400 2.13 GHz processor, 2 MB of cache, and 1 GB of RAM. 3.1. Example 1 The first example was used by Tu and Peskin to test their immersed boundary method [28], and by LeVeque and Li to test their immersed interface method [18]. The initial boundary is an ellipse with major and minor axes set to a ¼ 0:7 and b ¼ 0:5, respectively (see Fig. 2, dashed curve labelled ‘‘Initial”). The unstretched boundary was set to a circle with radius r 0 ¼ 0:5. The tension coefficient T 0 was set to 0.5. The fluid density q was set to 1. At t ¼ 0, the fluid
1
ð17Þ
0.5
Forward Euler is used to advance the boundary as follows: m Xmþ1 ¼ Xm k k þ DtuC k :
ð18Þ
Once a new boundary configuration has been obtained, the boundmþ1 can then be computed using (3) and (4), and a new ary force f forcing step Dt m can be taken. After N m force substeps Dt m , the fluid velocity and pressure can be computed over the entire domain X using the immersed interface method [19]. Jumps in the fluid variables and their derivatives are computed from the boundary forces and their derivatives; those jumps are then incorporated into a finite-difference scheme; for details see [19]. One technical detail that should be noted is that in this implementation, jumps in variables and their derivatives are computed along coordinate lines, as in a previous study [15] and in Mayo’s method [21,22] but unlike the immersed interface method which computes those jumps in the directions normal and tangential to C. (Nonetheless, we refer to the present method as a variation of the immersed interface
Final
0
Initial -0.5
-1
-1
-0.5
0
0.5
1
Fig. 2. Boundary configuration at initial time (dashed line) and at steady state (solid dark line).
270
A.T. Layton / Computers & Fluids 38 (2009) 266–272
was initialized to be at rest. The computational domain was set to ½2:5; 2:5 ½2:5; 2:5. At steady state, the initial ellipse converged to a circle with rapffiffiffiffiffi ffi dius re ¼ ab 0:6124, which is larger than the unstretched boundary but which has the same area as the initial ellipse, owing to the incompressibility of the enclosed fluid. At steady state, fluid velocity vanished everywhere; p attained constant values inside and outside the boundary, with a jump ½p < 0 (because the boundary was initialized to a stretched state and because fluid velocity vanished at steady state). We simulated the motion of the fluid and of the boundary for t 2 ½0; 5 using an overall time-step Dtn ¼ 2 103 , together with a smaller forcing substep of Dtm ¼ 2 104 . For a given Dtm , the xj depends in rate at which kGð^ xÞk1 decreases as a function of j^ large part on the fluid viscosity l. As an example, consider the case x2 > 0. The values of kGð^ xÞk1 for differing values of l in which ^ x1 ¼ ^ and r ¼ j^ xj are shown in Fig. 3. To further compare the boundary layers for differing viscosity values, for example, l ¼ 105 ;104 ; xj > 0:04473, 0.1414, 103 ; 102 , and 101 , kGk1 < 106 for j^ 0.4472, 1.1414, and 4.472, respectively. Thus, the boundary layer becomes increasingly thin as l decreases (see Fig. 3). Although the width of the band D over which the volume integral is computed can be decreased for smaller values of l, the gradient of G increases as l decreases. Thus, to maintain quadrature accuracy, we used approximately the same number of quadrature nodes N q , for a given spatial resolution N, for all l values. To assess the spatial accuracy of the method, we scaled N k , N q , and N simultaneously. Specifically, we let N k ¼ N, and lowered the error tolerance of the adaptive trapezoidal rule by a factor of 4 as N doubled, thereby increasing N q . Volume accuracy results for the case of l ¼ 104 are shown in Table 1 for N ¼ 80, 160, 320, and 640. The volume errors are given by the relative difference in the area bounded by C, between exact solution and the computed solution at approximate steady state. These results exhibit second-order convergence at sufficiently large N values. At N ¼ 80, there was substantial water loss in X , an error that can be attributed to the under-resolved boundary layer neighboring C. Nonetheless, with a spatial grid that is sufficiently refined to resolve the boundary layer, volume is approximately conserved. For example, for N P 160, the relative volume error is 6 0:3%. Analogous results were obtained also for l ¼ 103 and 105 . We also estimated relative errors in the fluid velocity u ðu; vÞ. Because the exact solution is not known for this problem, the approximation computed using a refined spatial discretization of N ¼ 1280 is used as a refer-
Table 1 Convergence results. A ð1Þ, steady-state area enclosed by C N
A ð1Þ
errorA
kerroru k2
80 160 320 640
0.71382151 0.78297285 0.78480984 0.78524711
9.113E2 3.088E3 7.491E4 1.923E4
1.261E1 4.686E2 1.384E2 3.771E3
The relative errors in A ð1Þ A ð1Þ ¼ 0:7853981635. Results convergence.
are computed using the exact solution show approximate second-order spatial
ence solution. The results, given in L2-norm, exhibit approximately second-order accuracy; see Table 1. 3.2. Example 2 The second test involves a more complicated initial boundary configuration, which is given in polar coordinates by r ¼ 0:8þ 0:3 sin 8h. The relaxed boundary is the circle with radius r0 ¼ 0:3. The initial and relaxed boundary configurations are shown in Fig. 4 (dashed and dotted curves, respectively). The larger curvature, compared to the first example, renders this problem significantly stiffer, inasmuch as the boundary tension force increases with the curvature. The density and viscosity coefficients of the fluid were taken to be q ¼ 1 and l ¼ 104 , respectively. The tension coefficient T 0 was set to 0.5. The fluid was initialized to be at rest. The spatial domain was set to ½1:2; 1:2 ½1:2; 1:2. The system was integrated using an overall time-step Dtn of 0.01, and 40 forcing substeps Dt m ¼ 2:5 104 , the largest forcing substep size that can be taken without incurring numerical instability. A spatial grid of N ¼ 640 subintervals was used. Nine quadrature nodes were used in each dimension to approximate the integrals, i.e., N q ¼ 13. With these parameters, we compared the computational costs of the present method to an explicit implementation of the traditional immersed interface method that does not involve fractional stepping. To maintain numerical stability, the immersed interface method requires a time-step Dt of 2:5 104 , which is equal to the forcing substep Dtm of our method. At each time-step, the immersed interface method solved the Navier–Stokes equations (7) and (8), and advanced fluid and boundary variables over the entire computational domain X. For this example, and with the time-steps specified, the present meth-
1.2 2
0.8
|G|
1.5
μ = 10
0.4
−1
μ = 10
1
−2
0
μ = 10
−3
μ = 10
0.5
-0.4
−4
μ = 10
−5
-0.8 0 0
0.0025
0.005
0.0075
0.01
0.0125
-1.2 -1.2
-0.8
-0.4
0
0.4
0.8
1.2
r Fig. 3. kGðrÞk1 as a function of r, for differing l values. Dtm ¼ 2 104 ; ^ x1 ¼ ^ x2 > 0.
Fig. 4. Boundary configuration at initial time (dashed line) and at t ¼ 0:1 (solid line). Dotted line shows relaxed configuration.
A.T. Layton / Computers & Fluids 38 (2009) 266–272
od was approximately 60% faster, measured in CPU time, than the immersed interface method. 4. Discussion In this work, we combined an integral approach with the immersed interface method for an efficient solution of immersed boundary problems with stiff forcings. The method integrates the boundary forces explicitly, an approach that requires small timesteps, compared to an implicit treatment of stiff force terms, but that avoids the need to solve a system of coupled non-linear equations. To reduce computational costs, fluid velocities values are computed only along the immersed boundary using an integral solution of the unsteady Stokes equations. Numerical test results suggest that the efficiency of the proposed method compares favorably to an explicit non-fractional stepping implementation of the original immersed interface method. A direct comparison between the present method (which we refer to below as the ‘‘integral immersed interface method”) and the immersed interface method is difficult to make for a general case, inasmuch as the relative merits of the two methods depend on a number of factors including spatial resolution, stiffness of the problem, accuracy requirement, etc. Moreover, using CPU time as a measure of computational costs has the drawback that metric depends not only on the problem at hand, but also on the particular implementation and the computer architecture. Thus, below we analyze the computational costs of the two methods as functions of spatial resolution. As previously noted, the evaluation of the volume integral at a forcing substep of the integral immersed interface method requires OðN 2q N k Þ work. An immersed interface method using fast Fourier transforms requires OðN 2 log NÞ work for each time-step. If both N q and N k scales as N, then technically the computational cost of the integral immersed interface method is OðN 3 Þ, which should be higher than the original immersed interface method at sufficiently large N’s. Nonetheless, in practice, N q N; thus, at moderate N values, the integral immersed interface method can be faster, as we have shown in Section 3.2. Furthermore, for large N q , specialized fast integration methods for nearly singular functions, e.g. [13,27], may be used. Compare to an implementation of the immersed interface method that uses fast Fourier transforms, the integral immersed interface method may be better suited for parallelization. In two- or three-dimensional computations, Fast Fourier transforms require transposing data, a step which gives rise to global communication among processors. In contrast, computing the integrals in the forcing substeps of the integral immersed interface method requires mostly data local to a processor. Because it requires fewer global inter-processor communications, a parallel implementation of the integral immersed interface method likely has a lower overall computational time, compared to a parallel implementation of the immersed interface method using fast Fourier transforms. The methodology of the present method can also be applied to the immersed boundary method, where the boundary forces are not singularly supported along C but are instead spread onto the fluid grid using approximate delta functions [26]. Assuming that the approximate delta function is compactly supported over D Dþ [ D , where Dþ ½C; C þ dn and D ½C dn; C, then the boundary integral takes the form Z Z X 1 1 Gij ð^ xÞfi ðxÞdx Gij ð^ xÞbi ðxÞdx ð19Þ uj ðx0 Þ ¼ 2pl D 2pl X i¼1;2 where fi above is the smoothed forcing term. A limitation of the integral immersed interface method might be that it requires Green’s function of the unsteady Stokes equa-
271
tion, and that the form of Green’s function depends on the number of spatial dimensions, the shape of the fluid domain, boundary conditions, etc. Thus, Green’s function and, as a result, the integral solution may not be available for arbitrary domain and boundary conditions. Nonetheless, an approximate integral solution may be obtained using the Green’s function that corresponds to a different domain or boundary conditions. The impact of using an approximate Green’s function on numerical stability and accuracy will need to be accessed. The present method is second-order in space, but only first-order in time. We opted for a first-order time-stepping scheme because of the simplicity of a first-order operator-splitting, and because of the likelihood (not verified) that order reduction, which is attributable to the stiffness of the boundary, would restrict solution accuracy to first order regardless. In addition, without incorporating temporal correction terms as in [29], a second-order time-stepping method such as that used in an implementation of the immersed interface method [19] would have yielded approximations with only first-order accuracy in time for problems with moving interfaces. Nonetheless, a version of the method that is second or higher order in time may be desirable. Strang splitting can be used to develop a second-order time-stepping method with operator-splitting. For third- or higher-order methods with operator-splitting, the multi-implicit spectral deferred correction (SDC) approach may be used to generate a high-order solution using a low-order time integration method and a low-order splitting, using an iterative procedure [6]. To attain high-order temporal accuracy for a stiff problem, we note that a version of the SDC method that uses the GMRES Krylov subspace method has been shown in the fully-explicit and fully-implicit cases to eliminate the order reduction phenomenon previously observed for stiff ODE systems [11]. That approach, together with appropriate temporal correction terms [29], may be used to generate high-order approximations in an immersed boundary problem with stiff forces. References [1] Arthurs KM, Moore LC, Peskin CS, Pitman EB, Layton HE. Modeling arteriolar flow and mass transport using the immersed boundary method. J Comput Phys 1998;147:402–40. [2] Beale T, Layton AT. On the accuracy of finite difference methods for elliptic problems with interfaces. Commun Appl Math Comput Sci 2006. [3] Beyer RP. A computational model of the cochlea using the immersed boundary method. J Comp Phys 1992;98:145–62. [4] Biros G, Ying L, Zorin D. An embedded boundary integral solvers for the Stokes equations. J Comput Phys 2004;193:317–48. [5] Biros G, Ying L, Zorin D. An embedded boundary integral solvers for the unsteady incompressible Navier–Stokes equations. Technical report, Ms-CIS04-42, 2004. [6] Bourlioux A, Layton AT, Minion ML. Higher-order multi-implicit spectral deferred correction methods for problems of reacting flow. J Comput Phys 2003;189:651–75. [7] Dillon R, Fauci L, Gaver III D. A microscale model of bacterial swimming, chemotaxis, and substrate transport. J Theor Biol 1995;177:325–40. [8] Fauci L, McDonald A. Sperm mobility in the presence of boundaries. Bull Math Biol 1995;57:679–99. [9] Fauci LJ, Folgelson AL. Truncated newton methods and the modeling of complex immersed elastic structures. Commun Pure Appl Math 1993;66:787–818. [10] Hou TY, Lowengrub JS, Shelley MJ. Removing the stiffness from interfacial flows with surface tension. J Comput Phys 1994;114:312–38. [11] Huang J-F, Jia J, Minion ML. Accelerating the convergence of spectral deferred correction methods. J Comput Phys 2006;214(2):633–56. [12] Ito K, Kyei Y, Li Z. Higher-order Cartesian grid based finite difference schemes for elliptic equations on irregular domains. SIAM J Sci Comput 2005;27(1):346–67. [13] Kapur S, Rokhlin V. High-order corrected trapezoidal quadrature rules for singular functions. SIAM J Numer Anal 1997;34(3):1331–56. [14] Layton AT. An efficient numerical method for the two-fluid stokes equations with a moving boundary. Comput Methods Appl Mech Eng [in press]. [15] Layton AT. Modeling water transport across elastic boundaries using an explicit jump method. SIAM J Sci Comput 2007;28(6):2189–207.
272
A.T. Layton / Computers & Fluids 38 (2009) 266–272
[16] Lee L, LeVeque RJ. An immersed interface method for the incompressible Navier–Stokes equations. SIAM J Sci Comp 2003;25:832–56. [17] LeVeque RJ, Li Z. The immersed interface method for elliptic equations with discontinuous coefficients and singular sources. SIAM J Numer Anal 1994;31:1019–44. [18] LeVeque RJ, Li Z. Immersed interface methods for Stokes flow with elastic boundaries or surface tension. SIAM J Sci Comput 1997;18(3):709–35. [19] Li Z, Lai M-C. The immersed interface method for the Navier–Stokes equations with singular forces. J Comput Phys 2001;171:822–42. [20] Linnick MN, Fasel HF. A high-order immersed interface method for simulating unsteady incompressible flows on irregular domains. J Comput Phys 2005;204:157–92. [21] Mayo A. The fast solution of Poisson’s and the biharmonic equations on irregular regions. SIAM J Numer Anal 1984;21:285–99. [22] Mayo A, Greenbaum A. Fast parallel iterative solution of Poisson’s and the biharmonic equations on irregular regions. SIAM J Sci Comput 1982;13:101–18.
[23] Mayo A, Peskin CS. An implicit numerical method for fluid dynamics problems with immersed elastic boundaries. In: Cheer AY, von Dam CP, editors. Fluid Dynamics in Biology. Providence, RI: AMS; 1993. p. 261–78. [24] Occhialini JM, Muldowney GP, Higdon JJL. Boundary integral/spectral element approaches to the Navier–Stokes equations. Int J Numer Methods Fluids 1992;15:1316–81. [25] Peskin CS. Numerical analysis of blood flow in the heart. J Comput Phys 1977;25:220–52. [26] Peskin CS. The immersed boundary method. Acta Numer 2002:1–39. [27] Strain J. Locally-corrected multidimensional quadrature rules for singular functions. SIAM J Sci Comput 1995;6(4):992–1017. [28] Tu C, Peskin CS. Stability and instability in the computation of flows with moving immersed boundaries: a comparison of three methods. SIAM J Sci Statist Comput 1992;13:1361–76. [29] Xu S, Wang ZJ. Systematic derivation of jump conditions for the immersed interface method in three-dimensional flow simulation. SIAM J Sci Comp 2006;27:1948–80.
Contents lists available at ScienceDirect
Computers & Fluids j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m p fl u i d
Using integral equations and the immersed interface method to solve immersed boundary problems with stiff forces q Anita T. Layton Department of Mathematics, Duke University, Durham, NC 27708, USA
a r t i c l e
i n f o
Article history: Received 17 August 2007 Received in revised form 20 February 2008 Accepted 22 February 2008 Available online 2 March 2008
a b s t r a c t We propose a fast, explicit numerical method for computing approximations for the immersed boundary problem in which the boundaries that separate the fluid into two regions are stiff. In the numerical computations of such problems, one frequently has to contend with numerical instability, as the stiff immersed boundaries exert large forces on the local fluid. When the boundary forces are treated explicitly, prohibitively small time-steps may be required to maintain numerical stability. On the other hand, when the boundary forces are treated implicitly, the restriction on the time-step size is reduced, but the solution of a large system of coupled non-linear equations may be required. In this work, we develop an efficient method that combines an integral equation approach with the immersed interface method. The present method treats the boundary forces explicitly. To reduce computational costs, the method uses an operator-splitting approach: large time-steps are used to update the non-stiff advection terms, and smaller substeps are used to advance the stiff boundary. At each substep, an integral equation is computed to yield fluid velocity local to the boundary; those velocity values are then used to update the boundary configuration. Fluid variables are computed over the entire domain, using the immersed interface method, only at the end of the large advection time-steps. Numerical results suggest that the present method compares favorably with an implementation of the immersed interface method that employs an explicit time-stepping and no fractional stepping. Ó 2008 Elsevier Ltd. All rights reserved.
1. Introduction Frequently in problems of fluid flow or wave propagation, a boundary between different regions exerts a force on the material, or an interface separates regions of different material properties. Substantial effort has been directed to developing computational techniques for simulating moving boundaries or interfaces within a viscous fluid domain. The most notable example is perhaps the immersed boundary method, developed by Peskin [26] for solving the full incompressible Navier–Stokes equations with moving boundaries. The immersed boundary method was originally developed for studying blood flow through a beating heart [25]. In the immersed boundary method, moving boundaries are represented by means of Lagrangian markers, at which boundary forces are computed. These forces are then transferred to an underlying Cartesian grid via a regularization using discrete (smoothed) delta functions, and the Navier–Stokes equations are solved on the Cartesian grid. In principle, the Lagrangian representation of the immersed boundary allows complicated boundaries to be incorporated without significant difficulty. The immersed boundary methq
This work was supported in part by the National Science Foundation, Grant DMS0715021. E-mail address: [email protected]
0045-7930/$ - see front matter Ó 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.compfluid.2008.02.003
od has since been applied to a wide variety of problems [1,3,7,8]. Because the approximate Dirac delta functions, which transfers the singular boundary forces onto the underlying fluid, typically have OðhÞ support, where h denotes the spatial mesh spacing, the immersed boundary method does not capture the jump discontinuity in solution (e.g., pressure) but rather approximates the solution as a continuous function with a large gradient. In general, the immersed boundary method computes approximations with firstorder spatial accuracy. An alternative approach that captures jumps in solution and its derivatives sharply and that generates approximations with second-order spatial accuracy is the immersed interface method developed by LeVeque and Li [17,18]. The immersed interface method is similar to a method developed earlier by Mayo for solving elliptic problems on an irregular domain [21,22]. Both methods are Cartesian grid methods with second-order spatial accuracy (though higher-order immersed interface methods have recently been developed [12,20]), with the key idea being the incorporation of known jumps in the solution or its derivatives into the finite difference schemes. One difficulty that researchers of immersed boundary problems have to contend with is that boundary forces arising from stiff immersed boundaries frequently impose a severe restriction on timestep size in order to maintain numerical stability. When the
A.T. Layton / Computers & Fluids 38 (2009) 266–272
267
boundary forces are treated explicitly, stiff forces impose a prohibitive restriction on the time-step, thereby increasing the computational costs. Thus, much effort has been invested in developing implicit and semi-implicit versions of the immersed boundary method and related methods, e.g. [9,16,18,23,28]. Nonetheless, the treatment of stiff forces remains a challenge. Because of the coupling between boundary forces and the underlying fluid, an implicit or semi-implicit treatment of the boundary forces requires the solution of a large system of non-linear equations via an iterative procedure, the convergence of which can be a concern. Perhaps owing to that difficulty, the majority (though by no means all) of the implementations of the immersed boundary and interface methods are explicit ones. As a result, when the computations become unstable, one is left to wonder whether the problem is attributable to a programming error, or whether one should try a yet smaller time-step, or in the case of an implicit method, a better initial guess for the iterations. The principal goal of this work is to develop an efficient timestepping approach for advancing solutions of immersed boundary problems for the case where the boundary forces are stiff. The proposed method offers the following advantages:
boundary is considered here, the generalization to multiple, closed immersed boundaries is straightforward. We seek to compute a solution on a bounded computational domain X ¼ C þ X þ Xþ , where Xþ 2 X1þ . The immersed boundary C is assumed to be elastic, so that when it is distorted from its relaxed state by stretching or compressing, it exerts a force on the local fluid. The model equations are the Navier–Stokes equations with forcings: ou ð1Þ þ u ru ¼ rp þ lr2 u þ f; q ot
Boundary forces are treated explicitly, thus avoiding the solution of a large system of coupled non-linear equations. With the forces treated explicitly, a small time-step may be needed, which can lead to large overall computational costs. The proposed method keeps the overall computational costs low by means of an integral equation, which provides a fast way of updating the boundary forces.
for i ¼ 1; 2. L0 is the length of C at its relaxed state; x 2 X, and Xðs; tÞ gives the location of C at time t, parametrized by the arclength s at the unstretched state for 0 6 s 6 L0 ; F i is the force strength at point s; and d is the two-dimensional delta function. Note that d is not the approximate (smoothed) delta function used in the immersed boundary method, and that f is singularly supported along C. We consider an elastic band, where the force strength F is given by
r u ¼ 0;
ð2Þ
where u ¼ ðu; vÞ denotes the fluid velocity; p is the pressure; q and l are the fluid density and viscosity, respectively, assumed to be constant; and f is a force term singularly support along C, with the xand y-components denoted by f1 and f2 , respectively. The forces f1 and f2 are elastic or tension forces that arise from C being stretched or compressed. These singular forces can be written as Z L0 fi ðx; tÞ ¼ F i ðs; tÞdðx Xðs; tÞÞds; ð3Þ 0
Our approach is based on a combination of the integral solution of the unsteady Stokes equations and the immersed interface method. Boundary integrals have been used to remove the stiffness from interfacial flows with surface tension [10]. Indeed, the idea of combining the boundary integral method and immersed interface method has been previously applied in other fluid problems. For example, using that approach, Beale and Layton solved the Stokes problem with discontinuous viscosity [2,14], and Biros et al. computed solutions for the Stokes equations [4] and for the Navier– Stokes equations [5] on irregular domains. In this work, the integral solution of the unsteady Stokes equations is used to approximate fluid velocity local to the immersed boundary. That boundary velocity is then used to update the boundary configuration and forces. Numerical results suggest that the efficiency of the proposed method compares favorably to an implementation of the immersed interface method that employs explicit time-stepping. 2. Numerical method We consider an evolving, closed boundary C that separates the two-dimensional domain X into two subdomains, X1þ (exterior) and X (interior); see Fig. 1. Note that although only one immersed
Ω+
Γ Ω−
Fig. 1. The immersed boundary model configuration.
Fðs; tÞ ¼
o ðTðs; tÞsðs; tÞÞ; os
with the tension Tðs; tÞ given by oX Tðs; tÞ ¼ T 0 1 : os
ð4Þ
ð5Þ
The tension coefficient T 0 describes the elastic properties of the immersed boundary and is assumed to be a constant in this model. The vector tangential to C is given by sðs; tÞ, where sðs; tÞ ¼
oX=os : joX=osj
ð6Þ
Thus, the force density can be computed directly from the location X of the boundary C. Note that at the relaxed state, oX=os ¼ 1 and tension vanishes. In the spatial discretization, let hx and hy denote the spatial grid interval along the x- and y-axis, respectively, and ðihx ; jhy Þ be the ith and jth spatial grid points along the x- and y-axis, respectively, where i; j ¼ 0; 1; 2; . . . ; N. For notational simplicity, let hx ¼ hy h. The position of the immersed boundary at t is represented by a set of boundary markers XðtÞ ðX k ðtÞ; Y k ðtÞÞ, for k ¼ 0; 1;2; .. .; N k , where ðX 0 ðtÞ;Y 0 ðtÞÞ ¼ ðX Nk ðtÞ; Y Nk ðtÞÞ because the boundary is assumed to be a simple closed curve. The kth boundary marker approximates ðXðsk ;tÞ; Yðsk ;tÞÞ, where sk ¼ kL0 =N k ; that is, the boundary markers are chosen such that in an unstretched or relaxed state, they are equally spaced. The present method updates boundary forces by means of the integral solution of the unsteady Stokes equations (see below). To apply that solution to the Navier–Stokes equation, we first discretize the momentum Eq. (1) in time using a semi-implicit approach, with the inertia and force terms treated explicitly. An operator-splitting approach is also applied, where the stiff forces are updated using a smaller time-step, which we refer to as the forcing substep, than the advection or overall time-step used to update the advection terms. Suppose we take N m (small) forcing sub-
268
A.T. Layton / Computers & Fluids 38 (2009) 266–272
steps Dt m per (large) advection time-step Dtn ; thus, Dt n ¼ N m Dt m . Let tm and t n denote the forcing and advection time-levels, respectively, where t m ¼ tn þ mDt m ¼ tn þ mDtn =N m , for m ¼ 1; 2; . . . ; N m . With a first-order time-differencing, the discretized form of Eq. (1) is given by frac1Dtm umþ1 þ mþ1
ru
1 l 1 m fm u un run þ ; rpmþ1 r2 umþ1 ¼ q q q Dtm
¼ 0:
i;k¼1;2
ð7Þ ð8Þ
n
We then select a point x0 2 Xþ , and integrate Eq. (11) over X . Using the divergence theorem to convert the surface integral over X into a boundary integral over C, we obtain X Z Gij ð^ xÞrij ðxÞ Tijk ð^ xÞui ðxÞ nk dx
n
Note that the advection term u ru is evaluated at t n . A simple first-order time-stepping is used, in part, because the operatorsplitting is first-order. Possible extensions to higher-order time discretization are noted in Section 4. Because the boundary forces are treated explicitly, the size of Dtm may be severely restricted to maintain numerical stability. For the proposed method to be efficient, the computation cost in advancing the solution by one Dtm must be sufficiently cheap. To achieve that, we note that to update the boundary forces, one must update the boundary position; to update the boundary position, one needs to update the fluid velocity local to the boundary—but nowhere else. With that observation, we propose an efficient time-stepping method that updates only the boundary velocity during each Dtm using boundary and volume integrals; fluid velocity is then computed (less frequently) in the entire domain at each advection step Dtn .
¼
C
X Z
i¼1;2
X
Gij ð^ xÞbi ðxÞdx;
ð14Þ
Next, we select a point x0 2 C, and integrate Eq. (11) over X [24]: X1 Z uj ðx0 Þ ¼ Gij ð^ xÞrij ðxÞ Tijk ð^ xÞui ðxÞ nk dx pl C i¼1;2 Z 1 Gij ð^ xÞbi ðxÞdx : ð15Þ þ 2pl X To express the integral solution uj ðx0 Þ in terms of the boundary forces f, we multiply Eq. (14) by 1=ðplÞ and take the difference between Eq. (15) and the resulting equation, which yields Z X 1 Z 1 uj ðx0 Þ ¼ Gij ð^ xÞfi ðxÞdx Gij ð^ xÞbi ðxÞdx : ð16Þ pl C 2pl X i¼1;2 Thus, given boundary forces f, boundary velocity, denoted uC , can be obtained by computing the boundary integrals and the volume integrals on the right side of (16). 2.2. Evaluation of boundary and volume integrals
2.1. Integral solution To derive the method, we drop the superscripts in Eq. (7) and let k2 q=lDtm . Using that notation, we rewrite Eq. (7) as ðr2 k2 Þu rp ¼ b;
ð9Þ
where p is a modified pressure that contains a 1=l factor. The left side contains unknown variables at t mþ1 , and the right side b um =Dt m þ un run contains variables at previous time-levels t n and t m . To develop an integral formulation for Eq. (9), which is also known as Brinkman’s equation, we rewrite Eq. (9) in terms of the stress tensor r: r r k2 u ¼ b:
ð10Þ
Let r and u be a solution of the inhomogeneous Eq. (10), and let r and u be a solution of the homogeneous equation with b ¼ 0. Then by Green’s identity, r, r , u, and u satisfy X o X ui rij ui rij ¼ ui bi : ð11Þ oxj i¼1;2 i¼1;2 We identify r and u with the fundamental solution of Brinkman’s equation, denoted T and G, respectively. In two-dimensional free space, Gij and Tijk are given by ! K 1 ðkrÞ 1 Gij ð^ xÞ ¼ dij K 0 ðkrÞ þ kr ðkrÞ2 ! ^ 2K 1 ðkrÞ 2 xi ^ xj ; ð12Þ þ 2 K 0 ðkrÞ þ r kr ðkrÞ2 ! 4K 1 ðkrÞ 4 xÞ ¼ ^ xj dik 1 þ 2K 0 ðkrÞ þ Tijk ð^ kr ðkrÞ2 ! 4K 1 ðkrÞ 4 xi djk Þ krK 1 ðkrÞ þ þ ð^ xk dij þ ^ kr ðkrÞ2 ! ^ 16K 1 ðkrÞ 16 xi ^ xj ^ xk ; ð13Þ 8K 0 ðkrÞ 2krK 1 ðkrÞ þ 2 þ r kr ðkrÞ2 where ^x x0 x, r j^xj, and K 0 and K 1 are modified Bessel functions of the second kind [24].
Normally, the evaluation of a volume integral can be computaR tionally expensive. However, in the case of X Gij bi dx, which is the second term in the summation on the right side of (16), one might note that Gij decays rapidly away from the poles (i.e., as a function of kr, see Eq. (12)). Thus, in practice, the volume integral can be approximated by a narrow band D neighboring C. In particular, note that k scales as 1=ðlDt m Þ, thus Gij decays faster at lower viscosity l and smaller time-step Dtm . For stiff forces, Dtm must be sufficiently small to maintain numerical stability, rendering the decay of Gij extremely fast. Consequently, the band D along which the volume integral is computed can be made narrow without significantly sacrificing accuracy. A challenge in computing both the boundary and volume integrals in (16) is that when k is large, Gij has large derivatives. (Spex, the distance from the pole, cifically, Gij decreases rapidly as ^ increases.) To reduce quadrature errors, the adaptive trapezoidal rule, which employs an increasingly refined grid near the poles, was used. One technical complication in using an adaptive quadrature grid is that the boundary force f is known only at the boundary markers but not at every quadrature nodes in the refined grid. Although spatial interpolation could be used to estimate f values at the quadrature nodes, one notes that, compared to Gij , f likely does not vary significantly from one quadrature node to its neighbor. Thus, in the evaluation of the boundary integral, we allow Gij to vary at every quadrature node (see below), but assume that f is constant within the interval ½ðXk1 þ Xk Þ=2; ðXk þ Xkþ1 Þ=2. This approach avoids the additional cost of interpolations without significantly sacrificing numerical accuracy. A similar approach is used in the approximation of the volume integral. An adaptively refined grid is used to compute the quadrature that approximates the volume integral over D. The values of b, which varies slowly compared to Gij , are known at grid points but not at each quadrature node. Following the approach used in the evaluation of the boundary integral, we assume that b is constant within the h h square centered at each grid point. Another issue here is that away from C, um (the first term in b) is not known. This is because u is updated at tm only along C, thus elsewhere u is known only at previous t n levels. The values of um away from C can be estimated by extrapolating in time using values at tn and
269
A.T. Layton / Computers & Fluids 38 (2009) 266–272
t n1 . Extrapolation may introduce inaccuracy, but because Gij decays rapidly away from C, many of those terms in the integral that require extrapolation, which are located away from C, likely have a small contribution. The most expensive step in the computation of uC is the evaluation of Green’s function using Eq. (12). Because Eq. (12) involves Bessel functions, the evaluation of which is computationally nontrivial, it is not clear that, in a naïve implementation, computing the boundary velocity using the integral Eq. (16) is less expensive than computing fluid velocity over X. To speed up the evaluation of Gij , we place quadrature nodes at fixed distances from the boundary markers; i.e., the adaptively refined grid used in the trapezoidal rule is the same for every time-step throughout the computation. Specifically, the quadrature nodes used in estimating boundary integral are chosen to be points along C that are located at distances of r q from each boundary marker Xk , for q ¼ 1; 2; . . . ; N q . The locations of the quadrature nodes in the computation of the volume integral are chosen analogously. Because the distances r q between the quadrature nodes and the boundary markers remain fixed, the two terms in Eq. (12) involving the Bessel functions, i.e., ! ! K 1 ðkr q Þ 1 1 2K 1 ðkrq Þ 2 ; ; K ðkr Þ þ K 0 ðkrq Þ þ 0 q kr q r 2q kr q ðkrq Þ2 ðkrq Þ2 can be precomputed. The evaluation of the boundary and volume integrals requires OðN q N k Þ and OðN 2q N k Þ work, respectively. On the other hand, if fluid velocity and pressure were computed everywhere to advance the boundary configuration, as in an implementation of the original immersed interface method (i.e., without the fractional stepping and the integral approach), then the computations would take OðN 2 log NÞ work if fast Fourier transforms are used. If one assumes that N q scales as N, then a comparison of the relative efficiency of those two approaches depends, in large part, on the values of N q and N. The efficiency of the two approaches will be compared using a numerical test below. 2.3. The algorithm Once the boundary velocity um C is obtained by evaluating Eq. (16), the immersed boundary position can be updated. As previously noted, the position X of C is determined by the location of a set of boundary markers along C; and at a given time t, C is represented by a cubic spline through those markers [18]. To advance the boundary by one time-step, we move the boundary markers at the same velocity as the local fluid: o Xðs; tÞ ¼ uðXðs; tÞ; tÞ: ot
method because the method is well-known among the immersed boundary community.) The boundary and volume integrals are not used to compute fluid variables everywhere at the advection time-level t n because such computations cost OðN 2q N2 Þ time and are likely more expensive compared to the immersed interface method. The algorithm for advancing the solution from tn to t nþ1 is summarized below: (1) For m ¼ 1; 2; . . . ; N m 1. (a) Compute boundary velocity uðXk ; t m Þ along C using the integral Eq. (16). using (18). (b) Update boundary position Xmþ1 k mþ1 using (3) and (4). (c) Update boundary forces f (2) Express jumps in solution and derivatives across C in terms N nþ1 of boundary forces f m ¼ f . Compute fluid solution unþ1 and pnþ1 over X using the immersed interface method [19].
3. Numerical examples In this section we present numerical results for two examples. Using these examples, we aim to demonstrate the second-order accuracy of the method in space, and to assess its efficiency. All calculations reported below were performed using Fortran programs, which were executed in double precision on a computer with an Intel Core 2 Duo 6400 2.13 GHz processor, 2 MB of cache, and 1 GB of RAM. 3.1. Example 1 The first example was used by Tu and Peskin to test their immersed boundary method [28], and by LeVeque and Li to test their immersed interface method [18]. The initial boundary is an ellipse with major and minor axes set to a ¼ 0:7 and b ¼ 0:5, respectively (see Fig. 2, dashed curve labelled ‘‘Initial”). The unstretched boundary was set to a circle with radius r 0 ¼ 0:5. The tension coefficient T 0 was set to 0.5. The fluid density q was set to 1. At t ¼ 0, the fluid
1
ð17Þ
0.5
Forward Euler is used to advance the boundary as follows: m Xmþ1 ¼ Xm k k þ DtuC k :
ð18Þ
Once a new boundary configuration has been obtained, the boundmþ1 can then be computed using (3) and (4), and a new ary force f forcing step Dt m can be taken. After N m force substeps Dt m , the fluid velocity and pressure can be computed over the entire domain X using the immersed interface method [19]. Jumps in the fluid variables and their derivatives are computed from the boundary forces and their derivatives; those jumps are then incorporated into a finite-difference scheme; for details see [19]. One technical detail that should be noted is that in this implementation, jumps in variables and their derivatives are computed along coordinate lines, as in a previous study [15] and in Mayo’s method [21,22] but unlike the immersed interface method which computes those jumps in the directions normal and tangential to C. (Nonetheless, we refer to the present method as a variation of the immersed interface
Final
0
Initial -0.5
-1
-1
-0.5
0
0.5
1
Fig. 2. Boundary configuration at initial time (dashed line) and at steady state (solid dark line).
270
A.T. Layton / Computers & Fluids 38 (2009) 266–272
was initialized to be at rest. The computational domain was set to ½2:5; 2:5 ½2:5; 2:5. At steady state, the initial ellipse converged to a circle with rapffiffiffiffiffi ffi dius re ¼ ab 0:6124, which is larger than the unstretched boundary but which has the same area as the initial ellipse, owing to the incompressibility of the enclosed fluid. At steady state, fluid velocity vanished everywhere; p attained constant values inside and outside the boundary, with a jump ½p < 0 (because the boundary was initialized to a stretched state and because fluid velocity vanished at steady state). We simulated the motion of the fluid and of the boundary for t 2 ½0; 5 using an overall time-step Dtn ¼ 2 103 , together with a smaller forcing substep of Dtm ¼ 2 104 . For a given Dtm , the xj depends in rate at which kGð^ xÞk1 decreases as a function of j^ large part on the fluid viscosity l. As an example, consider the case x2 > 0. The values of kGð^ xÞk1 for differing values of l in which ^ x1 ¼ ^ and r ¼ j^ xj are shown in Fig. 3. To further compare the boundary layers for differing viscosity values, for example, l ¼ 105 ;104 ; xj > 0:04473, 0.1414, 103 ; 102 , and 101 , kGk1 < 106 for j^ 0.4472, 1.1414, and 4.472, respectively. Thus, the boundary layer becomes increasingly thin as l decreases (see Fig. 3). Although the width of the band D over which the volume integral is computed can be decreased for smaller values of l, the gradient of G increases as l decreases. Thus, to maintain quadrature accuracy, we used approximately the same number of quadrature nodes N q , for a given spatial resolution N, for all l values. To assess the spatial accuracy of the method, we scaled N k , N q , and N simultaneously. Specifically, we let N k ¼ N, and lowered the error tolerance of the adaptive trapezoidal rule by a factor of 4 as N doubled, thereby increasing N q . Volume accuracy results for the case of l ¼ 104 are shown in Table 1 for N ¼ 80, 160, 320, and 640. The volume errors are given by the relative difference in the area bounded by C, between exact solution and the computed solution at approximate steady state. These results exhibit second-order convergence at sufficiently large N values. At N ¼ 80, there was substantial water loss in X , an error that can be attributed to the under-resolved boundary layer neighboring C. Nonetheless, with a spatial grid that is sufficiently refined to resolve the boundary layer, volume is approximately conserved. For example, for N P 160, the relative volume error is 6 0:3%. Analogous results were obtained also for l ¼ 103 and 105 . We also estimated relative errors in the fluid velocity u ðu; vÞ. Because the exact solution is not known for this problem, the approximation computed using a refined spatial discretization of N ¼ 1280 is used as a refer-
Table 1 Convergence results. A ð1Þ, steady-state area enclosed by C N
A ð1Þ
errorA
kerroru k2
80 160 320 640
0.71382151 0.78297285 0.78480984 0.78524711
9.113E2 3.088E3 7.491E4 1.923E4
1.261E1 4.686E2 1.384E2 3.771E3
The relative errors in A ð1Þ A ð1Þ ¼ 0:7853981635. Results convergence.
are computed using the exact solution show approximate second-order spatial
ence solution. The results, given in L2-norm, exhibit approximately second-order accuracy; see Table 1. 3.2. Example 2 The second test involves a more complicated initial boundary configuration, which is given in polar coordinates by r ¼ 0:8þ 0:3 sin 8h. The relaxed boundary is the circle with radius r0 ¼ 0:3. The initial and relaxed boundary configurations are shown in Fig. 4 (dashed and dotted curves, respectively). The larger curvature, compared to the first example, renders this problem significantly stiffer, inasmuch as the boundary tension force increases with the curvature. The density and viscosity coefficients of the fluid were taken to be q ¼ 1 and l ¼ 104 , respectively. The tension coefficient T 0 was set to 0.5. The fluid was initialized to be at rest. The spatial domain was set to ½1:2; 1:2 ½1:2; 1:2. The system was integrated using an overall time-step Dtn of 0.01, and 40 forcing substeps Dt m ¼ 2:5 104 , the largest forcing substep size that can be taken without incurring numerical instability. A spatial grid of N ¼ 640 subintervals was used. Nine quadrature nodes were used in each dimension to approximate the integrals, i.e., N q ¼ 13. With these parameters, we compared the computational costs of the present method to an explicit implementation of the traditional immersed interface method that does not involve fractional stepping. To maintain numerical stability, the immersed interface method requires a time-step Dt of 2:5 104 , which is equal to the forcing substep Dtm of our method. At each time-step, the immersed interface method solved the Navier–Stokes equations (7) and (8), and advanced fluid and boundary variables over the entire computational domain X. For this example, and with the time-steps specified, the present meth-
1.2 2
0.8
|G|
1.5
μ = 10
0.4
−1
μ = 10
1
−2
0
μ = 10
−3
μ = 10
0.5
-0.4
−4
μ = 10
−5
-0.8 0 0
0.0025
0.005
0.0075
0.01
0.0125
-1.2 -1.2
-0.8
-0.4
0
0.4
0.8
1.2
r Fig. 3. kGðrÞk1 as a function of r, for differing l values. Dtm ¼ 2 104 ; ^ x1 ¼ ^ x2 > 0.
Fig. 4. Boundary configuration at initial time (dashed line) and at t ¼ 0:1 (solid line). Dotted line shows relaxed configuration.
A.T. Layton / Computers & Fluids 38 (2009) 266–272
od was approximately 60% faster, measured in CPU time, than the immersed interface method. 4. Discussion In this work, we combined an integral approach with the immersed interface method for an efficient solution of immersed boundary problems with stiff forcings. The method integrates the boundary forces explicitly, an approach that requires small timesteps, compared to an implicit treatment of stiff force terms, but that avoids the need to solve a system of coupled non-linear equations. To reduce computational costs, fluid velocities values are computed only along the immersed boundary using an integral solution of the unsteady Stokes equations. Numerical test results suggest that the efficiency of the proposed method compares favorably to an explicit non-fractional stepping implementation of the original immersed interface method. A direct comparison between the present method (which we refer to below as the ‘‘integral immersed interface method”) and the immersed interface method is difficult to make for a general case, inasmuch as the relative merits of the two methods depend on a number of factors including spatial resolution, stiffness of the problem, accuracy requirement, etc. Moreover, using CPU time as a measure of computational costs has the drawback that metric depends not only on the problem at hand, but also on the particular implementation and the computer architecture. Thus, below we analyze the computational costs of the two methods as functions of spatial resolution. As previously noted, the evaluation of the volume integral at a forcing substep of the integral immersed interface method requires OðN 2q N k Þ work. An immersed interface method using fast Fourier transforms requires OðN 2 log NÞ work for each time-step. If both N q and N k scales as N, then technically the computational cost of the integral immersed interface method is OðN 3 Þ, which should be higher than the original immersed interface method at sufficiently large N’s. Nonetheless, in practice, N q N; thus, at moderate N values, the integral immersed interface method can be faster, as we have shown in Section 3.2. Furthermore, for large N q , specialized fast integration methods for nearly singular functions, e.g. [13,27], may be used. Compare to an implementation of the immersed interface method that uses fast Fourier transforms, the integral immersed interface method may be better suited for parallelization. In two- or three-dimensional computations, Fast Fourier transforms require transposing data, a step which gives rise to global communication among processors. In contrast, computing the integrals in the forcing substeps of the integral immersed interface method requires mostly data local to a processor. Because it requires fewer global inter-processor communications, a parallel implementation of the integral immersed interface method likely has a lower overall computational time, compared to a parallel implementation of the immersed interface method using fast Fourier transforms. The methodology of the present method can also be applied to the immersed boundary method, where the boundary forces are not singularly supported along C but are instead spread onto the fluid grid using approximate delta functions [26]. Assuming that the approximate delta function is compactly supported over D Dþ [ D , where Dþ ½C; C þ dn and D ½C dn; C, then the boundary integral takes the form Z Z X 1 1 Gij ð^ xÞfi ðxÞdx Gij ð^ xÞbi ðxÞdx ð19Þ uj ðx0 Þ ¼ 2pl D 2pl X i¼1;2 where fi above is the smoothed forcing term. A limitation of the integral immersed interface method might be that it requires Green’s function of the unsteady Stokes equa-
271
tion, and that the form of Green’s function depends on the number of spatial dimensions, the shape of the fluid domain, boundary conditions, etc. Thus, Green’s function and, as a result, the integral solution may not be available for arbitrary domain and boundary conditions. Nonetheless, an approximate integral solution may be obtained using the Green’s function that corresponds to a different domain or boundary conditions. The impact of using an approximate Green’s function on numerical stability and accuracy will need to be accessed. The present method is second-order in space, but only first-order in time. We opted for a first-order time-stepping scheme because of the simplicity of a first-order operator-splitting, and because of the likelihood (not verified) that order reduction, which is attributable to the stiffness of the boundary, would restrict solution accuracy to first order regardless. In addition, without incorporating temporal correction terms as in [29], a second-order time-stepping method such as that used in an implementation of the immersed interface method [19] would have yielded approximations with only first-order accuracy in time for problems with moving interfaces. Nonetheless, a version of the method that is second or higher order in time may be desirable. Strang splitting can be used to develop a second-order time-stepping method with operator-splitting. For third- or higher-order methods with operator-splitting, the multi-implicit spectral deferred correction (SDC) approach may be used to generate a high-order solution using a low-order time integration method and a low-order splitting, using an iterative procedure [6]. To attain high-order temporal accuracy for a stiff problem, we note that a version of the SDC method that uses the GMRES Krylov subspace method has been shown in the fully-explicit and fully-implicit cases to eliminate the order reduction phenomenon previously observed for stiff ODE systems [11]. That approach, together with appropriate temporal correction terms [29], may be used to generate high-order approximations in an immersed boundary problem with stiff forces. References [1] Arthurs KM, Moore LC, Peskin CS, Pitman EB, Layton HE. Modeling arteriolar flow and mass transport using the immersed boundary method. J Comput Phys 1998;147:402–40. [2] Beale T, Layton AT. On the accuracy of finite difference methods for elliptic problems with interfaces. Commun Appl Math Comput Sci 2006. [3] Beyer RP. A computational model of the cochlea using the immersed boundary method. J Comp Phys 1992;98:145–62. [4] Biros G, Ying L, Zorin D. An embedded boundary integral solvers for the Stokes equations. J Comput Phys 2004;193:317–48. [5] Biros G, Ying L, Zorin D. An embedded boundary integral solvers for the unsteady incompressible Navier–Stokes equations. Technical report, Ms-CIS04-42, 2004. [6] Bourlioux A, Layton AT, Minion ML. Higher-order multi-implicit spectral deferred correction methods for problems of reacting flow. J Comput Phys 2003;189:651–75. [7] Dillon R, Fauci L, Gaver III D. A microscale model of bacterial swimming, chemotaxis, and substrate transport. J Theor Biol 1995;177:325–40. [8] Fauci L, McDonald A. Sperm mobility in the presence of boundaries. Bull Math Biol 1995;57:679–99. [9] Fauci LJ, Folgelson AL. Truncated newton methods and the modeling of complex immersed elastic structures. Commun Pure Appl Math 1993;66:787–818. [10] Hou TY, Lowengrub JS, Shelley MJ. Removing the stiffness from interfacial flows with surface tension. J Comput Phys 1994;114:312–38. [11] Huang J-F, Jia J, Minion ML. Accelerating the convergence of spectral deferred correction methods. J Comput Phys 2006;214(2):633–56. [12] Ito K, Kyei Y, Li Z. Higher-order Cartesian grid based finite difference schemes for elliptic equations on irregular domains. SIAM J Sci Comput 2005;27(1):346–67. [13] Kapur S, Rokhlin V. High-order corrected trapezoidal quadrature rules for singular functions. SIAM J Numer Anal 1997;34(3):1331–56. [14] Layton AT. An efficient numerical method for the two-fluid stokes equations with a moving boundary. Comput Methods Appl Mech Eng [in press]. [15] Layton AT. Modeling water transport across elastic boundaries using an explicit jump method. SIAM J Sci Comput 2007;28(6):2189–207.
272
A.T. Layton / Computers & Fluids 38 (2009) 266–272
[16] Lee L, LeVeque RJ. An immersed interface method for the incompressible Navier–Stokes equations. SIAM J Sci Comp 2003;25:832–56. [17] LeVeque RJ, Li Z. The immersed interface method for elliptic equations with discontinuous coefficients and singular sources. SIAM J Numer Anal 1994;31:1019–44. [18] LeVeque RJ, Li Z. Immersed interface methods for Stokes flow with elastic boundaries or surface tension. SIAM J Sci Comput 1997;18(3):709–35. [19] Li Z, Lai M-C. The immersed interface method for the Navier–Stokes equations with singular forces. J Comput Phys 2001;171:822–42. [20] Linnick MN, Fasel HF. A high-order immersed interface method for simulating unsteady incompressible flows on irregular domains. J Comput Phys 2005;204:157–92. [21] Mayo A. The fast solution of Poisson’s and the biharmonic equations on irregular regions. SIAM J Numer Anal 1984;21:285–99. [22] Mayo A, Greenbaum A. Fast parallel iterative solution of Poisson’s and the biharmonic equations on irregular regions. SIAM J Sci Comput 1982;13:101–18.
[23] Mayo A, Peskin CS. An implicit numerical method for fluid dynamics problems with immersed elastic boundaries. In: Cheer AY, von Dam CP, editors. Fluid Dynamics in Biology. Providence, RI: AMS; 1993. p. 261–78. [24] Occhialini JM, Muldowney GP, Higdon JJL. Boundary integral/spectral element approaches to the Navier–Stokes equations. Int J Numer Methods Fluids 1992;15:1316–81. [25] Peskin CS. Numerical analysis of blood flow in the heart. J Comput Phys 1977;25:220–52. [26] Peskin CS. The immersed boundary method. Acta Numer 2002:1–39. [27] Strain J. Locally-corrected multidimensional quadrature rules for singular functions. SIAM J Sci Comput 1995;6(4):992–1017. [28] Tu C, Peskin CS. Stability and instability in the computation of flows with moving immersed boundaries: a comparison of three methods. SIAM J Sci Statist Comput 1992;13:1361–76. [29] Xu S, Wang ZJ. Systematic derivation of jump conditions for the immersed interface method in three-dimensional flow simulation. SIAM J Sci Comp 2006;27:1948–80.