A short note on Navier-Stokes flows with an incompressible interface and its approximations

(1)

Contents lists available atScienceDirect

Applied Mathematics Letters

www.elsevier.com/locate/aml

A short note on Navier–Stokes flows with an incompressible interface and its approximations

Ming-Chih Lai

^a

, Yunchang Seol

^b,^*

aDepartment of Applied Mathematics, National Chiao Tung University, 1001, Ta Hsueh Road, Hsinchu 300, Taiwan

bNational Center for Theoretical Sciences, No. 1, Sec. 4, Road. Roosevelt, National Taiwan University, Taipei 10617, Taiwan

a r t i c l e i n f o

Article history:

Received 4 August 2016 Received in revised form 25 September 2016

Accepted 25 September 2016 Available online 1 October 2016

Keywords:

Incompressible interface Navier–Stokes flow

Immersed boundary method

a b s t r a c t

In biological applications, a cell membrane consisting of a lipid bilayer usually behaves as fluid-like interface with surface incompressibility. Here we consider a mathematical formulation for an incompressible interface immersed in Navier–

Stokes flows and study the mathematical and physical features for this incompressible interface. The model formulation introduces an unknown tension which acts as a Lagrange’s multiplier to enforce such surface incompressibility. In this note, we show that the spreading operator of the tension and the surface divergence operator of the velocity are skew-adjoint with each other which indicates physically that the tension does not do extra work to the fluid under the condition of surface incompressibility. In order to avoid solving the unknown tension to enforce the surface incompressibility, we adopt a nearly surface incompressible approach (or penalty approach) by introducing two different modified elastic tensions which can be used efficiently in practical numerical simulations. Furthermore, we show that the resultant modified elastic forces have the same mathematical form as the original one derived from the unknown tension.

1. Introduction

In cell biology, a vesicle is a small liquid droplet with a radius of about 10 µm enclosed by a lipid bilayer membrane suspended in a viscous incompressible fluid media. This membrane is negligibly thin with the thickness around 6 nm and exhibits resistance against membrane area dilation and bending. Therefore, it is natural to consider this fluid-like membrane as an incompressible interface Σ (t) suspended in a three- dimensional fluid domain Ω . Here, for simplicity, we disregard the membrane bending effect and focus on the tension effect due to surface incompressibility to the fluid. A similar formulation for the deformation of liquid capsules with incompressible interfaces in Stokes flow was considered in [1]. A complete mathematical

*Corresponding author.

E-mail addresses:mclai@math.nctu.edu.tw(M.-C. Lai),ycseol@ntu.edu.tw(Y. Seol).

(2)

formulation of vesicle modeling consisting of incompressible Navier–Stokes flows with bending and elastic forces acting on the vesicle membrane simultaneously can be found in [2]. Using the immersed boundary formulation [3], the Navier–Stokes flow interacting with an incompressible interface can be written as follows.

ρ( ∂u

∂t + (u · ∇)u )

= −∇p + µ∆u +

∫

Σ

Fσ(α, β, t) δ(x − X(α, β, t)) dA in Ω , (1)

∇ · u = 0 in Ω , (2)

F_σ(α, β, t) = (∇_sσ − 2σHn)(α, β, t) on Σ , (3)

∂X

∂t(α, β, t) = U(α, β, t) =

∫

Ω

u(x, t) δ(x − X(α, β, t)) dx on Σ , (4)

∇s· U = 0 on Σ . (5)

Here, we assume that the fluids inside and outside of the membrane have the same density ρ and viscosity µ. Eqs.(1) and (2)are the incompressible Navier–Stokes equations with the fluid velocity u(x, t) and the pressure p(x, t), where x = (x, y, z) denotes Cartesian coordinates in Ω and t is the time. The membrane is assumed as a sufficiently differentiable surface represented by Σ (t) = {X(α, β, t)|0 ≤ α ≤ ℓα, 0 ≤ β ≤ ℓ_β}, where α and β are Lagrangian coordinates of the reference configuration. The last term in Eq.(1) is the forcing term arising from the membrane elastic force Fσ and represents the spreading of Lagrangian force Fσ to the Eulerian fluid via the linkage of the Dirac delta function δ(x) = δ(x)δ(y)δ(z). Similarly, the membrane surface Σ (t) moves with the Lagrangian velocity U which is also interpolated using the local fluid velocity via the delta function in Eq.(4). Since the rate of change of the local stretching factor satisfies

∂

∂t(|X_α×X_β|) = (∇_s· U) |X_α×X_β| (6)

as shown in [2], the surface incompressibility indicates the above rate of change should equal to zero locally which results in the condition of Eq.(5).

The membrane force F_σwritten in Eq.(3)can be derived by taking the variational derivative of the elastic energy E_σ=∫

Σ (t)σ(α, β, t) dA with respect to the configuration in which σ is the variable membrane tension.

We shall have the similar derivation later on. Unlike the pure droplet problem, the surface tension σ in above formulation is not a physical quantity but is self-adjustable and acts as a Lagrange’s multiplier to enforce the surface incompressibility condition Eq.(5)which plays the same role as the pressure p in Navier–Stokes equations to enforce the fluid incompressibility of Eq.(2). In Eq.(3), H is the mean curvature and n is the outward unit normal vector of the membrane surface. The detailed expressions for the surface gradient ∇_sσ and the surface divergence ∇_s· U in terms of the coefficients of the first fundamental form will be given in next section.

The contribution of this paper is twofold. Firstly, we show that the spreading operator of the tension and the surface divergence operator of the velocity are skew-adjoint with each other which indicates physically that the tension does not do extra work to the fluid under the condition of surface incompressibility. Secondly, to avoid solving the unknown tension to enforce the surface incompressibility, we adopt a nearly surface incompressible approach (or penalty approach) which is shown quite useful in practical numerical simulations [2,4]. We introduce two other elastic penalty energies which are different from the one in our previous work [2,4] whose resultant elastic forces turn out to have the same mathematical form as in Eq. (3). Notice that, the present penalty energy is derived using classical differential geometry while the existing model such as Skalak formulation [5] is derived from strain energy in mechanics. Furthermore, in practice, the numerical scheme based on our present approach does not require to compute two stretching factors defined in orthogonal directions on the surface so the implementation of our approach using local surface areas is simper than other schemes.

(3)

2. Surface incompressibility in flows

In this section, we shall prove that the spreading operator of the tension and the surface divergence operator of the velocity are skew-adjoint mathematically. Before to proceed, we first review some preliminaries in classical differential geometry and express the surface gradient ∇_sσ and the surface divergence ∇_s · U in terms of the coefficients of the first fundamental form. We then derive alternate formulas for above quantities.

Let us denote two linearly independent tangent vectors on the surface by Xα = ^∂X_∂α and Xβ = ^∂X_∂β, respectively; and assume that n = (Xα×Xβ)/ |Xα×Xβ| is the unit normal vector pointing outward.

Throughout the paper, the subscripts of a function denote its partial derivatives. Using these notations, we define the coefficients of the first fundamental form by

E = X_α· X_α, F = X_α· X_β, G = X_β· X_β, and those of the second fundamental form by

L = −Xα· nα= Xαα· n, M = −Xα· nβ = −Xβ· nα= Xαβ· n, N = −Xβ· nβ= Xββ· n.

Thus, the local stretching factor can be written as |Xα×Xβ| =√

EG − F². The surface gradient operator [6] for a scalar function σ is expressed by

∇sσ = GXα− F Xβ

EG − F² σα+EXβ− F Xα

EG − F² σβ (7)

and similarly, the surface divergence operator for a vector function U is

∇s· U = GXα− F Xβ

EG − F² · Uα+EXβ− F Xα

EG − F² · Uβ. (8)

The mean curvature H = ¹₂∇s· n can also be written in terms of those coefficients as H = −GL−EN +2F M 2(EG−F²) . One can also rewrite Eqs. (7) and (8), without using the coefficients of first fundamental form. By substituting n = (Xα×Xβ)/ |Xα×Xβ| and using the vector triple product formulas a×(b×c) = (a · c)b − (a · b)c and (a×b)×c = (c · a)b − (c · b)a, we can easily obtain

Xβ×n = Xβ×Xα×Xβ

|Xα×Xβ| = (Xβ· Xβ)Xα− (Xβ· Xα)Xβ

|Xα×Xβ| = GXα− F Xβ

√EG − F² , n×X_α= Xα×Xβ

|Xα×Xβ|×X_α=(Xα· Xα)Xβ− (Xα· Xβ)Xα

|Xα×Xβ| = EXβ− F Xα

√EG − F² .

Therefore, the surface gradient of σ in Eq.(7)and the surface divergence of U in Eq.(8)can be respectively expressed by

∇_sσ = (Xβ×n)σα+ (n×Xα)σβ

|Xα×Xβ| , ∇_s· U = (Xβ×n) · Uα+ (n×Xα) · Uβ

|Xα×Xβ| . (9)

The following lemma shows the forms of mean curvature vector Hn, and the elastic force Fσ in Eq.(3).

Lemma 1.

(i) − 2Hn = Xβ×nα+ nβ×Xα

|Xα×Xβ| , (ii) ∇sσ − 2σHn = (σ(Xβ×n))_α+ (σ(n×Xα))_β

|Xα×Xβ| .

(4)

Proof . From n · n = 1, we have n_α· n = n_β· n = 0 meaning that n_α and n_β are both tangent vectors on the surface. So the above vectors can be written in the form of Weingarten equations [6] as nα =

−S₁₁X_α− S₂₁X_β and n_β= −S₁₂X_α− S₂₂X_β. Using these representations, we have X_β×n_α= S₁₁X_α×X_β and nβ×Xα= S22Xα×Xβ which lead to Xβ×nα+ nβ×Xα= (S11+ S22)Xα×Xβ. This actually completes the proof of (i), since H = −trace(S_ij)/2. By simply combining the identity in (i) and the surface gradient of ∇sσ in Eq.(9), one can easily derive the equality of (ii). □

In the following theorem, we show that the spreading operator of the tension and the surface divergence operator of the velocity are skew-adjoint with each other. To proceed, let us denote the inner product of two scalar functions on the surface Σ by ⟨f, g⟩_Σ =∫

Σf (X)g(X) dA and the inner product of two vector functions in the domain Ω by ⟨u, v⟩_Ω =∫

Ωu(x) · v(x) dx.

Theorem 1. Let us denote the spreading operator of the tension by S[F_σ] = ∫

Σ(∇_sσ − 2σHn)δ(x − X(α, β, t))dA. Then the rate of change of work done by the elastic force, ^dW_dt, should be equal to ^dW_dt =

⟨S[F_σ], u⟩_Ω = −⟨σ, ∇_s· U⟩_Σ. Proof .

dW

dt = ⟨S[F_σ], u⟩_Ω=

∫

Ω

S[F_σ] · u(x) dx

=

∫

Ω

[∫

Σ

(∇sσ − 2σHn)δ(x − X(α, β, t)) dA ]

· u(x) dx

=

∫ ∫

(∇sσ − 2σHn) · U(α, β, t) |Xα×Xβ| dα dβ (since U(α, β, t) = u(X(α, β, t), t))

=

∫ ∫ [

(σ(Xβ×n))_α+ (σ(n×Xα))_β]

· U(α, β, t) dα dβ (from Lemma 1(ii))

= −

∫ ∫

σ(X_β×n) · Uα+ σ(n×Xα) · Uβdα dβ (by integration by parts)

= −

∫ ∫

σ (∇s· U) |Xα×Xβ| dα dβ (by Eq.(9))

= −

∫

Σ

σ (∇s· U) dA = −⟨σ, ∇s· U⟩_Σ. □

It is important to mention that when the surface is incompressible ∇s· U = 0, the rate of change of work done by the elastic force equals to zero which has no surprise since the tension here acts like a Lagrange’s multiplier to enforce the surface incompressibility constraint. However, if the bending force Fb(the detailed form can be seen in [2]) is added in Eq.(3), then the rate of change of work done by the elastic membrane inTheorem 1should become ^dW_dt = ⟨S[Fσ] + S[Fb], u⟩_Ω = −⟨σ, ∇s· U⟩_Σ+ ⟨Fb, U⟩_Σ = ⟨Fb, U⟩_Σ whenever

∇s· U = 0. So the elastic tension force does not do extra work to the fluid, but the bending force does.

3. Formulations of nearly incompressible surface

As mentioned earlier, the tension σ(α, β, t) in Eq.(3)is an unknown variable and must be solved as part of solutions which complicates the numerical procedures in practice. To avoid solving the unknown tension to enforce the surface incompressibility, we adopt a nearly surface incompressible approach and introduce two different elastic penalty energies that results in different tensions. Nevertheless, those resultant elastic forces have the same mathematical form as Eq.(3)despite the fact that the associated elastic tensions are different as we show next.

(5)

Theorem 2. For the following elastic penalty energies E₁ and E₂ with the corresponding tensions σ₁ and σ2, their resultant elastic forces have the exactly same mathematical form as in Eq.(3). Here, the penalty constant σ0 is chosen to be sufficiently large and the superscript X⁰ denotes the initial configuration of the surface.

(i) Relative difference penalty energy

E1= σ0

2

∫ ∫ (

|Xα×Xβ|

|X⁰_α×X⁰_β|− 1 )²

|X⁰_α×X⁰_β| dα dβ, σ1= σ0

(|Xα×Xβ|

|X⁰_α×X⁰_β|− 1 )

.

(ii) Logarithmic difference penalty energy

E₂=σ₀ 2

∫ ∫ ( ln

(|X_α×X_β|

|X⁰_α×X⁰_β| ))²

|X⁰_α×X⁰_β| dα dβ, σ₂= σ₀|X⁰_α× X⁰_β|

|Xα× Xβ|ln

(|X_α× X_β|

|X⁰_α× X⁰_β| )

.

Proof . We perturb the configuration X by Y so the energy becomes

E1(X + εY) = σ0

2

∫ ∫ (

|(X + εY)α×(X + εY)β|

|X⁰_α×X⁰_β| − 1 )2

|X⁰_α×X⁰_β| dα dβ,

where ε is a small parameter. By taking the derivative with respect to ε and setting ε = 0, we obtain d

dεE1(X + εY)

⏐

⏐_ε=0

=

∫ ∫ σ0

(|(X + εY)α×(X + εY)β|

|X⁰_α×X⁰_β| − 1

) 1

|X⁰_α×X⁰_β|

× (X + εY)α×(X + εY)β

|(X + εY)α×(X + εY)β|· [Yα×(X + εY)β+ (X + εY)α×Yβ]

⏐

⏐_ε=0

|X⁰_α×X⁰_β| dα dβ

=

∫ ∫ σ0

(|Xα×Xβ|

|X⁰_α×X⁰_β|− 1

) Xα×Xβ

|Xα×Xβ| · (Yα×Xβ+ Xα×Yβ) dα dβ

=

∫ ∫

σ₁n · (Yα×Xβ+ Xα×Yβ) dα dβ (by the definition of σ₁in (i))

=

∫ ∫

σ1(Xβ×n) · Yα+ σ1(n×Xα) · Yβdα dβ (by the scalar triple product formula)

= −

∫ ∫

(σ1(Xβ×n))α· Y + (σ1(n×Xα))β· Y dα dβ (by integration by parts)

= −

∫ ∫

(∇_sσ₁− 2σ₁Hn) · Y|X_α×X_β| dα dβ (fromLemma 1(ii))

= −

∫

Σ (t)

(∇sσ1− 2σ1Hn) · Y dA (since dA = |Xα×Xβ| dα dβ)

=

∫

Σ (t)

δE1

δX · Y dA =

∫

Σ (t)

−Fσ· Y dA. (10)

According to the principle of virtual work, this concludes that the resultant force has the same mathematical form as shown in Eq.(3).

Similarly, the perturbed energy of E₂ can be written as

E₂(X + εY) = σ₀ 2

∫ ∫ (

ln|(X + εY)α×(X + εY)β|

|X⁰_α×X⁰_β|

)²

|X⁰_α×X⁰_β| dα dβ.

(6)

By taking the derivative with respect to ε and setting ε = 0, it yields d

dεE₂(X + εY)

⏐

⏐_ε=0

=

∫ ∫ σ₀

(

ln|(X + εY)α×(X + εY)β|

|X⁰_α×X⁰_β|

) 1

|(X + εY)α×(X + εY)β|

× (X + εY)α×(X + εY)β

|(X + εY)_α×(X + εY)_β|· [Yα×(X + εY)β+ (X + εY)α×Yβ]

⏐

⏐_ε=0

|X⁰_α×X⁰_β| dα dβ,

=

∫ ∫ σ0

|X⁰_α×X⁰_β|

|Xα×Xβ| (

ln|Xα×Xβ|

|X⁰_α×X⁰_β|

) Xα×Xβ

|Xα×Xβ|· (Yα×Xβ+ Xα×Yβ) dα dβ

=

∫ ∫

σ₂n · (Y_α×Xβ+ X_α×Yβ) dα dβ (by the definition of σ₂ in (ii)) .

The last equation has the same form as in Eq. (10)so the rest of derivation follows the same procedure before. □

We conclude this note by discussing how those penalty energies are motivated. As shown in Eq.(6), the exact surface incompressibility means that the local stretching factor does not change as time evolves; that is,

Instead of finding the unknown tension, one can choose the elastic tension with sufficiently large penalty constant σ0as in (i) to keep the ratio |Xα×Xβ| /|X⁰_α×X⁰_β| close to one. This is exactly how the elastic energy E₁ and the tension σ1 are derived. Alternatively, one can rewrite Eq. (6) as _∂t^∂ (ln |Xα×Xβ|) = ∇s· U.

Similarly, one can derive the elastic energy E₂ as in (ii) so the corresponding elastic tension becomes σ₂. Again, this nearly incompressible approach is to make logarithmic difference of the surface stretching factors such as ln(|Xα×Xβ| /|X⁰_α×X⁰_β|) = ln |Xα×Xβ|−ln |X⁰_α×X⁰_β| ≈ 0. A similar formulation of (i) is used by the present authors in [2]; however, the logarithmic difference penalty energy is new to the best of our knowledge.

One should, however, mention that Hencky model is based on the logarithmic behavior of one dimensional strain and extended to three dimensions providing a comparison with the experimental results of an elastic material such as vulcanized rubber [7]. By using the elastic tension σ₁or σ₂, we can circumvent the difficulty of finding the unknown tension in vesicle simulations to reduce computational complexity. Meanwhile, both elastic tension show similar effectiveness to keep the nearly incompressibility. Although not shown here due to limited space, we have implemented the above both approaches in our 3D code and the simulation results are in good agreements.

Acknowledgments

The work of M.-C. Lai is supported in part by Ministry of Science and Technology of Taiwan under research grant NSC-104-2115-M-009-014-MY3 and NCTS.

References

[1] H. Zhou, C. Pozrikidis, Deformation of liquid capsules with incompressible interfaces in simple shear flow, J. Fluid Mech.

283 (1995) 175–200.

[2] Y. Seol, W.-F.Hu. Y. Kim, M.-C. Lai, An immersed boundary method for simulating vesicle dynamics in three dimensions, J. Comput. Phys. 322 (2016) 125–141.

[3] C.S. Peskin, The Immersed Boundary Method, Acta Numer. 11 (2002) 479–517.

[4] W.-F. Kim, Y. Hu, M.-C. Lai, An immersed boundary method for simulating the dynamics of three-dimensional axisymmetric vesicles in Naiver-Stokes flows, J. Comput. Phys. 257 (2014) 670–686.

[5] R. Skalak, A. Tozeren, R.P. Zarda, S. Chien, Strain energy function of red blood cell membranes, Biophys. J. 13 (1973) 245–264.

[6] M. DoCarmo, Differential Geometry of Curves and Surfaces, Pearson, 2009.

[7] H. Hencky, The elastic behavior of vulcanized rubber, Rubber Chem. Technol. 6 (1933) 217–224.