An immersed boundary method for simulating the dynamics of three-dimensional axisymmetric vesicles in Navier-Stokes flows

Contents lists available atScienceDirect

Journal of Computational Physics

www.elsevier.com/locate/jcp

An immersed boundary method for simulating the dynamics

of three-dimensional axisymmetric vesicles in Navier–Stokes

flows

Wei-Fan Hu

a

, Yongsam Kim

b

, Ming-Chih Lai

c

,

∗

a_{Department of Applied Mathematics, National Chiao Tung University, 1001, Ta Hsueh Road, Hsinchu 300, Taiwan} b_{Department of Mathematics, Chung-Ang University, Dongjak-Gu, Heukseok-Dong 221, Seoul 156-175, South Korea}

c_{Center of Mathematical Modeling and Scientific Computing & Department of Applied Mathematics, National Chiao Tung University, 1001,}

Ta Hsueh Road, Hsinchu 300, Taiwan

a r t i c l e

i n f o

a b s t r a c t

Article history:

Received 30 October 2012

Received in revised form 14 June 2013 Accepted 13 October 2013

Available online 18 October 2013

Keywords:

Immersed boundary method Inextensible interface Axisymmetric vesicle Navier–Stokes equations

In this paper, we develop a simple immersed boundary method to simulate the dynamics of three-dimensional axisymmetric inextensible vesicles in Navier–Stokes flows. 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 for the unknown tension can be avoided. We also derive a new elastic force from the modified vesicle energy and obtain exactly the same form as the originally unmodified one. In order to represent the vesicle boundary, we use Fourier spectral approximation so we can compute the geometrical quantities on the interface more accurately. A series of numerical tests on the present scheme have been conducted to illustrate the applicability and reliability of the method. We first perform the accuracy check of the geometrical quantities of the interface, and the convergence check for different stiffness numbers as well as fluid variables. Then we study the vesicle dynamics in quiescent flow and in gravity. Finally, the shapes of vesicles in Poiseuille flow are investigated in detail to study the effects of the reduced volume, the confinement, and the mean flow velocity. The numerical results are shown to be in good agreement with those obtained in literature.

©2013 Elsevier Inc. All rights reserved.

1. Introduction

Vesicles dynamics in fluid flow has become a quite active research area recently in the communities of soft matter physics and computational fluid mechanics. Indeed, the understanding of vesicle behaviors in fluids might lead to a better knowledge of red blood cells (RBCs) in bloods since they both share some similar mechanical behaviors. A vesicle can be visualized as a liquid droplet within another liquid enclosed by a lipid membrane with the size about 10 μm. Such lipid membrane consists of tightly packed lipid molecules with hydrophilic heads facing the exterior and interior fluids and with hydrophobic tails hiding in the middle and thus forms a bilayer phospholipid with thickness about 6 nm. It is well-known that the phospholipid membrane exhibits a resistance against area dilation and bending; therefore, it is quite natural to regard this membrane as an inextensible surface with mechanical properties defined by some energy functional, see the next section. Thus, the dynamics of moving vesicles can be determined by their membrane elasticity, inextensibility, and hydrodynamical force. For the past years, the vesicle problems have been extensively explored by the classic [14]and

*

Corresponding author.

E-mail addresses:weifanhu.am95g@g2.nctu.edu.tw(W.-F. Hu),kimy@cau.ac.kr(Y. Kim),mclai@math.nctu.edu.tw(M.-C. Lai). 0021-9991/$ – see front matter ©2013 Elsevier Inc. All rights reserved.

small deformation theories[25,21,9], flow experiments[12,7], and computer simulations[13,3,26,34,35,15,4,30,36,39]; just to name a few recent ones. Readers who are interested in more details on those relevant works can refer to the recent article [39]. In this paper, we will focus only on numerical computation of the vesicle problem; in particular, on the axisymmetric vesicles in Navier–Stokes flow.

The study of vesicle dynamics is a fluid-structure interaction problem. The numerical simulation of the problem consists of not only solving a two-phase flow but also requiring to enforce an inextensibility constraint along the surface, which makes the problem challenging to attack. Among the numerical methods for simulating vesicle problems in literature, one can characterize those methods by how the membrane surface (or interface) is represented and what the fluid solver is used. Based on this characterization, several methods have been developed such as boundary integral method[13,34,35,4,36,39], level set method [23,30,20,24,8], phase field method [6,3,24], particle collision method [26], immersed interface method [18,33], and immersed boundary method[15,16]. The boundary integral method[34–36,39] using spectral representation for the vesicle front tracking has been quite successfully applied to Stokes flow in recent years. Although most of physical applications of vesicles are in Stokes regime, it is interesting to observe that the inertial effect inhibits vesicles from tumbling [19]. Furthermore, as mentioned in[19], human arterioles with shear rate about 8000 s−1 results in the Reynolds number

O

(1)

which cannot be neglected in simulations. Therefore, it is quite natural to take the inertial effect into account, in which case Navier–Stokes equations should be solved.

Another numerical issue of simulating the vesicles in flows comes from the surface inextensibility. Unlike merely two-phase flow simulations, here, the surface tension has different physical meaning which is unknown and coupled with fluid dynamics via the membrane force. In fact, the tension acts like a Lagrange’s multiplier function to enforce the lo-cal inextensibility constraint along the surface which is exactly the same role played by the pressure to enforce the fluid incompressibility in Navier–Stokes equations. Therefore, the development of an efficient and accurate algorithm for vesicle dynamics remains quite challenging; not to mention the vesicle surface moves along with the surrounding fluid as well.

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. This 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 [34–36,4,39] or level set method [30,20] 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[16]. 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[15], the authors introduced a dual representation of Lagrangian immersed boundaries for 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 connected by the penalty forces which act on both boundaries by Newton’s third law. The tension must still be solved in this approach.

However, one can simplify the penalty approach even more 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 does not have to be solved for a part of the solution. Although the inextensibility constraint is not exactly satisfied, we can get physically reason-able results as long as the penalty number (stiffness number in this paper) is chosen large enough. This approach has been adopted successfully by incorporating the level set method to study 2D or 3D[22,23,8]vesicle dynamics in finite element framework. A partially implicit sharp interface method in[2]also used this approach to test a 2D inextensible interface. In this paper, we adopt a similar idea to study 3D axisymmetric vesicle dynamics in immersed boundary framework.

We summarize the contributions of the present work as follows.

•

We simplify our previous 2D penalty immersed boundary method for simulating vesicle dynamics developed in [15] in a way that no dual representation of immersed boundary points are needed. In addition, the tension or Lagrange’s multiplier for inextensibility is no need to be solved as part of the solution since it can be approximated by a spring-like elastic tension to keep the local surface dilating factor close to that of the original configurations (nearly inextensible).

•

We formulate the vesicle bending and elastic forces explicitly in axisymmetric setting. The new elastic force, which is derived by taking the variational derivative of the present modified elastic energy, has the exactly same form as the original tension force despite the fact that the definitions of tension

σ

are different.

•

We adopt the spectral Fourier discretization to represent the vesicle boundary and to compute the geometrical quanti-ties used in the computations of bending and elastic forces. Although the spectral boundary representation is previously used in other approaches such as the boundary integral method, it is newly applied to the immersed boundary method. Moreover, we have made quantitative comparisons with second-order finite difference schemes for the approximations of the mean and Gaussian curvatures, and the surface Laplacian of mean curvature.

•

We present several numerical results in Section 4to illustrate the applicability and reliability of the present scheme. In particular, we simulate the problem of vesicles in Poiseuille flows and study the effects of the reduced volume, the confinement, and the mean velocity in detail. Our numerical results are qualitatively consistent with experimental data. The rest of the paper are organized as follows. In Section 2, we describe the fluid-vesicle governing equations in ax-isymmetric 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 de-tailed numerical algorithm is described in Section3. A series of numerical tests to demonstrate the accuracy and reliability of the present scheme are presented in Section4which is followed by conclusions and future work in Section5.

2. Governing equations in axisymmetric formulation

In this paper, 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

)



,

(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.(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 de-fined 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 of 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

ρ



∂

u

∂

t

+ (

u

· ∇)

u



= −∇

p

+

μ



u

+

π



0 F

(

s

,

t

)δ

2



x

−

X

(

s

,

t

)



ds in

Ω,

(2)

∇ ·

u

=

0 in

Ω,

(3)

∇

s

·

U

=

0 on

Γ,

(4)

∂

X

∂

t

(

s

,

t

)

=

U

(

s

,

t

)

=



Ω u

(

x

,

t

)δ

2



x

−

X

(

s

,

t

)



dx on

Γ.

(5)

Here, we assume that the fluid has the same viscosity and density inside and outside of the vesicle boundary. For axisym-metry in the cylindrical coordinate system, the gradient and divergence operators are defined as

∇ =



∂

∂

r

,

∂

∂

z



,

∇· =



1 r

∂

∂

rr

,

∂

∂

z



·,

(6)

and thus the Laplace operator is

= ∇ · ∇ =

1 r

∂

∂

r



r

∂

∂

r



+

∂

2

∂

z2

,

(7) and



u

= (

u

−

_ru2

,

w

).

Eqs.(2)–(3)are the familiar incompressible Navier–Stokes equations with a singular force term F arising from the vesicle membrane force. Eq. (4) represents the inextensibility constraint of the vesicle surface which is equivalent to the zero surface divergence of the velocity along the interface, see the explanation later. Eq.(5)simply says that the interface moves along with the local fluid velocity (the interfacial velocity). Here, the interfacial velocity U is the interpolation of the fluid velocity at the interface defined as in traditional IB formulation. The interaction between the fluid and the interface is linked by the two-dimensional Dirac delta function

δ

2

₍

_x

₎

_{= δ(}

_r

_)δ(

_z

_).

2.1. Vesicle boundary force

A vesicle membrane is known to be inextensible and exhibits a resistance against bending. Thus, the membrane energy can be modeled by two parts; namely, the Helfrich type energy Eb [10] to resist the bending of the membrane, and the

tension energy Eσ to enforce the inextensibility constraint[40]. So the total energy is

E

(

t

)

=

Eb

(

t

)

+

Eσ

(

t

)

=



Γ



cb 2H 2

₊

_σ



dS

,

(8)

where cb is the bending rigidity, H is the surface mean curvature, and

σ

is the tension which acts as a Lagrange multiplier

to enforce the inextensibility constraint Eq.(4). By taking the variational derivative to the surface energy(8), one can derive the vesicle boundary force FΓ

=

Fb

+

Fσ consisting of the bending force Fb and the elastic force Fσ as

Fb

=

cb



sH

+

2H



H2

−

K



n

,

Fσ

= ∇

s

σ

−

2H

σ

n

,

(9)

where K is the membrane Gaussian curvature, n is the unit outward normal vector,

∇

sis the surface gradient operator, and

sis the surface Laplacian operator. The detailed derivation of Eq.(9)can be found in[37].

Under the axisymmetric assumption, one can express H , K , and

sH in terms of R and Z explicitly (seeAppendix Ain

detail) as H

=

1 2



Zs R

|

Xs

|

+

RsZss

−

RssZs

|

Xs

|

3



,

K

=

Zs

(

RsZss

−

RssZs

)

R

|

Xs

|

4

,

(10)

sH

=

1 R

|

Xs

|

∂

∂

s



R

|

Xs

|

∂

H

∂

s



.

(11)

The unit tangent

τ

and the unit outward normal vector n can be hereby defined as

τ

(

s

,

t

)

=

Xs

|

Xs

|

=

(

Rs

,

Zs

)



R2 s

+

Z2s

,

n

(

s

,

t

)

=



(

Zs

,

−

Rs

)

R2 s

+

Zs2

.

(12)

Thus, by substituting the geometrical quantities on a surface into Eq. (9), we can derive easily an explicit form for the bending force Fb. Similarly, the tension force Fσ can also be explicitly written as

Fσ

= ∇

s

σ

−

2Hσn

=



1

|

Xs

|

∂

σ

∂

s



τ

−



Zs R

|

Xs

|

+

RsZss

−

RssZs

|

Xs

|

3



σ

n

=



RssZs

−

RsZss

|

Xs

|

3

σ

n

+

1

|

Xs

|

∂

σ

∂

s

τ



−

Zs R

|

Xs

|

σ

n

=

1

|

Xs

|



σ

∂

τ

∂

s

+

∂

σ

∂

s

τ



−

Zs R

|

Xs

|

σ

n

=

1

|

Xs

|

∂

∂

s

(

σ τ

)

−

Zs R

|

Xs

|

σ

n

.

(13)

The delta function in 3D Cartesian coordinates and axisymmetric cylindrical coordinates can be converted using the relation

δ

3

= δ

2

/

r where

δ

2 is the two-dimensional delta function in r–z plane. Therefore, the singular force term can be written as



Γ FΓ

δ

3

(

x

−

X

)

dS

=

π



0 FΓ

δ

2

(

x

−

X

)

R R

|

Xs

|

ds

=

π



0 FΓ

|

Xs

|δ

2

(

x

−

X

)

ds

.

We thus obtain the immersed boundary force F(s

,

t

)

=

FΓ

|

Xs

|

. 2.2. Nearly inextensible approach

As mentioned before, to avoid solving the extra unknown tension

σ

(

s

,

t

), we adopt an idea that is simpler than our

previ-ous work[15]to study 2D vesicle dynamics. Let us start with the derivation of inextensibility constraint(4)in axisymmetric coordinates. The rate of change of the surface dilating factor R

|

Xs

|

can be derived as follows[17].

∂

∂

t



R

|

Xs

|



=

∂

R

∂

t

|

Xs

| +

R

∂

|

Xs

|

∂

t

=

U

|

Xs

| +

R Xs

·

Xst

|

Xs

|

=

U

|

Xs

| +

R Xs

·

Us

|

Xs

|

=

U

|

Xs

| +

R

∂

U

∂

s

·

τ

=



U R

+

∂

U

∂

τ

·

τ



R

|

Xs

| = (∇

s

·

U

)

R

|

Xs

|.

(14)

The above derivation shows the expression of the surface divergence

∇

s

·

U explicitly, and also shows that if the surface is

locally inextensible (the rate of change of surface dilating factor equals to zero), then we obtain

∇

s

·

U

=

0 which is exactly

the constraint shown in Eq.(4). In order to relax the constraint a bit, one can replace the unknown tension by a spring-like elastic force as

σ

=

σ

0



R

|

Xs

| −

R0

|

Xs

|

0



(15) with a large stiffness constant

σ

0, where R0

|

Xs

|

0 is the initial surface dilating factor. Therefore, the interface now is nearly

inextensible. Consequently, the elastic energy Eσ is modified by

Eσ

(

t

)

=

π



0

σ

0 2



R

|

Xs

| −

R0

|

Xs

|

0



2 ds

,

(16)

and the total membrane energy becomes

E

(

t

)

=

π



0



σ

0 2



R

|

Xs

| −

R0

|

Xs

|

0



2

+

cb 2H 2_R

_|

_X s

|



ds

.

(17)

2.3. Elastic boundary force

Since the elastic energy is modified, the resultant tension force Fσ must be modified as well. This can be done by taking the variational derivative of new energy Eσ as follows. Let



X

(

s

)

= (

R

(

s

), 

Z

(

s

))

be a perturbation of the interface X, and

ε

be a small number. Then the perturbed energy becomes

Eσ

(

X

+

ε



X

)

=

π



0

σ

0 2

(

R

+

ε



R

)

|

Xs

+

ε



Xs

| −

R0

|

Xs

|

0

2 ds

.

Taking the derivative of above energy with respect to

ε

, we obtain dEσ dε

(

X

+

ε



X

)

=

π



0

σ

0

(

R

+

ε



R

)

|

Xs

+

ε



Xs

| −

R0

|

Xs

|

0



R

|

Xs

+

ε



Xs

| + (

R

+

ε



R

)

Xs

+

ε



Xs

|

Xs

+

ε



Xs

|

· 

Xs



ds

.

Then we evaluate the above equation at

ε

=

0 to obtain dE_σ dε

(

X

+

ε



X

)

_ε= 0

=

π



0

σ

0



R

|

Xs

| −

R0

|

Xs

|

0



R

|

Xs

| +

R

τ

· 

Xs



ds

=

π



0

σ



R

|

Xs

|

ds

+

π



0

σ

R

τ

· 

Xsds



since

σ

=

σ

0



R

|

Xs

| −

R0

|

Xs

|

0



=

π



0

σ



R

|

Xs

|

ds

−

π



0

∂

∂

s

(

R

σ τ

)

· 

X ds

(

integration by parts and R

=

0 at end points

)

=

 

σ

Re1



· 

X dS

−

 

1

|

Xs

|

∂

∂

s

(

σ τ

)

+

Rs

σ τ

R

|

Xs

|



· 

X dS



using dS

=

R

|

Xs

|

ds



=

 

_σ

R

|

Xs

|

2



R2_s

+

Z_s2

,

0



−

σ

R

|

Xs

|

2



R2_s

,

RsZs



−

1

|

Xs

|

∂

∂

s

(

σ τ

)



· 

X dS

=

 

σ

Zs R

|

Xs

|



Zs

|

Xs

|

,

−

Rs

|

Xs

|



−

1

|

Xs

|

∂

∂

s

(

σ τ

)



· 

X dS

= −

 

1

|

Xs

|

∂

∂

s

(

σ τ

)

−

Zs R

|

Xs

|

σ

n



· 

X dS

= −



Fσ

· 

X dS

.

This leads to the corresponding elastic force Fσ

=

1

|

Xs

|

∂

∂

s

(

σ τ

)

−

Zs R

|

Xs

|

σ

n

.

(18)

It is interesting to see that this formula for the elastic force is exactly same as the one in Eq.(13)despite the fact that the definitions of tension

σ

are different.

3. Numerical scheme

3.1. Fourier representation of the interface and computation of geometrical quantities

Before we proceed, we need to describe how we perform the spatial discretization for the interface. Unlike our previous works in [15,16]which used finite difference discretizations, here we adopt the spectral Fourier discretization to achieve higher-order of accuracy as in[35]. Since the interface is axisymmetric and is defined in 2D meridian domain (r–z plane with r

>

0), one can extend the interface representation X(s

)

= (

R

(

s

),

Z

(

s

)),

s

∈ [

0,

π

]

to be in another half of meridian domain (r

<

0). More precisely, we perform the radial component R

(

s

)

by the odd extension and the axial component Z(s

)

by the even extension so that this extending interfacial representation is periodic in

[

0,2

π

]

in the r–z plane. Under such extension, one can choose an even number of grid points Mpsuch that the interface is discretized by X(sk

)

= (

R

(

sk

),

Z

(

sk

)),

sk

=

k

s

=

2k

π

/

Mp, k

=

0,1, . . . ,Mp. (Note that, the spectral representation for the interface adopted in[35]used the sine

and cosine expansions for the radial and axial components, respectively.) Thus, the interface is represented by the discrete Fourier series expansion as

R

(

s

)

=

Mp

_

/2−1 l=−Mp/2



Rlei l s

,

and Z

(

s

)

=

Mp

_

/2−1 l=−Mp/2



Zlei l s

,

(19)

whereR



l and



Zl are the corresponding Fourier coefficients for R

(

s

)

and Z(s

), respectively. These Fourier coefficients can be

performed very efficiently by using the Fast Fourier Transform (FFT).

As described in previous section, in order to obtain the interfacial bending force, one needs to compute the mean curvature H , its surface Laplacian

sH , and the Gaussian curvature K . From the explicit formulas for these geometrical

quantities in Eqs. (10)–(11), one can see that the computations of them require the calculations of up to fourth-order derivatives of interface components R(s

)

and Z

(

s

). However, these can be done quite easily by using the pseudospectral

method[32]in which the Fourier coefficients of p-th derivative of R(s

)

can be obtained by taking an inverse FFT to

(

il

)

pR



l.

(Similar procedures should be done for the derivatives of Z

(

s

).) Notice that we should expand explicitly all the derivatives

of R(s

)

and Z(s

)

used in Eqs.(10)–(11)and compute them to enhance the accuracy.

Remark. It is important to mention that in the spectral computations for Eqs.(10)–(11), special care must be taken at the points of s

=

0 and s

=

π

where the radial component R

=

0. However, due to the odd extension of R(s

)

across at s

=

0,

π

, one can conclude that the second derivative Rssand the fourth derivative Rssssare both zero at those two points. Similarly,

due to the even extension of the axial component Z

(

s

), the first derivative Z

s and the third derivative Zsss are both zero

at s

=

0,

π

. By applying the l’Hospital’s rule and using those conditions mentioned above, one can evaluate the geometrical quantities at s

=

0,

π

as follows. lim s→0,πH

(

s

)

=

Zss

|

Rs

|

Rs

s=0,π

,

lim s→0,πHs

(

s

)

=

0

,

lim s→0,πHss

(

s

)

=

2R2_sZssss

−

6Zss

(

RsRsss

+

Z2ss

)

−

2RsRsssZss 3Rs

|

Rs

|

3

s=0,π

,

lim s→0,πK

(

s

)

=

Z_ss2 R4s

s=0,π

,

lim s→0,π

sH

(

s

)

=

2Hss R2s

s=0,π

.

3.2. Time-stepping scheme

Once we know how to represent the interface and to compute the geometrical quantities on the interface which are needed for the computation of bending force, we now are ready to introduce the time integration for the solution proce-dures. As mentioned before, we discretize the governing equations by the IB method. The fluid variables are defined on the staggered marker-and-cell (MAC) grid introduced by Harlow and Welsh [11]; that is, the pressure is defined on the grid points labeled as x

= (

ri

,

zj

)

= ((

i

−

1/2) r

,

c

+ (

j

−

1/2) z

), i

=

1,2, . . . ,M and j

=

1,2, . . . ,N, while the velocity

components u and w are defined at

(

ri

,

zj−1/2

)

= ((

i

−

1) r

,

c

+ (

j

−

1/2) z

)

and

(

ri−1/2

,

zj

)

= ((

i

−

1/2) r

,

c

+ (

j

−

1) z

),

as described in previous subsection, we use Fourier collocation to represent the interface and compute all the necessary ge-ometrical quantities from it. It is important to note that we only need to choose sufficient collocation points Mp to resolve

the interface.

Let us denote the time step size by

t and the superscript of the variables as the time step index. At the beginning of

each time step n, the interface configuration Xn

(

s

)

=

X

(

s

,

n

t

)

= (

Rn

(

s

),

Zn

(

s

))

is given so the unit tangent vector

τ

n

₍

_s

₎

₌

Xn

s

/

|

Xs

|

n can be evaluated accordingly. Here and the following, the subscript s-index denotes the spectral derivative with

respect to s as described in previous subsection. The unit outward normal nn

₍

_s

₎

_{can be obtained easily from}

_τ

n

₍

_s

_{). The}

steps of numerical time integration are as follows. 1. Compute the interfacial tension and bending force.

σ

n

=

σ

0



Rn

|

Xs

|

n

−

R0

|

Xs

|

0



,

Fn_σ

=

σ

_sn

τ

n

+

σ

n

τ

n_s

−

Z n s Rn

σ

n_nn

_,

Fn_b

=

cb



sHn

+

2Hn



Hn



2

−

Kn



|

Xs

|

nnn

,

Fn

=

Fn_σ

+

Fn_b

,

where the mean curvature H , its surface Laplacian

sH , and the Gaussian curvature K are all computed in the way

described in previous subsection. One should also notice that the term Zn

s

/

Rn in Fσ at the collocation points s

=

0,

π

must be modified to Zn

ss

/

Rns since R

=

0 at those two points.

2. Distribute the interfacial force from the interface to the fluids.

fn

(

x

)

=



s

Fn

(

s

)δ

2_h



x

−

Xn

(

s

)



s

,

where the discrete delta function is

δ

2_h

(

x

)

=

_h12

φ (

r h

)φ (

z h

)

[38]with

φ (

r

)

=

⎧

⎪

⎪

⎪

⎪

⎨

⎪

⎪

⎪

⎪

⎩

3 8

+

32π

−

r2 4

,

|

r

| <

0

.

5

,

1 4

+

1−|r| 8



−

2

+

8

|

r

| −

4r2

₋

1 8arcsin

(

√

2

(

|

r

| −

1

)),

0

.

5

 |

r

| 

1

.

5

,

17 16

−

π 64

−

3|r| 4

+

r2 8

+

| r|−2 16



−

14

+

16

|

r

| −

4r2

₊

1 16arcsin

(

√

2

(

|

r

| −

2

)),

1

.

5

 |

r

| 

2

.

5

,

0

,

2

.

5

 |

r

|.

Note that, when a marker point for the interface Xn _{is close to the z-axis (r}

₌

_{0), the force spreading region could}

cover some grid points in another half of meridian domain (r

<

0), which we should make the same supplement to the symmetric grid points in the right half domain.

3. Solve the Navier–Stokes equations. This can be done by the following semi-implicit second-order projection method in which the nonlinear term is approximated by the Adams–Bashforth scheme.

ρ

3u∗

−

4u n

₊

_un−1 2

t

+



2



un

· ∇

h



un

−



un−1

· ∇

h



un−1



= −∇

hpn

+

μ



hu∗

+

fn

,

h

φ

=

3 2

t

∇

h

·

u ∗

_,

∂φ

∂

n

=

0 on

∂Ω,

un+1

=

u∗

−

2 3

t

∇

h

φ,

∇

hpn+1

=

ρ

∇

h

φ

+ ∇

hpn

−

2 3

t

μ



h

(

∇

h

φ).

The spatial operators

∇

h and

∇

h

·

are the standard second-order central difference approximations to the gradient and

divergence operators in Eq. (6). Also, the spatial operator

h is the standard second-order central difference

approxi-mation to the Laplace operator in Eq.(7).

Note that, due to the axisymmetric property, the boundary conditions for velocity and pressure at r

=

0 are u

=

0,

∂

w

/∂

r

=

0, and

∂

p

/∂

r

=

0. One can see that the above Navier–Stokes solver requires to solve two Helmholtz-type equations for the velocity u∗ and one Poisson equation for the pressure increment

φ. These elliptic equations can be

solved efficiently using the public software FISHPACK[1].

4. Interpolate the new velocity from the fluid lattice points to the Lagrangian markers and then advance the markers to new positions as

Xn+1

(

s

)

=

Xn_k

+

t



x

Fig. 1. Four different surfaces: spherical, prolate, oscillatory and peanut-like surface (from left to right).

4. Numerical results

In this section, we perform a series of numerical tests for the present scheme developed in the previous section. We first verify the spectral accuracy for the spatial discretization of the interfaces by calculating their geometrical quantities such as mean curvature, Gaussian curvature, and the surface Laplacian of mean curvature. Then we choose various stiffness numbers and bending rigidities to study their effects on the vesicle dynamics. The mesh refinement study is performed to check the accuracy of the scheme, and an equilibrium bi-concave shape of vesicle in quiescent flow will be demonstrated. We also study the vesicle deformations under the gravitational force and conclude that, even with the same reduced volume but different initial shapes, the gravity can cause very different final stationary shapes of vesicles. Our numerical results are consistent with those observed in either 3D axisymmetric[35]or 3D simulations[4,36]. We conclude our numerical tests by simulating the vesicles in Poiseuille flows and study the equilibrium shapes under the effects of the different reduced volume, the confinement, and the mean flow velocity.

4.1. Accuracy check for the interfacial geometrical quantities

In the following, we check the accuracy of the spectral scheme described in Section 3.1 for the computations of the geometrical quantities of the interface under the Fourier spatial representations. Here, we test four different axisymmetric surfaces depicted inFig. 1in which

(

R

(

s

),

Z

(

s

))

are listed as follows.

•

Spherical surface:

(0.5 cos s,

0.5 sin s),

•

Prolate surface:

(0.1 cos s,

0.5 sin s),

•

Oscillatory surface[35]:

(

₂₀3r

(

s

)

cos s,₄₀3r

(

s

)

sin s), where r(s

)

=



sin2s

+

9 cos2_s

₊

_cos2_4s,

•

Peanut-like surface[40]:

((0.0414

+

0.4006 sin2s

−

0.2246 sin4s

)

cos s,0.5 sin s).

Table 1 shows the maximum absolute errors (not relative errors) between the numerical results obtained by present spectral method and the analytical results which can be computed easily from their parametric forms. Here, we show the errors for mean curvature H , Gaussian curvature K , and the surface Laplacian of mean curvature

sH as well. One can

see that, for smooth surface such as spherical, prolate or peanut-like surface, using the grid number Mp

=

32 to compute H and K has already achieved the accuracy close to the machine error. As the mesh is refined, one can see the absolute

errors slightly increase due to the round-off error effect which is quite common in practical spectral computations. This effect becomes even more prominent when we compute

sH since higher derivatives cause larger conditioned numbers.

The spectral method used in[35]shows the same error behaviors.

For the oscillatory surface, as expected, the finer grid Mp

=

256 results in a better accuracy than the coarser grid as in

other cases. Meanwhile, in this case, one can see the error decreases exponentially, which shows the spectral accuracy of the method. We also show the errors obtained by using the second-order finite difference scheme for the oscillatory surface case inTable 2. We can see clearly that the spectral method outperforms the finite difference method for the computations of the geometrical quantities in terms of the magnitude of absolute errors. Nevertheless, the finite difference method still shows the convergence in a lower order way.

4.2. Study on different stiffness number

σ

0

As mentioned before, the stiffness parameter

σ

0 in Eq.(15)must be chosen large enough to keep the nearly inextensible property. In this subsection, we test three different values of

σ

0to see how the vesicle behaves on those choices. We put a vesicle with elliptical shape X0

(

s

)

= (

0.5 cos s,0.1 sin s)in quiescent flow. The bending coefficient is fixed as cb

=

2

×

10−2

while three different stiffness parameters

σ

0

=

2

×

103

,

2

×

104, and 2

×

105are chosen. In the present method, we compute the immersed boundary force in the beginning of each step level so a restrictive time step constraint is needed to ensure numerical stability. Since the bending coefficient cb is small and the elastic stiffness

σ

0is large, the time step constraint is mainly associated with

σ

0. Here, we choose the time step based on the heuristic constraint

t

=

√Ch_σ

0, where C is a number of O(1). A similar choice criterion for the time step is also given in our previous work[15]despite the different meanings

Table 1

The mesh refinement study for the computations of H , K , and sH by the spectral method. The maximum absolute error is defined asH−He∞, where

Heis exact value of mean curvature. The errors of the other quantities are similarly defined.

Mp Spherical surface Prolate surface

H K sH H K sH

32 1.643E−14 6.528E−14 1.284E−11 3.339E−13 3.365E−11 1.909E−08 64 6.662E−13 2.646E−13 1.787E−10 7.176E−13 6.207E−11 8.812E−08 128 3.586E−13 1.434E−12 6.741E−09 8.647E−12 8.649E−10 3.949E−06

Mp Oscillatory surface Peanut-like surface

H K sH H K sH

32 4.600E−02 1.440E+00 1.907E+03 7.993E−14 1.705E−12 5.575E−10 64 3.336E−04 5.181E−03 2.799E+01 2.398E−13 4.874E−12 4.561E−09 128 3.572E−09 5.304E−08 9.222E−04 1.877E−12 3.969E−11 1.823E−07 256 6.068E−12 1.511E−10 1.475E−06

Table 2

The mesh refinement study for the computations of H , K and sH by the second-order finite difference method. Only the results for the case of the oscillatory surface are listed.

Mp H K sH

32 6.359E+00 1.703E+02 1.956E+05

64 3.061E+00 2.529E+02 1.197E+05

128 1.787E+00 4.163E+01 7.090E+04

256 5.464E−01 1.145E+01 3.499E+04

Table 3

The errors of the area dilating factor, the total surface area, and the volume.

σ0 t R|Xs| −R0|Xs|0∞ |Ah−A0|/A0 |Vh−V0|/V0

2×103 _{h/16 2.988E−04 2.431E−03 9.391E−04}

2×104 _{h/32 6.551E−05 2.060E−04 2.865E−04}

2×105 _{h/256 2.903E−05 2.105E−05 2.657E−04}

of penalty numbers. A rigorous energetic stability analysis for the present scheme is beyond the scope of this work and shall be more carefully studied in the future.

We run the simulations up to time T

=

2. As shown inFig. 2(a), the three vesicles coincide with each others (better view from cross-section) at the same times and they all reach a bi-concave shape eventually. We also plot the time evolution of total energy computed by Eq.(17)inFig. 2(b) which shows once again that the total energies coincide well. This consistence indicates that, as long as the stiffness parameter

σ

0is large enough, the vesicle motion and the flow remain the same.

Table 3shows the maximum errors of area dilating factor, and the relative errors of the total surface area and the volume of the vesicle at T

=

2. Notice that, the total surface area and the volume are calculated in the axisymmetric way. We can see that, as the stiffness increases, the maximum error of area dilation factor and the relative error of the surface area decrease. Interestingly, the order of increment of the stiffness is the same as the order of decrement of the surface area error. This kind of error behavior is physically reasonable, since, the larger

σ

0 is, the stronger spring force is provided to keep Rn

|

Xs

|

n closer to the initial area dilating factor R0

|

Xs

|

0. Although the relative error of the volume does not show the

same asymptotic behavior as the surface area error, both relative errors are within 0.03% when

σ

0 is equal to or greater than 2

×

104_{. We also test the case of a vesicle under a Poiseuille flow and obtain almost the same error behaviors so we} omit here. Based on this test, we choose

σ

0

=

2

×

104 in the remainder of this paper.

4.3. Study on different bending rigidity cb

In this subsection, we study the effect of different bending rigidity on the equilibrium shape of a freely suspended vesicle. Again, we put an initially elliptical vesicle X0

(

s

)

= (

0.5 cos s,0.1 sin s)in a quiescent flow. Now, we fix the stiffness number

σ

0

=

2

×

104 but vary the bending rigidity as cb

=

2

×

10−2, 2

×

10−3, and 2

×

10−6. All simulations are done

until steady equilibrium shapes have been achieved. As the bending rigidity becomes larger, the vesicle bends to bi-concave shape more significantly, seeFig. 3 in detail. The corresponding total energies are also plotted inFig. 4 which shows the vesicle with a larger bending rigidity naturally has a larger energy. Note that, in these computations, the bending energy is the dominated one comparing with the elastic energy.

Fig. 2. (a) A freely suspended vesicle with three different stiffness numbers:σ0=2×103(‘–’),σ0=2×104(‘×’), andσ0=2×105 (‘·’). (b) The corre-sponding time evolutions of total energy.

Fig. 3. Snapshots of the freely suspended vesicle with different bending rigidity. (a) cb=2×10−2; (b) cb=2×10−3; (c) cb=2×10−6. 4.4. Convergence test

In this subsection, we perform the convergence study for the present numerical algorithm. Again, we put a vesicle with the elliptical shape X0

(

s

)

= (

0.5 cos s,0.15 sin s) in quiescent flow initially. The computational domain is chosen as

[

0,1

] × [−

0.5,0.5

]

. We choose the different mesh sizes as m

=

64,128,256 and 512 so the corresponding mesh is

r

=

z

=

h

=

1/m. Since the analytical solution is not available in this test, we choose the result obtained from the finest

mesh m

=

512 as our reference solution and compute the maximum error between the reference solution and the other numerical solutions. The bending coefficient is chosen as cb

=

2

×

10−2, and all the numerical solutions are computed up to

time T

=

0.5.

Table 4shows the relative errors of the vesicle surface area and the vesicle volume and the maximum error of the vesicle interface configuration and the fluid velocity field. Since the fluid variables are defined at the staggered grid, when we refine the mesh, the numerical solutions on the refined mesh do not exist at the grid locations of the coarse mesh. In these runs, we simply use a cubic spline interpolation to compute the solutions at the necessary locations. Since the immersed

Fig. 4. The energy plots of the vesicle with different bending rigidity. cb=2×10−2(solid line); cb=2×10−3(dashed line); cb=2×10−6(dash-dotted line).

Table 4

The mesh refinement results for the surface area Ah, the enclosed volume Vh, the configuration Xh, and the velocity field uhand wh.

m=64 m=128 Rate m=256 Rate

|Ah−A0|/A0 4.032E−04 2.024E−04 0.99 1.009E−04 1.00

|Vh−V0|/V0 6.434E−04 1.505E−04 2.10 3.170E−05 2.25

Xh−Xref∞ 1.689E−03 4.625E−04 1.87 9.937E−05 2.22

uh−uref∞ 2.363E−03 1.162E−03 1.02 5.321E−04 1.13

wh−wref∞ 1.876E−03 7.652E−04 1.29 7.069E−04 0.11

boundary formulation has the singular forcing term in the equations, regularizing the singular term by smoothing discrete delta function causes the method to be first-order accurate. The numerical results shown in Table 4are consistent with what we expect from theory.

4.5. A suspended vesicle in quiescent flow

In this test, we put a freely suspended vesicle with an oscillatory initial configuration as described in Section 4.1 in a quiescent flow. The computational domain is the same as the previous tests, and the bending coefficient is chosen as

cb

=

5

×

10−2.Fig. 5shows the snapshots of the vesicle shapes in both cross-sectional and three-dimensional visualizations.

One can see that the initially oscillatory interface reaches a stationary bi-concave shape, which is consistent with the observations in[4,35]. The evolutional plot of the corresponding energy described in Eq.(17)is shown inFig. 6from which one can see that the total membrane energy is decreasing and tends to a constant value, indicating that the vesicle is in an equilibrium shape eventually.

4.6. Vesicles under the gravity

As in [4,35], we study the gravitational effect on the shape of the vesicle by considering different density across the interface. To simulate this problem, we simply add an extra interfacial force Fg as

Fg

=



ρ

i

−

ρ

o



(

g

·

X

)

|

Xs

|

n

,

into the vesicle boundary force [27]. Here,

ρ

i _and

_ρ

o _{represent the fluid densities inside and outside of the membrane,}

respectively, and g is the gravitational field. In the presence of the gravitation, a vesicle drops (since

ρ

i

_>

_ρ

o_{) and reaches}

a terminal velocity with an equilibrium shape.Fig. 7shows the snapshots of the vesicles with two different initial shapes; namely, an oblate shape X0

(

s

)

= (

0.25 cos s,0.125 sin s)and a prolate shape X0

(

s

)

= (

0.1531 cos s,0.3333 sin s). Both vesicles have the same volume and surface area. The oblate one (top) converges toward to a “parachute-like” shape while the prolate vesicle (bottom) converges toward to a “pear-like” shape. These results are qualitatively consistent with the results observed in[4,35].

Fig. 5. (a) Snapshots of a freely suspended vesicle in quiescent flow. (b) The corresponding 3D shapes of the vesicles.

Fig. 6. The time evolution of the membrane energy of the freely suspended vesicle.

Fig. 7. Snapshots of the vesicles under the gravity. The gravitational force is pointing into the negative z direction. Top: initially oblate shape; bottom:

Fig. 8. The velocity profile of Poiseuille flow. 4.7. Vesicles in Poiseuille flow

As mentioned before, the vesicles and RBCs share some similar mechanical behaviors; thus, the understanding of the former ones will lead to a better knowledge of the latter. The rheology of axisymmetric RBCs passing through capillaries in Poiseuille flows has been widely investigated numerically in[31,28]and the references therein. In[34], Veerapaneni et al. use the boundary integral method to study the dynamics of 3D axisymmetric vesicles suspended in Stokes flow. Recently, Coupier et al.[5]study the problem through experiments, numerical and theoretical computations to characterize the shape diagram of vesicles in Poiseuille flows. The authors reported that the most commonly stationary vesicle shapes observed in axisymmetric flow are bullet-like or parachute-like shapes. Here, we study the shapes of vesicles in 3D axisymmetric Poiseuille flow and compare the results with the ones obtained in[5].

In this subsection, we set up the Poiseuille flow as illustrated inFig. 8by

w

=

Wm



1

−

r 2 L2



,

u

=

0

,

where L is the capillary radius, and Wmis mean flow velocity indicating the centerline velocity. Throughout this subsection,

we draw the z-axis along the horizontal direction as inFig. 8(different from those in the previous subsections).

It is known that, in Poiseuille flow, a vesicle reaches its equilibrium shape and then moves at a constant velocity. As discussed in [5], the reduced volume

ν

, which measures the deflation of the vesicle, plays a significant role in the vesicle dynamics. It is defined as

ν

=

₄ V0

3

π

(

A0

/

4π

)

3/2

,

where V0 and A0 represent the volume and the surface area, respectively. This dimensionless number is nothing but the volume ratio of the vesicle to a sphere with the same surface area, and, thus, is equal to 1 for a sphere.

In addition, we define a characteristic non-dimensional parameter R to indicate the vesicle confinement in the flow by

ˆ

ˆ

R

=

R0 L

,

where the effective radius can be computed by R0

= (

3V0

/4

π

)

1/3. Thus, a larger R means a stronger confinement. In this

ˆ

subsection, we shall investigate three different effects by varying the reduced volume

ν

, the confinement R, and the mean

ˆ

velocity Wm individually. The computational r–z domain is chosen as

Ω

= [

0,L

] × [−

5L,5L

]

. Unless otherwise stated, we

use oblate vesicles as our initial shapes in the simulations.

•

Effect of the reduced volume. In this test, we choose an oblate vesicle with three different reduced volumes

ν

=

0.48,0.75,0.9 but the same volume in a Poiseuille flow with a weak confinement R

ˆ

=

0.3 and Wm

=

1. We run the

simulations until the equilibrium shapes are obtained.Fig. 9(a) shows the snapshots of the vesicles toward equilibrium shapes for the three different reduced volumes. One can see parachute-like shapes are observed in all three cases while the smaller reduced volume deforms significantly more since it is thinner. This behavior is consistent with the results obtained in [5]. On the other hand, it is interesting to see that, if we choose the prolate vesicle initially (with same surface area and volume as the oblate one with

ν

=

0.9), we obtain a bullet-like shape as shown inFig. 9(b). Notice that, the parachute-like shape has a concave rear part while the bullet-like shape has a convex one instead.

•

Effect of the confinement. To study the confinement effect, we keep the mean velocity Wm

=

1 and reduced volume

ν

=

0.95 fixed but vary the confinement as R

ˆ

=

0.3,0.375,0.5, respectively.Fig. 10shows the cross-sectional view of the equilibrium shapes with different confinements in which the confinement gets stronger from left to right. One can see that the equilibrium shape turns from a parachute-like shape into a bullet-like one as the confinement gets stronger. This is a physically interesting result which can be simply explained as the confinement is sufficiently large, its effect dominates the other flow effects. This is also in good agreement with the result obtained in[5].

Fig. 9. (a) Snapshots of the vesicles in Poiseuille flow with different reduced volume:ν=0.48 (top),ν=0.75 (middle), andν=0.9 (bottom). The flow comes from left to right. (b) A vesicle with initially prolate shape and reduced volumeν=0.9 results in a bullet-like shape eventually.

Fig. 10. Cross-sectional view of the equilibrium shapes with three different confinements: left:Rˆ=0.3; middle:Rˆ=0.375 and right: Rˆ=0.5. As the confinement increases from left to right, the shape turns from parachute-like to bullet-like.

•

Effect of the mean velocity. In the final test, we keep the reduced volume

ν

=

0.95 and a weak confinement R

ˆ

=

0.3 fixed but vary the mean velocity as Wm

=

1,10,100. Notice that, varying the mean velocity alone (but keeping the

other parameters fixed) means to vary the capillary number Ca

=

μR 4 0Wm

cbL2 . Thus, the effect of the mean velocity can

be regarded as the effect of capillary number.Fig. 11 shows the cross-sectional view of the equilibrium shapes with the different mean velocities. As the mean flow velocity increases, the equilibrium vesicle turns from a parachute-like shape to an unexpected bullet-like one. This interesting result has also been obtained in the experimental observations [5](the authors used

ν

between 0.95 and 0.97), and our result inFig. 11is qualitatively consistent with theirs.

Fig. 11. Cross-sectional view of the equilibrium shape with three different mean velocities: left: Wm=1; middle: Wm=10 and right: Wm=100. As the mean flow velocity increases, the shape turns from parachute-like to bullet-like.

5. Conclusions

In this paper, we have developed a simple numerical scheme to simulate 3D axisymmetric vesicle dynamics in Navier– Stokes flows. The mathematical model is based on the immersed boundary formulation in which the vesicle membrane forces are explicitly written as singular force terms in the equations. Instead of introducing a Lagrange’s multiplier (mem-brane tension) to enforce the vesicle inextensibility constraint, a spring-like tension is adopted to make the vesicle boundary nearly inextensible so that solving the unknown tension can be avoided. By choosing the spring stiffness large enough, we verify numerically that different choices of the stiffness do not affect the dynamics of the vesicles. Unlike the traditional IB method using finite difference approximations to discretize the interface (vesicle boundary), we instead use Fourier spectral approximations so that a higher-order of accuracy to compute the interfacially geometrical quantities and vesicle boundary forces can be achieved. A series of tests on the present scheme have been conducted to illustrate the applicability of the method. In particular, the shapes of vesicles in Poiseuille flow are investigated in detail to study the effects of the reduced volume, the confinement, and the mean flow velocity.

Our future work include developing an efficiently implicit or semi-implicit scheme so that the restriction on time step can be alleviated, and extending the present method to fully three dimensions.

Acknowledgements

This work has been finished during the time that M.-C. Lai was visiting the Research Institute for Mathematical Sci-ences (RIMS) at Kyoto University, Japan. He acknowledges the support and hospitality by the staff at RIMS and especially Professor Hisashi Okamoto. M.-C. Lai is also supported in part by National Science Council of Taiwan under research grant NSC-101-2115-M-009-014-MY3 and NCTS. Y. Kim was supported by National Research Foundation of Korea Grant funded by the Korean Government (2010-0006165).

Appendix A

In this appendix, we derive the mean curvature H , Gaussian curvature K , and the surface Laplacian of mean curvature

sH in axisymmetric coordinates as shown in Eq. (10). The mean curvature H and Gaussian curvature K can be

com-puted by the first and second fundamental forms [29]. As before, under axisymmetric assumption, the interface can be parameterized as

X

(

s

, θ )

=



R

(

s

)

cos

θ,

R

(

s

)

sin

θ,

Z

(

s

)



.

Here, we omit the dependence of time t for simplicity. The coefficients E, F , G of the first fundamental form for the above surface, and the surface area dilating factor W can be obtained as

E

=

Xs

·

Xs

= |

Xs

|

2

,

F

=

Xs

·

Xθ

=

0

,

G

=

Xθ

·

Xθ

=

R2

,

W

=



E G

−

F2

₌

_R

_|

_X

s

|.

The unit outward normal vector is

n

=

Xθ

×

Xs

W

=

1

|

Xs

|

(

Zscos

θ,

Zssin

θ,

−

Rs

).