DOI 10.1007/s10915-014-9933-y
An Unconditionally Energy Stable Penalty Immersed Boundary Method for Simulating the Dynamics of an Inextensible Interface Interacting with a Solid Particle
Po-Wen Hsieh · Ming-Chih Lai · Suh-Yuh Yang · Cheng-Shu You
Received: 11 August 2013 / Revised: 29 May 2014 / Accepted: 9 October 2014 / Published online: 21 October 2014
© Springer Science+Business Media New York 2014
Abstract In this paper, a novel penalty method based on the immersed boundary formulation is proposed for simulating the transient Stokes flow with an inextensible interface enclosing a suspended solid particle. The main idea of this approach relies on the penalty techniques by modifying the constitutive equation of Stokes flow to weaken the incompressibility condition, relating the surface divergence to the elastic tensionσ to relax the interface’s inextensibility, and connecting the particle surface-velocity with the particle surface force F to regularize the particle’s rigid motion. The advantage of these regularized governing equations is that when they are discretized by the standard centered difference scheme on a staggered grid, the resulting linear system can easily be reduced by eliminating the unknowns ph, σh and Fh directly, so that we just need to solve a smaller linear system of the velocity approximation uh. This advantage is preserved and even enhanced when such approach is applied to the
Po-Wen Hsieh was partially supported by the National Science Council of Taiwan under the Grant NSC 102-2115-M-033-007-MY2.
Ming-Chih Lai was partially supported by the National Science Council of Taiwan under the Grant NSC 101-2115-M-009-014-MY3.
Suh-Yuh Yang was partially supported by the National Science Council of Taiwan under the Grant NSC 101-2115-M-008-008-MY2.
P.-W. Hsieh
Department of Applied Mathematics, Chung Yuan Christian University, Jhongli City 32023, Taoyuan County, Taiwan
e-mail: pwhsieh0209@gmail.com M.-C. Lai
Center of Mathematical Modeling and Scientific Computing, Department of Applied Mathematics, National Chiao Tung University, Hsinchu 30010, Taiwan
e-mail: mclai@math.nctu.edu.tw S.-Y. Yang (
B
)· C.-S. YouDepartment of Mathematics, National Central University, Jhongli City 32001, Taoyuan County, Taiwan e-mail: syyang@math.ncu.edu.tw
C.-S. You
e-mail: csdyou@gmail.com
transient Stokes flow with multiple compound vesicles. Moreover, this smaller linear system is symmetric and negative-definite, which enables us to use efficient linear solvers. Another important feature of the proposed method is that the discretization scheme is unconditionally stable in the sense that an appropriately defined energy functional associated with the discrete system is decreasing and hence bounded in time. We numerically test the accuracy and stability of the immersed boundary discretization scheme. The tank-treading and tumbling motions of inextensible interface with a suspended solid particle in the simple shear flow will be studied extensively. The simulation of the motion of multiple compound vesicles will be performed as well. Numerical results illustrate the superior performance of the proposed penalty method.
Keywords Immersed boundary method· Penalty method · Stokes flow · Inextensible interface· Solid particle · Stability
Mathematics Subject Classification 65M06· 65M12 · 76D07 · 76M20
1 Introduction
In recent years, the study of vesicle dynamics has been the focus of intense research. Vesi- cles provide the simplest model system for simulating complex behavior of biomembranes suspended in the biological fluids. For example, vesicle can be considered as a simplified model for red blood cell in blood flow, since they share many similar characteristic behavior, such as the tank-treading, tumbling and vacillating-breathing (swinging) motions [15,17].
Basically, such a system consists of two fluids separated by an elastic membrane having the inextensible property, which may be deformed due to the interaction with the biological fluids. Besides, vesicles can also be considered as promising drug carriers for the delivery at specific locations in the organisms [1]. These explain the increasing interest and impor- tance of understanding the vesicle dynamics. In essence, the dynamics of vesicles can be determined by their boundary rigidity, inextensibility, and the hydrodynamical forces.
On the other hand, the immersed boundary (IB) method was introduced by Peskin in the early seventies to model blood flow in the heart and through heart valves. Currently, it has evolved into a simple but powerful method for formulating the coupled equations of motion of a viscous, incompressible fluid with one or more elastic surfaces immersed in the fluid; see [18] and many references therein. In the IB method, an Eulerian description is used for the fluid dynamics, while a Lagrangian form is used for each immersed object. The key idea of the IB method is replacing each suspended object by a suitable contribution to a force density term in the fluid dynamics equations. This allows a single set of fluid dynamics equations to hold in the entire spatial domain without any internal boundary conditions. Moreover, without generating an interface-fitting grid for both exterior and interior regions of each surface at each time step, instead, the IB discretization schemes can be implemented by employing a uniform Cartesian grid over the entire domain and the immersed boundaries are discretized by a set of points that are not constrained to lie on the grid. In the meantime, the Eulerian and Lagrangian variables are linked by interaction equations that involve a smoothed approximation to the Dirac delta function.
Recently, Lai et al. [11] developed a fractional step IB method for Stokes flow with an inextensible interface enclosing a suspended solid particle. In addition to solving for the fluid variables, velocity and pressure, the proposed system in [11] involves finding an extra unknown, elastic tensionσ , such that the surface divergence of the velocity is vanished along
the interface, and an extra unknown, particle surface force F, such that the velocity satisfies the no-slip boundary condition along the particle surface. The interface moves with local fluid velocity, while the enclosed particle undergoes a rigid body motion. Finally, the force-free and torque-free conditions along the particle surface are imposed to close the system. They showed that the nullity of the linear algebraic system arising from the centered discretizations of the IB equations over a staggered grid is nonzero, and thus the existence of a solution is guaranteed. They then applied the idea of the fractional step method developed in [20] to solve the resultant linear system of equations. Although the proposed system in [11] seems not to be a satisfactory model for certain biological cells such as the human leukocytes, a lipid bilayer membrane enclosing a fluid with a core, it can be viewed as a heuristic model for developing efficient numerical schemes for analyzing the dynamics of compound vesicles.
In this paper, still based on the formulation of [11], we will propose a novel penalty IB method for the transient Stokes flow problem with an inextensible interface enclosing a suspended solid particle. The main idea of the proposed approach relies on the three penalty techniques. First, we modify the constitutive equation of Stokes flow (cf. [5]) to weaken the incompressibility condition∇ · u = 0 by ∇ · u + εp = 0 with a small penalty parameterε > 0. This technique has been well studied in the finite element computation for the incompressible viscous flow problems to circumvent the cumbersome constraint of incompressibility [21]. Second, we assume that the elastic tensionσ is given in a specific form to follow the usual Hooke’s law for an elastic body [8,9]. This enables us to relate the surface divergence to the elastic tension such that the interface is relaxed to nearly inextensible.
Finally, we connect the particle surface-velocity with the particle surface force F through the use of a virtual particle to regularize the particle’s rigid motion [9]. The advantage of the regularized governing equations is that when they are further discretized by the standard centered difference scheme on a staggered grid, the resulting linear system can easily be reduced by eliminating the unknowns ph, σh and Fh directly so that we just need to solve a smaller linear system of the velocity approximation uh. Moreover, the linear system of uh is symmetric and negative-definite. This enables us to use efficient solvers, such as the preconditioned conjugate gradient and the algebraic multigrid methods, to find solution of the linear system. We emphasize that this advantage will be preserved and even enhanced when such penalty IB approach is applied to the transient Stokes flow with multiple compound vesicles. Indeed, for multiple vesicle problems, the resulting large linear system of equations is still tractable to block elimination, no matter how many compound vesicles are involved in the fluid, and we only need to solve a reduced system of the velocity. Also, the present approach can be easily extended to Navier–Stokes flow by treating the nonlinear advection terms explicitly in the time integration.
It is well known that the discretizations of the IB equations suffer from a severe time step restriction for maintaining the stability and this time step restriction is typically much more stringent than the one that would be imposed from using explicit differencing techniques; see [2,3,13,18,19]. In recent years, considerable effort has been devoted to developing implicit and semi-implicit schemes to alleviate this severe restriction [3,13,14,22]. The instability of these schemes is known not to be a problem related to the advection terms in the incompress- ible Navier–Stokes equations. Also, it has been pointed out in [16] that discretization schemes need not be fully-implicit in order to achieve unconditional stability. One of important fea- tures of the present penalty IB approach is that our semi-implicit discretization scheme for the penalty IB equations is unconditionally stable in the sense that when an appropriately energy functional associated with the discrete system is defined, we can prove that the energy is decreasing, and hence bounded, in time. This result reflects the assertion of [16] in some measure.
Fig. 1 A diagram of an inextensible interfaceΓ enclosing a suspended solid particle P
Γ Ω
P
Since our penalty IB approach relaxes the interface to be nearly inextensible, it is important to us to measure how extent the inextensibility is conserved as time step advances. In this paper, we will prove that the difference of local stretching factors for two successive time steps is still of first order in time step sizeΔt, i.e., O(Δt). This result is similar to that of [11]. In addition, we find that though the difference of local stretching factors is of first order in time step size, the length of the elastic interface may not be always increasing as much as expected when time step advances. In other words, we could keep the inextensibility of the interface very well. This assertion is very different with [11] and the numerical results reported in Sect.5will support this observation.
We will test the accuracy and stability of the semi-implicit discretization scheme through a number of numerical experiments. The first one, which is quoted from [11], is a Stokes problem without the effects of the inextensible interface and the suspended solid particle. The exact solution of this problem is given and so we can easily compute the errors and estimate the convergence rates. Then we perform a series of numerical simulations for inextensible interface with a suspended solid particle in the simple shear flow. In particular, the tank- treading and tumbling motions will be studied extensively. The simulation of the motion of multiple compound vesicles will be performed as well. Numerical results illustrate the superior performance of the proposed penalty IB approach.
The remainder of this paper is organized as follows. In Sect.2, we regularize the governing equations by the penalty techniques. In Sect.3, we discretize the IB equations using the standard centered difference scheme on a staggered grid, and the unique solvability of the resulting linear system is also studied. Some advantageous properties of the penalty IB method, including the unconditional stability and the first-order error of inextensibility for two successive time steps, will be derived in Sect.4. Numerical results are presented in Sect.5. Some concluding remarks are given in Sect.6.
2 Regularization of the Governing Equations of Motion
In this paper, we consider a moving immersed inextensible interfaceΓ enclosing a sus- pended solid particle P in a two-dimensional fluid domainΩ; see Fig.1. We assume that the fluids inside and outside of the interface are the same and governed by the unsteady incompressible Stokes equations, and the interface is massless and the particle’s grav- ity is neglected. Let the elastic interfaceΓ and the particle surface ∂ P be parameterized by X(s, t) = (X1(s, t), X2(s, t)) and Y(α, t) = (Y1(α, t), Y2(α, t)), respectively, where s∈ [0, L1] and α ∈ [0, L2] are the corresponding Lagrangian parameters. The governing equations of motion in dimensionless form can be formulated as follows (cf. [11]): for t> 0,
∂u
∂t + ∇ p = μΔu +
Γ
∂(σ τ)(s, t)
∂s δ(x − X(s, t))ds +
∂ PF(α, t)δ(x − Y(α, t))dα in Ω, (1)
∇ · u = 0 in Ω, (2)
∇s· U =
∂U
∂s · τ
∂ X
∂s
−1= 0 on Γ, (3)
∂ X
∂t (s, t) = U(s, t) =
Ωu(x, t)δ(x − X(s, t))dx on Γ, (4)
∂Y
∂t(α, t) = V(α, t) =
Ωu(x, t)δ(x − Y(α, t))dx = Vc+ ωr on ∂ P, (5)
∂ P F(α, t)dα = 0, (6)
∂ P F(α, t) · r(α, t)dα = 0, (7)
whereμ > 0 is the kinematic viscosity, u = (u, v) and p are the velocity field and the pressure, respectively,σ is the elastic tension and F = (F1, F2) is the particle surface force, Vc(t) is the translational velocity of the center of particle and ω(t) is the angular velocity component of the particle,τ is the unit tangential vector on Γ, r(α, t) = (−(Y2(α, t) − Y2c(t)), Y1(α, t) − Y1c(t)) is a tangential vector on ∂ P, (Y1c, Y2c) is the center of mass of P, δ(x) := δ(x)δ(y) is a two-dimensional Dirac delta function, and ∇s· denotes the surface divergence operator.
We note that Eqs. (1) and (2) are the unsteady incompressible Stokes equations with singular force terms arising from the elastic interface and the particle surface. Equations (3) and (4) represent that the interface is inextensible and moves along with the local fluid velocity so that the velocity U is the interpolation of the fluid velocity at the interface. Equation (5) describes the particle surface velocity which is consisted of the translational velocity Vcand the angular velocityω making the particle moving like a rigid body. Finally, the system of equations will be closed by coupling (6) and (7), which mean the force-free and torque-free conditions for the rigid body motion.
We now introduce the penalty techniques [5,8,9,21] to regularize Eqs. (2), (3) and (5).
First, by modifying the constitutive equation of Stokes flow (cf. [5]), we are able to weaken the incompressibility condition∇ · u = 0 in Ω by
∇ · u + εp = 0 in Ω, (8)
whereε > 0 is a small penalty parameter. This penalty technique has been well studied in the finite element computation for the incompressible Stokes equations. For example, it has been shown that the solution(uε, pε) of the penalty approach will converge to the Stokes solution(u, p) in some suitable norms as ε → 0+. More precisely, we have
uε→ u in the H1(Ω) norm and pε→ p in the L2(Ω) norm as ε → 0+.
For the theoretical analysis, we refer the reader to [21] for more details. In the practical computation, an appropriate choice of the parameterε has been suggested by Hughes et al.
[5] thatε = 1/(cμ), where μ is the dynamic viscosity and the constant c is taken as c = 107 for Stokes flow. They pointed out that the choice of c seems to be problem independent. This is consistent with our numerical experience.
Fig. 2 The reacting force F pulls Y back to target position Y
Second, based on the Hooke’s law, we assume that the elastic tensionσ is given in the form
σ (s, t) = σ0
∂ X
∂s(s, t)
−
∂ X
∂s (s, 0)
, (9)
with a sufficiently large elastic coefficientσ0. Sinceσ0is large, the length of the tangential vector∂ X(s, t)/∂s will be always close to the initial length, i.e.,
∂ X
∂s (s, t)
≈
∂ X
∂s (s, 0)
for all s ∈ [0, L1] and t > 0.
In other words, the interface is no longer exactly inextensible. Instead, the elastic interface is allowed to be nearly inextensible. Indeed, from [12], we have
∂
∂t
∂ X
∂s
= (∇s· U)
∂ X
∂s
. (10)
Combining (9) with (10), the surface divergence free Eq. (3) should be replaced by
∂U
∂s · τ
∂ X
∂s
−1= ∇s· U = 1 σ0
∂σ
∂t
∂ X
∂s
−1 onΓ. (11) In other words, we obtain
∂σ
∂t = σ0
∂U
∂s · τ
onΓ. (12)
To sum up, assumption (9) relates the surface divergence with the elastic tensionσ , as that given in (11), such that the interface is allowed to be nearly inextensible.
Finally, following the Hooke’s law again, we assume that
F= k0(Y − Y), (13)
with a sufficiently large stiffness constant k0, where Y is the parametrization of the virtual solid particle P and Y is assumed to have the same parametric variableα with Y in the same interval[0, L2]. The significance of this assumption can be interpreted as follows. When the flow hits the solid particle, we image that the particle surface may be deformed a little bit to Y but at the same moment, the reacting force F pulls it back to Y , a real position of the particle surface; see Fig.2. Since k0is large, Y(α, t) ≈ Y(α, t) for all α ∈ [0, L2] and t > 0.
Under assumption (13) with the equation of rigid motion of Y ,
∂Y
∂t = V = Vc+ ωr, (14)
we have, combining with (5),
∂ F
∂t = k0
∂Y
∂t −∂Y
∂t
= k0(V − V) = k0(Vc+ ωr − V) on ∂ P. (15) Consequently, we connect the particle surface-velocity V with the particle surface force F by (15). This enables us to regularize the particle’s rigid motion, see (19) and (23) below.
In summary, we obtain the following system of equations, which models the dynamics of a nearly inextensible vesicle enclosing a suspended solid particle in Stokes flow: for t> 0,
∂u
∂t + ∇ p = μΔu +
Γ
∂(σ τ)(s, t)
∂s δ(x − X(s, t))ds +
∂ P F(α, t)δ(x − Y(α, t))dα in Ω, (16)
∇ · u + εp = 0 in Ω, (17)
∂σ
∂t = σ0
∂U
∂s · τ
onΓ, (18)
∂ F
∂t = k0(Vc+ ωr − V) on ∂ P, (19)
∂ P F(α, t)dα = 0, (20)
∂ P F(α, t) · r(α, t)dα = 0, (21)
∂ X
∂t (s, t) = U(s, t) =
Ωu(x, t)δ(x − X(s, t))dx on Γ, (22)
V(α, t) =
Ωu(x, t)δ
x− (Y(α, t) − 1
k0F(α, t))
d x on∂ P, (23) whereε is small but σ0and k0are large enough, and∂ P denotes the surface of the virtual solid particle P, which is parameterized by Y(α, t) = (Y1(α, t), Y2(α, t)) = Y(α, t) − (1/k0)F(α, t) for α ∈ [0, L2].
3 Discretizations of the Immersed Boundary Equations
In this section, we are going to apply the second-order centered difference scheme to discretize the immersed boundary Eqs. (16)–(23) to reach a linear algebraic system. For simplicity, we assume that the computational domain is a rectangular regionΩ = [a, b]×[c, d] and that the fluid variables are defined on the staggered marker-and-cell (MAC) grids [4]. As shown in Fig.3, we define the pressure on the grid points(xi, yj) = (a+(i −1/2)Δx, c+( j −1/2)Δy) for 1 ≤ i ≤ mx and 1≤ j ≤ my, while the velocity components u andv are defined at (xi−1/2, yj) = (a + (i − 1)Δx, c + ( j − 1/2)Δy) for 1 ≤ i ≤ mx+ 1 and 1 ≤ j ≤ myand (xi, yj−1/2) = (a + (i − 1/2)Δx, c + ( j − 1)Δy) for 1 ≤ i ≤ mx and 1≤ j ≤ my+ 1, respectively. Here, we use a uniform mesh with mesh size h= Δx = Δy.
For the immersed nearly inextensible elastic interface X at a given time, we use the Lagrangian markers Xk= X(sk), where sk= kΔs for k = 0, 1, . . . , Me, andΔs is roughly chosen as a half of fluid mesh size h. Since the elastic interface is closed, we have X0 = XMe. We remark that these points will be used in the computation of the values of the discrete delta functions. The elastic tension is defined on the “half-integer” points given by
Fig. 3 A schematic diagram of the computational domainΩ with staggered grid, where the unknowns u, v and p are approximated at the grid points marked by cross, open circle and filled circle, respectively
sk−1/2= (k − 1/2)Δs, so we denote it by σk−1/2. Similarly, we use the Lagrangian markers Yk= Y(αk) on particle surface with αk= kΔα, k = 0, 1, . . . , Mp, andΔα is also roughly chosen as a half of fluid mesh size h.
In what follows, the discrete spatial operators∇h, Δh, and∇sh· denote the standard second- order centered difference approximations to the gradient, Laplacian and surface divergence operators, respectively. For any functionψ(s) defined on the interface Γ , we approximate the partial derivative∂ψ/∂s by the centered difference scheme as
Dsψ(s) = ψ(s + Δs/2) − ψ(s − Δs/2)
Δs .
Thus, the interface stretching factor|∂ X/∂s| and the unit tangent τ can be approximated by|DsX| and DsX/|DsX|, which are defined at the half-integer points. We denote them by
|DsX|k−1/2andτk−1/2, respectively.
LetΔt be the time step size, and the superscript of the variables denote the time step index. At the beginning of each time step n, the elastic tensionσkn−1/2, the particle surface force Fnk, the interface configuration Xnk, the particle surface Ynk= (Y1kn, Y2kn), and its center Ync = (Y1cn, Y2cn) are all given. Despite the face that the problem is non-linear, here we propose a linearly semi-implicit difference scheme for the system of Eqs. (16)–(23). The time step can be advanced as follows:
−un+1− un
Δt + μΔhun+1− ∇hpn+1+
Me
k=1
Ds(σn+1τn)kδh(x − Xnk)Δs
+
Mp
k=1
Fn+1k δh(x − Ynk)Δα = 0, (24)
∇h· un+1+ εpn+1= 0, (25)
σk−1/2n+1 − σk−1/2n
Δt − σ0
Un+1k − Un+1k−1
Δs · τnk−1/2= 0, (26)
Fn+1k − Fnk Δt − k0
Vnc+1+ ωn+1
−(Y2kn − Y2cn) Y1kn − Y1cn
− Vn+1k
= 0, (27)
Mp
k=1
Fnk+1Δα = 0, (28)
Mp
k=1
F1kn+1(Y2kn − Y2cn) − F2kn+1(Y1kn − Y1cl )
Δα = 0, (29)
Xn+1k − Xnk
Δt = Unk+1=
x
un+1(x)δh
x− Xnk
h2, (30)
Vn+1k =
x
un+1(x)δh
x− Ynk+ 1 k0
Fnk
h2. (31)
We remark that since the stiffness constant k0 is large, in practical computation, we may useδh(x − Ynk) to replace the term δh(x − Ynk+k10Fnk) in (31). For the computation of the discrete delta functionδh(x), we adopt the following smoother version developed in [25],
δh(x, y) = 1 h2φx
h
φy h
, (32)
φ(ξ) =
⎧⎪
⎪⎪
⎪⎪
⎪⎪
⎨
⎪⎪
⎪⎪
⎪⎪
⎪⎩
3
8+32π −ξ42, |ξ| < 0.5,
1
4+1−|ξ|8
−2 + 8|ξ| − 4ξ2−18sin−1√
2(|ξ| − 1)
, 0.5 ≤ |ξ| ≤ 1.5,
17
16−64π −3|ξ|4 +ξ82 +|ξ|−216
−14 + 16|ξ| − 4ξ2 +161 sin−1√
2(|ξ| − 2)
, 1.5 ≤ |ξ| ≤ 2.5,
0, 2.5 ≤ |ξ|.
(33)
Now, using the staggered grid for the fluid variables, the matrix obtained by the discrete divergence of the fluid velocity can be written as the transpose of the discrete gradient of the pressure. The matrix produced by the discrete surface divergence of the velocity can be written as the transpose of the matrix obtained by the discrete spreading operator of the tension with some suitable rescaling. Similarly, the matrix obtained by the discrete spreading operator of the force arising from the suspended solid particle boundary can be written as the transpose of the matrix obtained by the discrete interpolating operator of velocity with some suitable rescaling; see [11] and also the discrete skew-adjoint and adjoint properties given in Sect.4below. Thus, the resulting linear system assembled from (24)–(31) is symmetric and can be written as
⎡
⎢⎢
⎢⎢
⎣
A B S E 0
B εI 0 0 0
S 0 ε1I 0 0
E 0 0 ε2I R
0 0 0 R 0
⎤
⎥⎥
⎥⎥
⎦
⎡
⎢⎢
⎢⎢
⎣ unh+1 phn+1 σhn+1 Fn+1h Θn+1h
⎤
⎥⎥
⎥⎥
⎦=
⎡
⎢⎢
⎢⎢
⎣ c1 c2 mσ mF
0
⎤
⎥⎥
⎥⎥
⎦, (34)
whereε1 = Δs/(σ0h2Δt), ε2= Δα/(k0h2Δt), ψhn+1denotes the unknown vector value of ψ at the corresponding grid points, and Θ := (V1c, V2c, ω). The submatrices A, B, S, E, and R respectively represent the discrete Laplacian-like operator,μΔh− (1/Δt)I , the discrete gradient operator,−∇h, the discrete spreading operator on tension, the discrete spreading operator on the surface force, and the discrete rigid body motion equation. In the right-hand side of (34), c1comes from the boundary condition and unh, c2consists only of the velocity boundary conditions since the pressure is evaluated at each cell center, mσ and mF are depending onσhnand Fnh, respectively.
In general, (34) is a large linear system and it is cost-ineffective to solve it directly even if it is a sparse one. We propose an alternative way to efficiently solve the system as follows. From (34), we obviously have phn+1= (1/ε)(c2− Bunh+1) and σhn+1= (1/ε1)(mσ− Sunh+1).
In addition, RFnh+1= 0 and Eunh+1+ ε2Fnh+1+ RΘnh+1= mFeasily imply
REun+1h + RRΘn+1h = RmF. (35) We will show the invertibility of the 3× 3 matrix RR later and after that we obtain
Θnh+1= (RR)−1R
mF− Eunh+1
, (36)
and then
Fn+1h = 1 ε2
(I − R(RR)−1R)mF− (I − R(RR)−1R)Eun+1h
. (37)
Therefore, the linear system (34) can be reduced into a smaller one:
Aunh+1= b, (38)
where the matrix A and the right-hand side vector b are defined as
A:= A −1
εB B− 1 ε1
S S− 1
ε2E(I − R(RR)−1R)E, (39) b:= c1−1
εBc2− 1 ε1
Smσ− 1 ε2
E(I − R(RR)−1R)mF. (40) We now show the invertibility of the 3× 3 matrix RR. First, the matrix R assembled from (27) can be written as
R= −Δα h2
⎡
⎢⎢
⎢⎢
⎢⎢
⎢⎢
⎢⎣
1 0 −(Y21n − Y2cn) ... ... ...
1 0−(Y2Mn p− Y2cn) 0 1 (Y11n − Y1cn)
... ... ...
0 1 (Y1Mn p− Y1cn)
⎤
⎥⎥
⎥⎥
⎥⎥
⎥⎥
⎥⎦ .
To simplify the notation, we define ak= Y1kn − Y1cn and bk= Y2kn − Y2cn. Then we have
h4 Δα2
3
det(RR) = det
⎡
⎢⎢
⎢⎢
⎢⎢
⎢⎢
⎣
Mp 0 −Mp
k=1bk
0 Mp
Mp
k=1ak
−Mp
k=1bk Mp
k=1ak Mp
k=1ak2+Mp
k=1b2k
⎤
⎥⎥
⎥⎥
⎥⎥
⎥⎥
⎦
= Mp
⎧⎪
⎨
⎪⎩Mp Mp
k=1
ak2+ Mp Mp
k=1
b2k−
⎛
⎝
Mp
k=1
ak
⎞
⎠
2
−
⎛
⎝
Mp
k=1
bk
⎞
⎠
2⎫
⎪⎬
⎪⎭ (by Cauchy–Schwarz inequality)
≥ Mp
⎧⎨
⎩Mp
Mp
k=1
ak2+ Mp Mp
k=1
b2k−
⎛
⎝
Mp
k=1
ak2
⎞
⎠
⎛
⎝
Mp
k=1
12
⎞
⎠ −
⎛
⎝
Mp
k=1
b2k
⎞
⎠
⎛
⎝
Mp
k=1
12
⎞
⎠
⎫⎬
⎭
= Mp
⎧⎨
⎩Mp
Mp
k=1
ak2+ Mp Mp
k=1
b2k− Mp Mp
k=1
ak2− Mp Mp
k=1
b2k
⎫⎬
⎭= 0,
where the equality holds in Cauchy–Schwarz inequality if and only if ak= C1and bk = C2
for all k= 1, 2, . . . , Mpand for some constants C1and C2, both independent of k. However, it it will never occur in our case. This leads to det(RR) > 0, and thus RR is invertible.
Note that we can efficiently compute the inverse of RR because it is only a 3× 3 matrix.
Furthermore, we can show the following result:
Theorem 3.1 The matrix A is symmetric and negative-definite.
Proof It is obvious that matrix A is symmetric. Using the facts that A is negative-definite and both B Band S Sare positive semi-definite, we only need to show that I− R(RR)−1R is a positive semi-definite matrix. A direct computation,
I− R(RR)−1R
I− R(RR)−1R
= I − R(RR)−1R, (41)
shows that I− R(RR)−1Ris a projection matrix, and then its eigenvalues should be 0 or 1. This leads to the assertion and thus completes the proof.
Thanks to this property of matrix A, we can solve (38) by using some efficient linear solvers, such as the preconditioned conjugate gradient and the algebraic multigrid methods.
Once we obtain the new velocity field un+1on the fluid grid, we can interpolate the new velocity to the marker points by (30) and move the Lagrangian markers to new positions by
Xn+1k = Xnk+ ΔtUn+1k . (42)
For the solid particle motion, the particle center velocity Vcand the angular velocity com- ponentω will be obtained by (36). Instead of first using the interpolation formula to find the new velocity and then updating the particle surface points, we simply adopt the idea of rigid body motion to determine the particle surface position (cf. [11]). Indeed, we compute the new particle center and the rotational angle by
Ync+1= Ync+ ΔtVnc+1, θn+1= θn+ Δtωn+1 (43) and use them to determine the new particle surface position as
Yn+1k = Ync+1+
cosθn+1− sin θn+1 sinθn+1 cosθn+1
Y0k. (44)
Finally, we remark that the proposed penalty IB approach can be applied to the transient Stokes flow with multiple compound vesicles as well. Similar to the above procedure, the resulting large linear system of equations is still tractable to block elimination, no matter how many compound vesicles are involved in the fluid, so that we only need to solve a reduced system with the velocity being the single unknown. Consequently, for multiple vesicle problems, the computational cost of the proposed penalty method is mainly devoted to solving a reduced linear system of the velocity approximation at each time step. This is one of the advantageous features of the proposed method.
4 Properties of the Penalty Immersed Boundary Method
In this section, we will discuss some properties of the penalty IB method.
4.1 Skew-adjointness and Adjointness in the Continuous Case
Let us first define the following spreading operators S1(σ ) and S2(F) and the inner products of functions onΩ, Γ and ∂ P:
S1(σ ) =
Γ
∂
∂s(σ τ)δ(x − X(s, t))ds, S2(F) =
∂ P Fδ(x − Y(α, t))dα, T1(u) = ∂U
∂s · τ, T2(u) =
Ωu(x)δ(x − Y(α, t))dx, u, v Ω =
Ωu· vdx, f, g Γ =
Γ f(s)g(s)ds, φ, ϕ ∂ P =
∂ Pφ(α)ϕ(α)dα.
Then we have
S1(σ ), u Ω =
Ω
Γ
∂
∂s(σ τ)δ(x − X(s, t))ds
· u(x)dx
=
Γ
∂
∂s(σ τ) ·
Ωu(x)δ(x − X(s, t))dx
ds
= −
Γσ
τ ·∂U
∂s
ds= σ, −T1(u) Γ, (45) which means that S1andT1are skew-adjoint. We also have
S2(F), u Ω =
Ω
∂ PFδ(x − Y(α, t))dα
· u(x)dx
=
∂ PF·
Ωu(x)δ(x − Y(α, t))dx
dα
= F, V ∂ P = F,T2(u) ∂ P, (46)
that is, S2andT2are adjoint. Moreover, since we assume the elastic tension and the particle surface force in the form
σ (s, t) = σ0
∂ X
∂s(s, t)
−
∂ X
∂s(s, 0)
, F = k0(Y − Y),
the associated potential energies are given by Eσ(t) =
Γ
σ0
2
∂ X
∂s (s, t)
−
∂ X
∂s(s, 0)
2
ds, (47)
EF(t) =
∂ P
k0
2|Y(α, t) − Y(α, t)|2dα. (48) Therefore, combining (47) with (45) and (18), and (46) with (20) and (21), we have
dEσ
dt =
Γ
1 2σ0
∂σ2
∂t ds=
Γσ1 σ0
∂σ
∂t
ds
=
Γσ τ ·∂U
∂s
ds= − S1(σ ), u Ω (49)
and
S2(F), u Ω = F, V ∂ P = F, Vc+ ωr ∂ P = Vc· F, 1 ∂ P+ ω F, r ∂ P = 0. (50) From (50) and (19), we obtain
!
F, V + 1 k0
∂ F
∂t
"
∂ P = 0, which leads to
dEF
dt =
! F, 1
k0
∂ F
∂t
"
∂ P = −# F, V$
∂ P. (51)
According to (49) and (51), we can conclude that the negative rate of change of potential energy of the elastic interface is equal to the work done by the interface on the fluid, while the work done by surface force F pulling the virtual solid particle, characterized by Y , back to the true particle position Y equals the negative rate of change of potential energy of the particle.
4.2 Skew-adjointness and Adjointness in the Discrete Case
We first define the following discrete operators S1h, S2h,T1h,T2h, and discrete inner products:
S1h(σn+1) =
k
Ds(σn+1τn)kδh(x − Xnk)Δs, S2h(Fn+1) =
k
Fn+1k δh(x − Ynk)Δα,
T1h(un+1k−1/2) = Un+1k − Un+1k−1
Δs · τnk−1/2, T2h(un+1k ) =
x
un+1(x)δh(x − Ynk)h2, w, v Ωh =
x
w(x) · v(x)h2, f, g Γh =
k
fk−1/2gk−1/2Δs, φ, ϕ ∂ Ph =
k
φk−1/2ϕk−1/2Δα.
Then we can verify that S1handT1hare skew-adjoint, while S2h andT2h are adjoint as follows:
S1h(σn+1), un+1 Ωh =
x
k
Ds(σn+1τn)kδh(x − Xnk)Δs
· un+1(x)h2
=
k
Ds(σn+1τn)k·
x
un+1(x)δh(x − Xnk)h2 Δs
=
k
Ds(σn+1τn)k· Un+1k Δs
=
k
σkn+1+1/2τnk+1/2− σkn+1−1/2τnk−1/2
Δs · Un+1k Δs
= −
k
σk−1/2n+1 Un+1k − Un+1k−1
Δs · τnk−1/2Δs
= σn+1, −T1h(un+1) Γh (52)
and
S2h(Fn+1), un+1 Ωh =
x
k
Fnk+1δh(x − Ynk)Δα
· un+1(x)h2
=
k
Fn+1k ·
x
un+1(x)δh(x − Ynk)h2 Δα
= Fn+1,T2h(un+1) ∂ Ph. (53) Indeed, based on these two properties (52) and (53), we are able to rescale the discretiza- tions (24)–(31) to be a symmetric system of linear equations.
4.3 The Difference of Local Stretching Factors of Two Successive Time Steps is of O(Δt) From the moving equation of the elastic interface (42), we have
Xnk+1− Xnk−1+1
Δs = Xnk− Xnk−1
Δs + ΔtUnk+1− Unk−1+1
Δs (54)
and then
DsXn+1k−1/2= DsXnk−1/2+ Δt DsUn+1k−1/2. (55) By the triangle inequality, we obtain
|DsXn+1k−1/2| − |DsXnk−1/2| ≤ |DsXn+1k−1/2− DsXnk−1/2| ≤ Δt|DsUn+1k−1/2|, (56) which shows the difference of local stretching factors of two successive time steps is of O(Δt) in the discrete case. Now letting Δs → 0, we obtain from (55) that
∂
∂sXn+1k−1/2= ∂
∂sXnk−1/2+ Δt ∂
∂sUn+1k−1/2, (57)
and this leads to
∂
∂sXn+1k−1/2
2=
∂
∂sXnk−1/2
2+ 2Δt∂
∂sXnk−1/2· ∂
∂sUn+1k−1/2+ (Δt)2
∂
∂sUn+1k−1/2
2. (58) On the other hand, by virtue of (26), we have
DsUnk−1/2+1 · DsXnk−1/2= 1 σ0Δt
σk−1/2n+1 − σkn−1/2
|DsXnk−1/2|
= 1
Δt
∂
∂sXn+1k−1/2
−
∂
∂sXnk−1/2
|DsXnk−1/2|, (59) which implies
Δt ∂
∂sXnk−1/2· ∂
∂sUnk−1/2+1 =
∂
∂sXnk−1/2+1
−
∂
∂sXnk−1/2
∂
∂sXnk−1/2
. (60) Finally, combining (58) with (60), we have
∂
∂sXn+1k−1/2
−
∂
∂sXnk−1/2
2
= (Δt)2
∂
∂sUn+1k−1/2
2. (61)
That is, for the spatial continuous case, the difference of local stretching factors is of O(Δt), too.
