Maximal time step comparison

Now, let us investigate the numerical stability of the present schemes by testing the maximal time steps for the explicit scheme (EP), the modified BE, and CN scheme. Here, we use three different elastic coefficients σ₀ = 10⁷, 10⁸, 10⁹ to study the numerical stability by comparing the maximal time step that can be used in each scheme. To determine the maximal time step, we make sure that the vesicle boundary behaviors are reasonable in which the relative errors of area and perimeter are both within 1% and no numerical instability occurs. The numerical parameters are chosen as same as in previous convergence study except the initial vesicle configuration is chosen as X(s) = (0.1 sin(s), 0.5 cos(s)). Table 5.2 shows the maximal time steps for three different schemes. One can see that the time step for both modified BE and CN schemes can be chosen as 3-4 order larger than the time step used in the explicit scheme. For explicit scheme, as the grid becomes finer or the elastic coefficient σ₀ becomes larger, the time step should become smaller accordingly in order to maintain numerical stability. However, for both BE or CN schemes, we can always set ∆t = O(h) to maintain the desired numerical stability.

Table 5.3 shows the average CPU time (in seconds) for each time step and total CPU time for the computation up to T = 1. One can see that, the present modified BE and CN schemes outperform the explicit scheme in terms of total CPU time.

m, n σ0= 10⁷ σ0= 10⁸ σ0= 10⁹

EP BE CN EP BE CN EP BE CN

128 2.22e-5 1.56e-2 1.56e-2 4.36e-6 1.56e-2 1.56e-2 1.02e-6 1.56e-2 1.56e-2 256 1.03e-5 7.81e-3 7.81e-3 1.74e-6 7.81e-3 7.81e-3 3.81e-7 7.81e-3 7.81e-3 512 4.74e-6 3.91e-3 3.91e-3 7.44e-7 3.91e-3 3.91e-3 1.53e-7 3.91e-3 3.91e-3

Table 5.2: Maximum time steps for the explicit(EP), modified BE and CN schemes with different elastic coefficients σ₀.

m, n EP BE CN

∆t Average Total ∆t Average Total ∆t Average Total

128 1.02e-6 0.07 68627 1.56e-2 0.13 8 1.56e-2 0.13 8

256 3.81e-7 0.16 419947^∗ 7.81e-3 0.39 50 7.81e-3 0.43 55

512 1.53e-7 0.45 2941176^∗ 3.91e-3 1.25 320 3.91e-3 1.45 371

Table 5.3: The average CPU time (in seconds) of each time step for different schemes.

The total time with ”*” means the estimated value.

5.6.3 Tank-treading motion under shear flow

To be more physically realistic, we now consider the transient motion of a vesicle with the initial configuration X(s) = (0.1 sin(s), 0.5 cos(s)) under shear flow u = (γy, 0). The elastic coefficient is chosen as σ₀ = 10⁵ and the bending coefficient c_b = 0.01. It is well-known that the equilibrium dynamics of vesicle under a simple shear flow undergoes tanking-treading motion if the viscosity contrast is under a certain threshold [24]. Here, by tank-treading motion we mean that the configuration of the vesicle remains stationary while there is a tangential motion along the vesicle boundary. Figure 5.1 shows the evolutional motion of the vesicle at different times.

One can see that after some time, the vesicle shape does not change; however, the Lagrangian point (marked by ∗) along the interface moves along with its tangential velocity. This tank-treading motion is in good agreement with previous studies [71, 66, 28, 35].

As discussed in literature [71, 26, 25, 66, 28], the motion of a steady vesicle under shear flow can be characterized by the inclination angle θ between the long axis of vesicle and the flow direction, and the tank-treading frequency f = 2π± R

Γ dl

uτ, where u_τ is the tangential velocity component. The inclination angle has been found to be strongly dependent on the reduced area ν. As the reduced area increases, the inclination angle will increase too. However, the angle is independent of the shear rate γ. This behavior has been verified in the left panel of Figure 5.2 which shows the steady inclination angle (θ/π) versus the reduced area (ν) with different shear rates γ = 1, 5, 10. Our numerical results are again in good agreement with those

T = 0.125 T = 0.375 T = 0.625 T = 0.875

T = 13.5 T = 16.5 T = 19.5 T = 22.5

T = 25.5 T = 28.5 T = 31.5 T = 34.5

Figure 5.1: The tank-treading motion of a vesicle under shear flow.

0.5 0.6 0.7 0.8 0.9

0.1 0.12 0.14 0.16 0.18

reduced area ν

θ/π

γ = 1 γ = 5 γ = 10

0.5 0.6 0.7 0.8 0.9

0 1 2 3 4 5

reduced area ν

f

γ = 1 γ = 5 γ = 10

Figure 5.2: 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 motion of an inextensible vesicle in a shear flow.

obtained in previous studies [66, 28]. 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 rate increases, the tangential motion becomes stronger; thus, the frequency becomes larger. Moreover, by fixing the shear rate, if the vesicle has larger reduced area then it has larger frequency as well. Again, our numerical results are in good agreement with those obtained in previous studies [66, 28].

Chapter 6

Simulating three-dimensional axisymmetric vesicles

There are two different approaches to enforce the local inextensibility constraint in literature. The first one needs to discretize the whole equations first (regardless of using boundary integral, finite element or finite difference method) and then to solve the tension and fluid variables simultaneously. Based on different time stepping schemes, this kind of approach can be explicit or semi-implicit depending on how we treat the force computations. There usually exists a trade-off between the time-step stability and efficiency in those algorithms. Some recent works using boundary integral method [66, 67, 68, 4, 72], level set method [60, 36] et. al. fall into this category, just to name a few. Recently, we have also introduced a linearly semi-implicit scheme to solve a 2D Stokes flow with an inextensible interface enclosing a solid particle [37]. The scheme solves the fluid variables, the tension and the particle surface force altogether. There are two advantages of the method as follows: (1) the scheme preserves some operators symmetrically so the resultant linear system is symmetric; (2) the local inextensibility error can be estimated analytically. Although one can increase the time step significantly, more improvements have to be done on the iterative methods for solving such resultant linear system, especially finding good pre-conditioners.

Another approach is called penalty method. In [28], the authors introduced a dual Lagrangian immersed boundaries to represent the vesicle boundary. One vesicle boundary interacts with the fluid dynamics directly and the other vesicle boundary follows the equations of the vesicle dynamics, including the inextensibility constraint, without a direct interaction with the fluid dynamics. The two boundaries are con-nected by the penalty forces which act on both boundaries by Newton’s third law.

The tension must still be solved in above approach.

However, one can even simplify the above penalty approach by introducing a spring-like force to replace the tension and to keep the surface area from dilating.

Thus, the tension or Lagrange’s multiplier for inextensibility is no need to be solved

as part of the solution. Although the inextensibility constraint is not exactly satisfied, we can get physically reasonable results as long as the penalty number (stiffness number) is chosen large enough. This approach has been adopted successfully by incorporating the level set method to study 2D or 3D [39, 41, 11] vesicle dynamics in finite element framework. A partially implicit sharp interface method in [5] also used the above approach to test a 2D case with an inextensible interface. In this chapter, we use the above similar idea to study 3D axisymmetric vesicle dynamics in immersed boundary framework. One should notice that all the calculations about the spatial derivatives, geometrical quantities (mean and Gaussian curvatures), the bending and elastic forces, and the Navier-Stokes solver have non-trivial differences comparing with our previous 2D works.

The rest of the chapter are organized as follows. In Section 2, we will describe the fluid-vesicle governing equations in axisymmetric form based on immersed boundary formulation. Instead of introducing a Lagrange’s multiplier to enforce the vesicle inextensibility constraint, we modify the model by adopting a spring-like tension to make the vesicle boundary nearly inextensible so that solving the unknown tension can be avoided. Since higher-order spatial derivatives are needed in computations of interfacial geometrical quantities, we use Fourier spectral approximation to represent the interface. A detailed numerical algorithm is described in Section 6.2. A series of numerical tests to demonstrate the accuracy and reliability of the present scheme are presented in Section 6.3.

6.1 Governing equations

We consider a single 3D axisymmetric inextensible vesicle Γ(t) suspended in a viscous incompressible Navier-Stokes fluid domain Ω. Here, the vesicle is time-dependent and is symmetry in the θ direction so that the vesicle surface has the parametric form as X(s, θ, t) = (R(s, t) cos θ, R(s, t) sin θ, Z(s, t)), (6.1) where the parameters (s, θ) are in [0, π] × [0, 2π]. One should also note that the vesicle membrane force which is defined later has the same symmetric form as in Eq. (6.1). Under such symmetry, one can regard the vesicle boundary as the surface of revolution around the z-axis with the radius R(s, t) so that it can be simplified into a two-dimensional immersed boundary X(s, t) = (R(s, t), Z(s, t)), s ∈ [0, π] in the 2D meridian r − z domain.

Under the same assumption of axisymmetry, the 3D Navier-Stokes equations can be simply written in a 2D manner using the axisymmetric cylindrical coordinates x = (r, z). In the following, we denote u = (u, w) as the velocity field defined on a 2D meridian domain Ω = {(r, z)|0 < r ≤ a, c ≤ z ≤ d}, where u and w are the radial (r-coordinate) and axial (z-coordinate) velocity components, respectively.

We also denote U = (U, W ) as the corresponding velocity components on the vesicle

boundary. The mathematical equations of motion consist of a viscous incompressible fluid in a domain Ω containing an immersed, inextensible, massless vesicle boundary (or interface) Γ which can be written in an immersed boundary formulation as

ρ Here, we assume the fluid has the same viscosity and density inside and outside of the vesicle boundary. For axisymmetry, the gradient and divergence operators are defined as

thus, the Laplace operator is

∆ = ∇ · ∇ = 1

Eqs. (6.2)-(6.3) are the familiar incompressible Navier-Stokes equations with a singular force term F arising from the vesicle membrane force. Eq. (6.4) represents the inextensibility constraint of the vesicle surface which is equivalent to the zero sur-face divergence of the velocity along the intersur-face, see the explanation later. Eq. (6.5) simply explains that the interface moves along with the local fluid velocity (the in-terfacial velocity). Here, the inin-terfacial velocity U is the interpolation of the fluid velocity at the interface defined as in traditional IB formulation. The interaction be-tween the fluid and the interface is linked by the two-dimensional Dirac delta function δ(x) = δ(r)δ(z).