Tank-treading motion under shear flow

Unlike the previous subsection that we focus on the numerical convergence test for our present scheme, in this subsection, we consider the physical transient deformation of an inextensible interface subject to a simple shear flow. As mentioned before, the motivation of this test is to mimic the simulation of the vesicle dynamics which has a lot of applications in bio-fluid problems.

As in previous test, we put an inextensible interface Γ with initial configuration (X(s), Y (s)) = (0.18 cos(s), 0.5 sin(s)) under a shear flow (u, v) = (γ y, 0) in a fluid domain Ω = [−1, 1] × [−1, 1], see Figure 4.1 in detail. The shear rate γ is chosen to be γ = 1 and the fluid viscosity is µ = 1. The mesh used is 128 × 128 and the residual tolerance for is 10⁻⁴. It is worth mentioning that the elastic tension σ computed in previous numerical experiments [71, 28, 35] tends to oscillate along the interface which makes the conjugate gradient method for solving in Eq. (4.25) difficult to converge if the residual tolerance is too large. It will be more appealing if we can find a good preconditioner for the matrix in Eq. (4.25) so that PCG method converges faster.

This issue should be investigated further in later work. Figure 4.2 shows the interface configuration at three different times t = 0.0625, 1.25, 3.125, while Figure 4.3 shows the corresponding elastic tension σ plotted counterclockwise along the interface. (The starting point of zero length in Figure 4.3 is marked by x in Figure 4.2.) One can indeed see some slight oscillations of the tension along the interface as seen in previous literature [71, 28, 35]. Furthermore, one can see that the tension has smallest values at the interface positions where the curvatures are largest (both tips) which also in agreement with those in previous literature.

To be more physically realistic, we now run the problem to a longer time. It is well-known that the equilibrium dynamics of inextensible interface or vesicle under a simple shear flow undergoes a tanking-treading motion if the viscosity contrast under a certain threshold [24]. Here, by tank-treading motion we mean that the

−1 −0.5 0 0.5 1

−1

−0.8

−0.6

−0.4

−0.2 0 0.2 0.4 0.6 0.8 1

x−axis

y−axis

t = 0.0625 t = 1.25 t = 3.125

Figure 4.2: The motion of an inextensible interface in a shear flow with initial con-figuration (X(s), Y (s)) = (0.18 cos(s), 0.5 sin(s)).

0 0.5 1 1.5 2

−0.5 0 0.5 1 1.5 2

vesicle length

tension σ

t = 0.0625 t = 1.25 t = 3.125

Figure 4.3: The plots of tension σ along the interface.

configuration of the interface or vesicle remains stationary while there is a tangential motion along the interface. Figure 4.4 shows the evolutional motion of the interface (bi-concave shape initially) at different times. We can see that after some time, the interface shape does not seem to change at all; however, the Lagrangian point along the interface marked by ”∗” moves along with its tangential velocity. The streamlines at three different chosen times are shown in Figure 4.5 in which we can see no normal motion in equilibrium. This tank-treading motion is good agreement in previous studies [71, 66, 28, 35].

T = 0 T = 0.015625 T = 0.51563 T = 1.0156

T = 7.5781 T = 10.0781 T = 12.5781 T = 15.0781

T = 17.5781 T = 20.0781 T = 22.5781 T = 25.0781

Figure 4.4: The tank-treading motion of an inextensible interface under a shear flow.

Figure 4.5: The streamlines of the flow along with an interface.

As discussed in previous literature [71, 26, 25, 66, 28], the motion of a steady or an equilibrium inextensible interface (or vesicle) can be characterized by both the inclination angle θ between the long axis of interface and the flow direction, and the tank-treading frequency f = 2π± R

Γ dl

uτ of the revolution, where u_τ is the tangential velocity component. The inclination angle has been founded to be strongly dependent on the reduced area ν. However, the angle is independent of the shear rate γ. This behavior has been verified in the left panel of Figure 5.2 and is in good agreement with previous studies [71, 26, 25, 66, 28] which shows the steady inclination angle (θ/π) versus the reduced area (ν) with different shear rates γ = 1, 5, 10. We have observed that the inclination angle increases with the reduced area but is nearly independent of the shear rate. The right panel of Figure 5.2 shows that the tank-treading frequency f versus the reduced area ν for different shear rates. One can see that as the shear

0.5 0.6 0.7 0.8 0.9 0.1

0.12 0.14 0.16 0.18 0.2

reduced area ν

θ/π

γ = 1 γ = 5 γ = 10

0.5 0.6 0.7 0.8 0.9

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

reduced area ν

f

γ = 1 γ = 5 γ = 10

Figure 4.6: The inclination angles θ/π (left) and the tank-treading frequency f = 2π±R

Γ dl u · τ

(right) versus reduced areas ν with different shear rates for the tank-treading motions of an inextensible interface in a shear flow.

rate increases, the tangential motion becomes stronger; thus, the frequency becomes larger. Moreover, by fixing the shear rate, if the interface has larger reduced area then it has larger frequency as well. Again, our numerical results are in a good agreement with previous studies in literature.

Chapter 5

Unconditionally energy stable immersed boundary method

In this chapter, we propose an unconditionally stable IB method to simulate vesicle dynamics. The vesicle is assumed to be nearly inextensible, that is, the magnitude of time changing of stretching factor is not exactly zero but quite small. Under this penalty method, the tension of vesicle has an explicit form so we do not need to solve fluid variables and surface tension simultaneously (differ form previous chapter).

However, since the cost of penalty method contributes to interfacial stiffness, there is a harsh restriction on choosing time step ∆t to maintain numerical stability.

As known in literature [62, 38, 58, 50], the IB method suffers a time step restriction to maintain numerical stability. This restriction becomes more stringent when the elastic force is stiff and the force spreading occurs at the beginning of each time step (an explicit scheme). Note that, such time step restriction cannot be elevated even the fluid solver is discretized in a semi-implicit manner (the explicit differencing of the advection term and implicit differencing of the diffusion term). Rather than performing the force spreading at the beginning of the time step, one might consider to perform the procedure at the intermediate (semi-implicit scheme) or even at the end of time step (implicit scheme). There have been many attempts to reduce the stiffness or to overcome this difficulty of time step restriction in the past decade [44, 46, 20, 21, 6, 7, 16]. However, there is always a trade-off between the stability and efficiency of those algorithms. In this chapter, we shall propose a new semi-implicit scheme that can be solved quite efficiently (the resultant linear system is symmetric positive definite) and the time step size can be significantly increased.

This chapter is organized as follows. In Section 5.1, we introduce the formulation of incompressible Navier-Stokes equations with an ordinary immersed elastic interface (tension has explicit formulation). We then develop semi-implicit IB schemes based on the projection method for the fluid solver in Section 5.2. We also show those developed schemes to be unconditionally energy stable. Then we modify those semi-implicit schemes to be more efficient so that the resultant linear system is symmetric

positive definite. Finally we apply present method to simulation of vesicle dynamics, the numerical results are given in Section 5.6.

5.1 Governing equations

We begin by stating the mathematical formulation of the Navier-Stokes flow with an immersed boundary (or an interface). We consider a moving, immersed, elastic boundary Γ(t) which exerts forces into an incompressible fluid in a fixed fluid domain Ω. We assume that the fluids inside and outside of the boundary are the same so the governing equations in immersed boundary formulation can be written as follows.

ρ µ∂u

∂t + (u · ∇)u

¶

+ ∇p = µ∆u + Z

Γ

F(s, t)δ(x − X(s, t)) ds in Ω, (5.1)

∇ · u = 0 in Ω, (5.2)

∂X

∂t (s, t) = U(s, t) = Z

Ω

u(x, t)δ(x − X(s, t)) dx on Γ. (5.3) Eqs. (5.1) and (5.2) are the familiar incompressible Navier-Stokes equations with a singular force arising from the immersed boundary. Eq. (5.3) simply represents that the immersed boundary moves along with the local fluid velocity (the interfacial veloc-ity). Here, the interfacial velocity U is simply an interpolation of the fluid velocity at the immersed boundary. The interaction between the fluid and the immersed bound-ary is linked by the two-dimensional Dirac delta function δ(x) = δ(x)δ(y). In this chapter, the elastic force mainly comes from the tension which satisfies the Hooke’s law as

F(s, t) = ∂

∂s(στ ), σ = σ₀(|X_s| − r₀), (5.4)

where σ₀ is the elastic coefficient, and r₀ is the rest length. Note that, the rate of local stretching factor can be derived to satisfy the following equation [33]

∂

∂t|X_s| = (∇_s· U) |X_s| = ∂U

∂s · τ . (5.5)

5.2 Numerical algorithm

Let ∆t be the time step size and the superscript index be the time step level. At the beginning of each time level n, the boundary configuration Xⁿ_k, and the unit tangent τⁿ_k−1/2 are all given. To simplify our notations, we define the discrete spreading

operator acting on the tension and discrete surface divergence operator acting on the velocity as follows.

S_h^m[σ](x) = S_h^m[σ₀(|X_s|−r₀)](x) = XM k=1

D_s(σ₀(|X_s| − r₀)τ^m)_kδ_h(x−X^m_k)∆s, (5.6)

T_h^m[u](s_k−1/2) = U_k− U_k−1

∆s · τ^m_k−1/2 U_k=X

x

u(x)δ_h(x − X^m_k)h². (5.7) The superscript m stands for the information of immersed boundary position adopted at time level m∆t. Note that, the above discrete operators are skew-adjoint as hS_h^m[σ], uiΩh = hσ, −T_h^m[u]iΓh. It is also worthy mentioning that the skew-adjointness property plays an important role in our energy stability analysis later.

In present chapter, our numerical discretization and energy stability analysis are based on unsteady Stokes equations instead of Navier-Stokes. For the latter, the nonlinear advection term can be treated explicitly during the time evolution with moderate CFL condition. Alternatively, one can split the Navier-Stoke equations into an advection part and unsteady Stokes part in which the advection equation is solved by an alternating direction implicit (ADI) method to maintain the unconditionally numerical stability [21].

Here we introduce two time-integration schemes; namely, the backward Euler (BE) and Crank-Nicholson (CN) scheme in conjunction with the projection method [31] as follows.

• Backward Euler (BE) scheme ρu^∗− uⁿ

∆t = µ∆hu^∗ + S_hⁿ[σ0(|Xs|ⁿ⁺¹− r0)] in Ω_h (5.8)

∆_hφ = 1

∆t∇_h· u^∗ in Ω_h; ∂φ

∂n

¯¯

¯¯

∂Ωh

= 0, (5.9)

uⁿ⁺¹ = u^∗− ∆t∇hφ in Ωh (5.10)

∇_hpⁿ⁺¹= ρ∇_hφ − µ∇_h(∇_h· u^∗) in Ω_h (5.11)

|X_s|ⁿ⁺¹− |X_s|ⁿ

∆t = T_hⁿ[uⁿ⁺¹] on Γ_h (5.12)

Uⁿ⁺¹_k =X

x

u(x)ⁿ⁺¹δh(x − Xⁿ_k) h² on Γh, (5.13)

Xⁿ⁺¹_k = Xⁿ_k+ ∆t Uⁿ⁺¹_k on Γh. (5.14)

• Crank-Nicholson (CN) scheme The spatial operators ∇_h and ∆_h are the standard second-order central difference approximations to the gradient and Laplacian. Here, δh is a smoother version of discrete delta function developed in [70]. Notice that, unlike the explicit scheme, both backward Euler and Crank-Nicholson schemes treat the stretching factor |X_s|ⁿ⁺¹ as an unknown solution which is closed by the discretization of Eq. (6.14). Thus, in terms of immersed boundary force computation, the above two schemes are semi-implicit.