key: cord-0045956-2f4mi8hb
authors: Zhu, Guangpu; Li, Aifen
title: Decoupled and Energy Stable Time-Marching Scheme for the Interfacial Flow with Soluble Surfactants
date: 2020-05-25
journal: Computational Science - ICCS 2020
DOI: 10.1007/978-3-030-50436-6_1
sha: 76d6928d8a8e2bc67386f475fda1d4a863a6cc8a
doc_id: 45956
cord_uid: 2f4mi8hb

In this work, we develop an efficient energy stable scheme for the hydrodynamics coupled phase-field surfactant model with variable densities. The thermodynamically consistent model consists of two Cahn–Hilliard–type equations and incompressible Navier–Stokes equation. We use two scalar auxiliary variables to transform nonlinear parts in the free energy functional into quadratic forms, and then they can be treated efficiently and semi-implicitly. A splitting method based on pressure stabilization is used to solve the Navier–Stokes equation. By some subtle explicit-implicit treatments to nonlinear convection and stress terms, we construct a first-order energy stable scheme for the two-phase system with soluble surfactants. The developed scheme is efficient and easy-to-implement. At each time step, computations of phase-field variables, the velocity and pressure are decoupled. We rigorously prove that the proposed scheme is unconditionally energy stable. Numerical results confirm that our scheme is accurate and energy stable.

Surfactants, interface active agents, are known to lower the interfacial tension and allow for the formation of emulsion [1, 2] . Commonly-used surfactants are amphiphilic compounds, meaning they contain both hydrophilic heads and hydrophobic tails [1, 3] . This special molecular composition enables surfactants to selectively absorb on fluid interfaces. Surfactants play a crucial role in everyday life and many industrial processes, such as the cleanser essence, the crude oil recovery and pharmaceutical materials, thus having an understanding of their behavior is a necessity. Numerical simulation is taking an increasingly significant position in investigating the interfacial phenomena, as it can provide easier access to some quantities such as surfactant concentration, pressure and velocity, which are difficult to measure experimentally. However, the computational modeling of interfacial dynamics with surfactants remains a challenging task. where F / ð Þ is the double well potential and G w ð Þ the logarithmic Flory-Huggins potential,

Two phase-field variables are used in the free energy functional. The first phase-field variable / uses two constants (-1 and 1) to distinguish two phases, and it varies continuously across the interface between -1 and 1. The other phase-field variable w is used to represent the surfactant concentration. The parameter Cn determines the interfacial thickness and Pi is a temperature-dependent parameter. More details of the free energy functional can refer to [6] and [19] . Although both the double well potential and the Flory-Huggins potential are bounded from below, the latter is not always positive in the whole domain. Thus, we add a zero term PiB À PiB to the free energy functional, and rewrite (2.1) into

where the positive constant B ensures G w ð Þþ B [ 0; and B = 1 is adopted in this study. Note that the free energy is not changed due to the introduction of the zero term PiB À PiB: We now use the scalar auxiliary variable (SAV) approach [12, 20] to transform the free functional into a new form. Through the simple substitution of scalar variables, the nonlinear parts of the free energy are transformed into quadratic forms of new scalar variables. More precisely, we define two scalar variables

Then the free energy can be transformed into

ð2:4Þ

Through the functional derivatives of E f with respect to phase-field variables / and w, we can obtain chemical potentials w / and w w

ð2:6Þ

Note that / 2 À 1 are denoted as Win (2.5) and (2.6).

Evolutions of phase-field variables / and w can be described by the conserved Cahn-Hilliard-type equations [14, 21] ,

where Pe / and Pe w are Péclet numbers. A degenerate mobility M w ¼ w 1 À w ð Þ; which vanishes at the extreme points w ¼ 0 and w ¼ 1; is adopted to combine with the logarithmic chemical potential w w : Eqs. (2.6)-(2.9) are coupled to the Navier-Stokes equation in the form [4, 14] qu t þ qu Á ru þ J Á ru À 1 Re

where DðuÞ ¼ ru þ r T u; and J ¼ k q À 1 À Á rw / 2Pe / : u is the velocity field, p is the pressure, Re is the Reynolds number and Ca is the Capillary number. We usually assume the density q and viscosity g has the following linear relations,

where k q and k g are density and viscosity ratios, respectively. In particular, if we consider the body force, e.g., the gravitational force, the dimensionless momentum equation read

where Bo ¼ ReCa is the Bond number, and g is the unit vector denoting the direction of body force. Periodic boundary conditions or the following boundary conditions

can be used to close the above governing system. Here C denotes boundaries of the domain. The total energy E tot of the hydrodynamic system (2.5)-(2.10) is the sum of kinetic energy E k and free energy E f

where W e ¼ ReCaCn, and we can easily derive the following energy dissipation law.

Next, we will develop an efficient time-marching scheme for the above governing system and carry out the energy estimate. To simplify the presentation, in the next section, we will take (2.9) as an example to construct the desired scheme.

We now present a first-order time-marching scheme to solve the governing system in Sect. 2. To deal with the case of nonmatching density, a cut-off function [4] is defined as

( Given w n ;/ n ; u n and p n , the scheme (3.1) calculates w n þ 1 ;/ n þ 1 ; u n+1 and p n+1 for n ! 0 in three steps.

In step 1, we update w n þ 1 and / n þ 1 by solving

ð3:1cÞ

ð3:1eÞ

with periodic boundary conditions or the following boundary conditions

where u n à is the intermediate velocity

In step 2, we update u n+1 by solving [22] q n u n þ 1 À u n dt þ q n u n

ð3:1iÞ

In step 3, we update p n+1 by solving the pressure Poisson equation with a constant coefficient [4, 23] 

where v ¼ 1 2 min 1 1 ; k q À Á :

; u n+1 and p n+1 are decoupled, which indicate that the scheme (3.1) is efficient and easy-to-implement. At each time step, u n+1 and p n+1 can be obtained by solving only two elliptic equations; Moreover, V n+1 and U n+1 do not involve any extra computational cost, since they can be calculated explicitly once we obtain w n þ 1 and / n þ 1 : (2) In the explicit convective velocity u n à ; we introduce a firstorder stabilization term [24] , which plays a dominant role in decoupling the computation of w n þ 1 ; / n þ 1 À Á from u n+1 and constructing the unconditionally energy stable scheme.

Theorem 3.1. The scheme (3.1) is unconditionally energy stable, and satisfies the following discrete energy dissipation law:

here Á k k denotes the L 2 -norm in X. Now we will rigorously prove the discrete energy dissipation law in (3.2). We first introduce an intermediate kinetic energy [25] as

The difference between E n þ 1 k and E n k;Ã is estimated as

We 2 q n þ 1 À q n ; u n þ 1 2 :

ð3:4Þ

Substituting (3.1i) into (3.1a), we obtain the following identity

We can also easily derive from (3.1g) that W eq n u n þ 1 À u n À Á þ dt w n rw n þ 1 w þ /rw n þ 1

Using the identity (3.6), we have W eq n u n þ 1 À u n

By taking the L 2 inner product of (3.7) with u n+1 , and using (3.4) and the following identities

ð3:8Þ

Using the Eq. (3.1g), we obtain

ð3:9Þ

Summing up Eqs. (3.8) and (3.9), we get

À W e 2 q n ; u n à À u n 2 :

By taking the L 2 inner product of (3.1j) with dt 2 We(p n+1 − 2p n + p n − 1)/v and withdt 2 Wep n+1 /v separately, we obtain À dt 2 W e 2v r p n þ 1 À p n À Á 2 À r p n À p nÀ1 À Á 2 þ r p n þ 1 À 2p n þ p nÀ1 À Á 2

¼ dtW e p n þ 1 À 2p n þ p nÀ1 ; r Á u n þ 1 À Á ;

ð3:11Þ and dt 2 W e 2v rp n þ 1 2 À rp n k k 2 þ r p n þ 1 À p n À Á 2 ¼ ÀdtW e p n þ 1 ; r Á u n þ 1 À Á : ð3:12Þ

Combining (3.11) and (3.12), yields À dtW e p n þ 1 À 2p n þ p nÀ1 ; r Á u n þ 1 À Á þ dtW e p n þ 1 ; r Á u n þ 1 À Á

ð3:13Þ

We take the difference of (3.1j) at step t n+1 and t n , pair the resulting equation with dt 2 We(p n+1 − 2p n + p n − 1)/(2v) then take integration by parts for both sides to derive dt 2 W e 2v r p n þ 1 À 2p n þ p nÀ1 À Á 2 v 2 W e u n þ 1 À u n 2 W e 4 q n ; u n þ 1 À u n À Á 2 : ð3:14Þ

Summing up Eqs. (3.10), (3.13) and (3.14), and using the triangle inequality W e 2 q n ; u n þ 1 À u n à 2 þ W e 2 q n ; u n à À u n 2 ! W e 4 q n ; u n þ 1 À u n À Á 2 :

ð3:15Þ

we can derive that

ð3:16Þ

By taking the inner product of (3.1a) with dtw n þ 1 w ; we can easily derive that

By taking the inner product of (3.1b) with À w n þ 1 À w n À Á ; we can derive that À w n þ 1 À w n ; w n þ 1 w ¼ À Pi V n þ 1 a n ; w n þ 1 À w n À Á À 1 4Ex / n j j 2 ; w n þ 1 À w n þ 1 4 W n j j 2 ; w n þ 1 À w n :

ð3:18Þ

where a n ¼ G 0 w n ð Þ ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi E v w n ð Þ p : Taking the inner product of (3.1c) with 2dtPiV n+1 to obtain

Summing up Eqs. (3.17)-(3.19), we get

/ n j j 2 ; w n þ 1 À w n þ 1 4 W n j j 2 ; w n þ 1 À w n :

By taking the inner product of (3.1d) with dtw n þ 1 / ; we have

By taking the inner product of (3.1e) with À / n þ 1 À / n À Á ; we can derive that À / n þ 1 À / n ; w n þ 1

w n þ 1 / n þ 1 ; / n þ 1 À / n À Á þ 1 2 w n þ 1 W n / n þ 1 þ / n À Á ; / n þ 1 À / n À Á ¼ À Cn 2 4 r/ n þ 1 2 À r/ n k k 2 þ r/ n þ 1 À r/ n 2 À U n þ 1 b n ; / n þ 1 À / n À Á À 1 4Ex w n þ 1 ; / n þ 1 2 À w n þ 1 ; / n j j 2 þ w n þ 1 ; / n þ 1 À / n 2 h i þ 1 4 w n þ 1 ; W n þ 1 2 À w n þ 1 ; W n j j 2 À w n þ 1 ; W n þ 1 À W n 2 h i :

: Taking the inner product of (3.1f) with 2dtU n+1 to obtain ð3:24Þ

Finally, combining (3.16) and (3.24) , we arrive at the desired result.

To implement the scheme (3.1), we use a finite difference method on staggered grids to discretize space. We pay special attention to the discretization of the convection terms in the Cahn-Hilliard and Navier-Stokes equations. A composite high resolution scheme, known as the MINMOD scheme, is used to reduce the undershoot and overshoot around the interface. The computations of w n þ 1 ;/ n þ 1 ; u n+1 and p n+1 can be totally decoupled if we replace w n þ 1 in (3.1e) with w n : The simplified scheme is extremely efficient and easy-to-implement. However, this simplification will definitely destroy the unconditional energy stability of our scheme. The implementation of such a simplified scheme requires small time step-sizes to obtain the desired accuracy and energy stability. The above scheme is adopted in [26] and numerical results demonstrate the energy stability of the proposed scheme. Here we will not present these results due to the limit of article length. We simulate the droplet deformation under the horizontal body force and a shear flow in a computational domain X = [0, 3] Â [0, 1]. Periodic boundary conditions are applied on the left and right sides. A circular droplet with the radius of r = 0.3 is initially placed at (1, 0.5). Other simulation parameters are listed as follows:

Pe / ¼ 10; Pe w ¼ 100; Re ¼ 10; Bo ¼ 1; Cn ¼ 0:01; Ex ¼ 1; Pi ¼ 0:1227; k q ¼ k v ¼ 10: Figure 1 shows the time evolution plots of droplet deformation and surfactant concentration. The droplet continuously deforms and moves forward under the action of the shear flow and the body force. We can divide the whole process into two stages based on the droplet deformation and surfactant migration. At the first stage, the body force has limited effect on the droplet deformation compared with the shear flow. Surfactants gradually migrate toward droplet tips, as shown in Fig. 1(b) , resulting in the non-uniformity of interfacial tension along the interface. As we mentioned before, the surfactant concentration gradient induces the Marangoni stress, which will resist the further migration of surfactants. However, the Marangoni stress is not large enough to balance the effect of shear flow, and surfactants continue to move toward tips. In Fig. 1 (c), surfactants are swept into the bulk phases when concentration reaches the maximum at the droplet tips. At the second stage, the body force plays an important role in the droplet deformation and surfactant migration. In Fig. 1(d) , surfactants on the tip A are slowly swept towards the ABC segment under the effect of the body force. Surfactants along the ADC segment continuously move to the tips under the combined action of the shear flow and the body force. Figure 2 demonstrates the profiles of phase-field variable / at three different w b values. A more prolate profile of / is observed for a higher surfactant bulk concentration, which confirms the effect of surfactants in reducing the interfacial tension.

The numerical approximation of incompressible and immiscible two-phase flows with soluble surfactants is the main topic in this paper. An efficient, accurate and energy stable time-marching scheme is constructed using the SAV approach for the hydrodynamics coupled phase-field surfactant model with variable densities. We rigorously prove the unconditional energy stability of the semi-implicit scheme. Numerical results demonstrate the energy stability of the proposed scheme.

An embedded boundary method for soluble surfactants with interface tracking for two-phase flows

Numerical approximations for the Cahn-Hilliard phase field model of the binary fluid-surfactant system

Phase-field modeling droplet dynamics with soluble surfactants

Decoupled, energy stable schemes for phase-field models of two-phase incompressible flows

A diffuse-interface method for simulating two-phase flows of complex fluids

A phase-field moving contact line model with soluble surfactants

Sharp-interface limits of a phase-field model with a generalized Navier slip boundary condition for moving contact lines

Efficient energy-stable schemes for the hydrodynamics coupled phase-field model

The effect of surfactants on the dynamics of phase separation

A comparison study of phase-field models for an immiscible binary mixture with surfactant

Linear and unconditionally energy stable schemes for the binary fluidsurfactant phase field model

Decoupled, energy stable schemes for a phase-field surfactant model

Diffuse interface model of surfactant adsorption onto flat and droplet interfaces

On diffuse interface modeling and simulation of surfactants in two-phase fluid flow

A hybrid numerical method for interfacial fluid flow with soluble surfactant

Numerical simulation of phase separation in the presence of surfactants and hydrodynamics

A diffuse-interface method for two-phase flows with soluble surfactants

A new phase-field model for a water-oil-surfactant system

Numerical approximation of a phase-field surfactant model with fluid flow

The scalar auxiliary variable (SAV) approach for gradient flows

Thermodynamically consistent modelling of two-phase flows with moving contact line and soluble surfactants

A novel energy stable numerical scheme for Navier-Stokes-Cahn-Hilliard two-phase flow model with variable densities and viscosities

An efficient scheme for a phase field model for the moving contact line problem with variable density and viscosity

Decoupled energy stable schemes for phase-field models of two-phase complex fluids

Thermodynamically consistent modeling and simulation of multicomponent two-phase flow with partial miscibility

Interfacial dynamics with soluble surfactants: a phase-field two-phase flow model with variable densities

Àdt r u n à w n À Á ; w n þ 1 w À dt r u n à / n