• 沒有找到結果。

# 模擬帶有表面活性劑之軸對稱界面流問題

N/A
N/A
Protected

Share "模擬帶有表面活性劑之軸對稱界面流問題"

Copied!
40
0
0

(1)

## 學

(2)

### 碩 士 論 文

A Thesis

Submitted to Department of Applied Mathematics College of Science

National Chiao Tung University in partial Fulfillment of the Requirements

for the Degree of Master

in

Applied Mathematics

June 2009

Hsinchu, Taiwan, Republic of China

(3)

：賴明治

(4)

Yi-Min Huang

Ming-Chih Lai

### Department﹙Institute﹚of

Applied Mathematics

(5)

(6)

s

s

(7)

### 1 Introduction

The main purpose of this paper is to develop an immersed boundary method for Navier-Stokes equations in cylindrical geometries, and to use this solver to simulate the three-dimensional fluid interfaces with insoluble surfactant. We shall restrict ourselves in this paper to the axisymmetric case. Surfactant are surface active agents such as soap, deter-gent and food grease, which form a kind of skin on the surface of liquid and affect the interface surface tension. Effects of the surfactant are also important in bio-mechanical flows. For example, some insects drift on the water by injecting a chemical at the rear. This chemical reduces the surface tension behind their bodies so that the insects are pulled forward.

We refer to [5] and [6] in order to transform the immersed boundary formulation [7] in three-dimensional coordinates into the cylindrical coordinates. For our axisymmetric

case, we must set the boundary conditions on r = 0 such as ur = 0 and ∂uz/∂r = 0 [5].

In [7], the surfactant concentration equation in Lagrangian coordinates on the interface is presented, but it is in Cartesian coordinates. For testing a variety of axisymmetric cases, we should transform the equation in Cartesian coordinates into cylindrical coordinates according to [8].

Since our initial velocity field has to satisfy Navier-Stoke equations and the axisym-metric boundary conditions, the uniaxial extensional flow [9] is a good option. Surfactant on the interface drifts along the direction of the fluid flowing. Similarly, the markers in Lagrangian coordinates on the interface are affected by the flow. The following section will introduce a artificial method about marker redistribution.

For testing the case about moving contact line, put a drop on a solid on the slip bound-ary and control its equilibrium contact angle. Then, we impose the surfactant on the interface to observe that the drop slips with the time evolution. This simulation result is similar to that we use the detergent to remove grease. According to [12], we have Robin boundary condition on the slip boundary. We have to modify our solver for Navier-Stokes equations since this special boundary condition.

The remaining sections of this paper is organized as follows. In Section 2, we present the governing equations which includes the immersed boundary formulas and the sur-factant concentration equation in Lagrangian coordinates on the interface. Note that we impose artificial velocities on some equations to modify the distance between markers. The numerical method is described in Section 3 which includes an algorithm of solving the Navier-Stokes equations and a mass conservative scheme for the surfactant equation. We show our results about the effect of surfactant on drop deformation in some velocity fields in Section 4. In Section 5, we summarize some conclusions and our future work.

(8)

### 2 The governing equations

In a fixed three-dimensional cubic domain, we consider an incompressible two-phase flow problem and assume the interface is a simple closed surface immersed in the fluid

domain. The interface of the drop Σ separates two parts, Ω1 inside the interface and Ω2

outside the interface from our domain. i.e. Ω = [a, b] × [c, d] × [e, f ] = Ω1

S

Ω2. As

the interface is contaminated by surfactant, the distribution of the surfactant changes the surface tension accordingly. In each fluid region, the Navier-Stokes equations are satisfied as ρi ∂ui ∂t + (ui· ∇) ui ! = ∇ · Ti+ ρig, in Ωi, (2.1) ∇ · ui = 0, in Ωi, (2.2) u = ub, in ∂Ω, (2.3)

where for i = 1, 2 in each fluid domain, Ti = −piI + µi(∇ui + ∇uTi) is the stress tensor,

pi is the pressure, ui is the fluid velocity, ρi is the density, µi is the viscosity, and g is the

gravitational constant.

The velocity is continuous across the interface Σ. Thus,

[u]Σ = u|Σ,2− u|Σ,1 = 0 (2.4)

and the normal stress jump is balanced by the interfacial force F (defined only on Σ) as

[Tn]Σ+ F = 0 (2.5)

where n is the unit normal vector on Σ directed towards fluid 2. In this paper, we use the immersed boundary method to approach the velocity solutions of the Navier-Stokes equations. The immersed boundary ( the interface ) exerts force F to the fluids and moves with local fluid velocity. For computing easily, we first consider the case of equal viscosity

µ1 = µ2 = µ, equal density ρ1 = ρ2 = ρ, and neglect gravity. Hence, the interfacial force

term in delta function formulation and the surfactant concentration equation are the same problem in the single phase.

To be convenient to compute our numerical experiment in the single phase, we first non-dimensionlize the following three dimensional equations,

ρ ∂u ∂t + (u · ∇) u ! = −∇p + ∇ · (2µE) + f, (2.6) ∇ · u = 0, (2.7) ∂Γ ∂t + (∇s· u) Γ = Ds∇ 2 sΓ, (2.8)

where u is the velocity, p is the pressure, ρ is the density, µ is the viscosity, Γ is the

surfactant concentration, Dsis diffusivity, and E = (∇u+∇uT)/2 is the rate of deformation

(9)

Assume the dimensionless numbers, x = r0x, u = U∞u∗, t = (r0/U)t, p = ρU∞2 p∗,

f = f0f∗, and Γ = Γ∞Γ∗. Then, the non-dimensional equations are presented as follows.

∂u∗ ∂t∗ + (u ∗· ∇) u= −∇p+ 1 Re∇ ∗· (2E) + r0f0 ρU2 ∞ f∗, (2.9) ∇∗· u= 0, (2.10) ∂Γ∗ ∂t∗ + ∇ ∗ s· u∗  Γ∗ = 1 Pes ∇∗2 s Γ∗, (2.11)

where Re = ρUr0/µ and Pes = r0U/Ds. In the axisymmetric cylindrical coordinates,

we spread the interfacial force into fluid by Dirac delta function δ(x) = δ(r)δ(z) as follows. f (x, t) =

Z

Σ

∂στ

∂s δ (x − X (s, t)) ds,

where x is the fluid position, X is the interfacial position, σ is the surface tension, τ is the tangential vector on the interface, s is the interfacial parameter, and Σ is the domain of parameter.

Assume the dimensionless number σ = σ0σ∗, then we have

Z Σ ∂στ ∂s δ (x − X (s, t)) ds = σ0 r2 0 Z Σ ∂σ∗τ ∂s δ (x ∗− X) ds.

We suppose Ca = µU∞/σ0and use ∇ · u = 0, then

∂u∗ ∂t∗ + (u ∗· ∇) u= −∇p+ 1 Re∆ ∗(u) + 1 ReCa Z Σ ∂σ∗τ ∂s δ (x ∗− X) ds, (2.12) ∇∗· u= 0, (2.13) ∂Γ∗ ∂t∗ + ∇ ∗ s· u∗  Γ∗ = 1 Pes ∇∗2 s Γ∗. (2.14)

In fact, the characteristic velocity U= Gr0where G is the principle strain rate and 2r0is

the equivalent diameter of the bubble. Hence, we have Re = ρ(Gr0)r0/µ, Ca = µGr0/σ0,

and Pes = r0U/Ds. The following subsections will introduce the operator forms of

governing equations in the axisymmetric cylindrical coordinates.

### 2.1 Immersed boundary formulation

Throughout this paper, we consider the Navier-Stokes equations in cylindrical coordi-nates instead of Cartesian coordicoordi-nates, and the domain is axisymmetric. It is convenient to compute our data in a cylindrical coordinate system (r, θ, z). Due to the azimuthal sym-metry, the flow depends spatially on only two cylindrical coordinates(r, z). The interface

(10)

Σ is represented by a parametric form (r(s, t), z(s, t)), 0 ≤ s ≤ Lb, where s is the parameter

of the initial configuration of the interface, which is not necessarily the arc-length. By using the non-dimensionlization process, we can write down our governing equa-tions in the usual immersed boundary formulation as follows.

∂ur ∂t + ur ∂ur ∂r + uz ∂ur ∂z + ∂p ∂r = 1 Re ∂2u r ∂r2 + 1 r ∂ur ∂r + ∂2u r ∂z2 − ur r2 ! + 1 ReCafr, (2.15) ∂uz ∂t + ur ∂uz ∂r + uz ∂uz ∂z + ∂p ∂z = 1 Re ∂2u z ∂r2 + 1 r ∂uz ∂r + ∂2u z ∂z2 ! + 1 ReCafz, (2.16) 1 r ∂rur ∂r + ∂uz ∂z = 0, (2.17) f (x, t) = Z Σ F (s, t) δ (x − X (s, t)) ds, (2.18) ∂X (s, t) ∂t = U (s, t) = u (X (s, t) , t) = Z Ω u (x, t) δ (x − X (s, t)) dx, (2.19) F (s, t) =∂s(σ (s, t) τ (s, t)) , (2.20) τ (s, t) = ∂X ∂s ∂X ∂s , (2.21)

The dimensionless numbers are the Reynolds number (Re) describing the ratio between the inertial force and the viscous force, and the capillary number (Ca) describing the strength of the surface tension. The Eq.(2.18) represents the interfacial force F spreads to

the fluid such that we get the force f = ( fr, fz) on the fluid, and the force F on the interface

is balanced by the normal stress. Here, u = (ur, uz) is the velocity ( uris the radial velocity

and uzis the axial velocity ), σ is the surface tension, and τ is the unit tangent vector on

the interface. The Eq.(2.19) represents the interface moves with the fluid velocity. The Eulerian (x) and Lagrangian (X) variable are used in the three-dimensional Dirac delta function δ(x) = δ(r)δ(z).

The interfacial force F arises from the surface tension and its form is derived from Laplace-Young condition. If σ is not a constant, then we further take the derivatives in Eq.(2.20) as follows. F (s, t) =∂s(στ) = ∂σ ∂sτ + σ ∂τ ∂s = ∂σ ∂sτ + σκn ∂X∂s , (2.22)

where κ is the curvature of the interface and n is the unit outward normal. In Eq.(2.22), (∂σ/∂s) τ is the Marangoni force (the tangential force) and σκn |∂X/∂s| is the capillary force (the normal force). Note that if the surface tension is a constant, then the force only exerts in the normal direction. Let the interface be contaminated by the surfactant

(11)

which reduces the surface tension, we can find the phenomenon that the surface tension is smaller as the surfactant concentration is higher. The relation between surface tension and surfactant concentration can be described by the non-linear Langmuir equation. The following is the approximation of non-linear Langmuir equation.

σ (Γ) = σc(1 + ln (1 − βΓ)) , (2.23)

where Γ is the surfactant concentration, σcis the surface tension of a clean interface, and

β is a dimensionless number that measures the sensitivity of surface tension to changes in surfactant concentration.

Since we want to close the system, the surfactant is insoluble to the buck fluids such that it is convected and diffused along the interface. Hence, there is no exchange between the interface and the bulk fluids. By those reasons, the total mass of the surfactant must be conserved. The equation of surfactant concentration is derived in next subsection.

### 2.2 Surfactant concentration equation

In order to make the total mass of the surfactant conserved, we present a slightly different derivation from Stone [10] for the surfactant transport equation in the three-dimensional cylindrical coordinates.

In axisymmetric cylindrical coordinates, we can regard the surface as the curve in r-z coordinates. Thus, let L(t) be an interfacial segment where the surfactant concentration (the mass of the surfactant per unit surface area) is defined. Since the surfactant remain on the material element and do not transport or diffuse to the surrounding bulk fluids, the mass on the segment is conserved

d dt

Z

L(t)

Γ (l, t) rdl = 0, (2.24)

where dl is the arc-length element. To apply the time derivative more easily, we rewrite the above equation in terms of the initial parameter s as

d dt Z L(0) Γ (s, t) ∂X ∂s rds = 0. (2.25)

By taking the time derivative inside the integral, we obtain Z L(0) ∂Γ ∂t ∂X∂s r + Γ∂t∂ ∂X∂s r ! ! ds = 0. (2.26)

Note that, in our present formulation, both the interface and surfactant concentration are tracked in a Lagrangian manner. Thus, the time derivative of the first term in Eq.(2.26) is

(12)

exactly the material derivative of Stones derivation. Now we need to compute the rate of the stretching factor, and using Eq.(2.19), we have

∂t ∂X∂s r ! = ∂r ∂t ∂X∂s + r∂t∂ ∂X∂s = ur r r ∂X∂s + r ∂r ∂s∂s∂  ∂r ∂t  + ∂z ∂s∂s∂  ∂z ∂t  ∂X ∂s = ur r r ∂X∂s + r ∂r ∂s ∂ur ∂s + ∂z ∂s ∂uz ∂s ∂X ∂s (2.27) = ur r r ∂X∂s + r ∂r ∂s  ∇ur· ∂X∂s  + ∂z ∂s  ∇uz· ∂X∂s  ∂X ∂s = ur r r ∂X∂s + ∂u∂τ · τ ! r ∂X ∂s = (∇s· u) r ∂X∂s .

Here, the notation ∇s·u = ur/r+(∂u/∂τ)·τ means the surface divergence in axisymmetric

cylindrical coordinates. Since the material segment is arbitrary, we thus have r ∂X ∂s ∂Γ∂t + r ∂X∂s ( ∇s· u ) Γ = 0. (2.28)

If we allow surfactant diffusion along the interface, we obtain the surfactant transport-diffusion equation as r ∂X ∂s ∂Γ∂t + r ∂X∂s ( 5s· u ) Γ = 1 Pes∂s " r∂Γ ∂s ! , ∂X∂s # (2.29)

where Pesis the surface Peclet number and the boundary condition of Γ is ∂Γ/∂s = 0 on

∂L (t) . We note that surface diffusion is also written in terms of initial parameter s.

### 2.3 Controlling the interfacial markers uniformly

When there is the velocity field in the domain, the interface may be deformed ,and our computational markers on the interface concentrate on some part of interface with the direction of the fluid flowing. This phenomenon makes the larger error for the drop volume than the one which let the interfacial markers be distributed uniformly.

Since drop deformation is only related to the force in the normal direction, we impose an artificial velocity on Eq.(2.19) to redistribute the interfacial markers uniformly. Hence,

(13)

substitute Eq.(2.19) into ∂X (s, t) ∂t = Z Ω u (x, t) δ (x − X (s, t)) dx + UA(s, t) τ = u (X (s, t) , t) + UA(s, t) τ = U (s, t) + UA(s, t) τ, (2.30)

where the scalar UA is the artificial velocity, and it can modify the marker positions such

that they avoid being moved by the tangential force. This condition ∂ |∂X/∂s| /∂s = 0 means that the tangential displacement of markers remains unchanged. According to those conditions, we obtain

UA(s, t) − UA(0, t) = s 2π Z 0 ∂U ∂s0 · τ 0ds0 Z s 0 ∂U ∂s0 · τ 0ds0, (2.31)

where UA(0, t) = 0, since we should fix one marker to modify their relative positions.

After adding the artificial velocity term, the new marker positions have changed, so Eq.(2.29) has to add some term to preserve that the equation holds. The following is the modified surfactant transport-diffusion equation,

r ∂X ∂s ∂Γ∂t + r ∂X∂s ( 5s· u ) Γ − ∂ ∂s  rUAΓ= 1 Pes∂s " r∂Γ ∂s ! , ∂X∂s # . (2.32)

### 3 Numerical method

In this paper, the fluid flow variables are defined on a center of the cell which is divided by uniform grid mesh; that is, the pressure is defined on the grid points labelled as x = (ri−1/2, zj−1/2) = ((i − 1/2)h, ( j − 1/2)h) for i, j = 1, 2, . . . , N, the velocity components

ur and uz are defined at (ri−1/2, zj−1/2) = ((i − 1/2)h, ( j − 1/2)h), where the spacing h =

∆r = ∆z . For the immersed interface, we use a collection of discrete points sk = k∆s, k =

1, 2, . . . , M such that the Lagrangian markers are denoted by X = X(sk) = (rk, zk). The

surfactant concentration Γk and surface tension σk are defined at the half-integer points

given by sk−1/2 = (k − 1/2)∆s. Without loss of generality, for any function defined on the

interface φ(s), we approximate the partial derivative ∂φ/∂s by

Dsφ (s) =

φ (s + ∆s/2) − φ (s − ∆s/2)

∆s . (3.1)

Since |DsXk| can approximate the interface stretching factor by using this finite difference

convention, the unit tangent vector τkare defined at the half-integer points.

Let ∆t be the time step size, and n be the superscript time step index. At the beginning

of each time step, e.g., step n, the variables Xn

k = X (sk, n∆t) , Γnk = Γ sk−1/2, n∆t



, un =

u (x, n∆t), and pn = p (x, n∆t) are all given. The details of the numerical time integration

(14)

1. Compute the surface tension and unit tangent vector on the interface as σnk = σc 1 + ln 1 − βΓnk  , (3.2) τn k = DsXnk DsXnk , (3.3)

both of which hold for sk−1/2 = (k − 1/2)∆s. Then we define the interfacial force as

Fn

k = Ds σnkτnk



, (3.4)

at point Xk.

2. Distribute the force from the markers to the fluid by

fn(x) =X

k

Fnkδh x − Xnk



∆s, (3.5)

where the smooth version of Dirac delta function,

δh(x) = ( 1 4h  1 + cosπx 2h  , if −2h ≤ x ≤ 2h 0, otherwise, (3.6) is used.

3. Solve the Navier-Stokes equations. This can be done by the following semi-implicit second-order projection scheme, where the nonlinear term is approximated by a second-order extrapolation to avoid solving a nonlinear system at each time step.

3 ˜un+1− 4un+ un−1 2∆t +  2un· ˜∇ h  unun−1· ˜∇ h  un−1+ ˜∇pn = 1 Re˜∆˜u n+1+ 1 ReCaf n, (3.7) ˜u = ub, on ∂Ω, (3.8) ˜∇h· un+1= 0, (3.9) ˜∇2 hφn+1= 3 ˜∇h· ˜u 2∆t , ∂φ ∂n = 0, on ∂Ω (3.10) ˜∇hφn+1= 3˜un+1− un+1 2∆t , (3.11)

(15)

pn+1= pn+ φn+1 1

Re˜∇h· ˜u

n+1. (3.12)

The boundary conditions of the velocities are ur = 0 and ∂uz/∂r = 0 on r = 0. Here

the gradient operator ˜∇h = (∂/∂r, ∂/∂z) , the divergent operator ˜∇h· = (1r∂rr,∂z∂)· and

the Laplace operator ˜ ∆h˜un+1= ∂2un+1 r ∂r2 + 1 r ∂un+1 r ∂r + ∂2un+1 r ∂z2 − un+1 r r2 , ∂2un+1 z ∂r2 + 1 r ∂un+1 z ∂r + ∂2un+1 z ∂z2 ! ˜∇2 hφn+1 = ∂2φn+1 ∂r2 + 1 r ∂φn+1 ∂r + ∂2φn+1 ∂z2

are the standard centered difference operator on the grid which is defined at the center of the uniform grid mesh. One can see that the above Navier-Stokes solver

involves solving two Helmholtz equations for velocity ˜un+1and one Poisson

equa-tion for pressure. These elliptic equaequa-tions are solved by using the fast Poisson solver which is provided by the public software package Fishpack.

4. Interpolate the new velocity on the fluid lattice points onto the marker points and move the marker points to new positions.

Un+1 k = X x un+1δ h(x − Xnk)h2, (3.13) Xn+1k = Xnk + ∆tUn+1k . (3.14)

In order to redistribute the interfacial markers more uniformly, we should add the artificial velocity, UAn+1 k = k∆sM X k0=0 DsUn+1k · τn+1k ∆s − k X k0=0 DsUn+1k · τn+1k ∆s. (3.15)

From the above equations, we have the modified marker positions,

Xn+1

k = Xnk + ∆tUn+1k + ∆tUAn+1k τn+1k . (3.16)

5. Update surfactant concentration distribution Γn+1

k . Since the surfactant is insoluble,

the total mass on the interface must be conserved. Thus, it is important to develop a numerical scheme for the surfactant concentration equation to preserve the total mass. This can be done as follows. Substitute Eq.(2.27) of rate of stretching factor into the equation(2.29), we have

r ∂X ∂s ∂Γ∂t + ∂t∂ ∂X∂s r ! Γ = 1 Pes∂s " r∂Γ ∂s ! , ∂X∂s # (3.17)

(16)

Now we discretize the above equation by the Crank-Nicholson scheme in a sym-metric way as Γn+1 k − Γnk ∆t rn+1 k DsXn+1k + rnk DsXnk 2 + rn+1 k DsXn+1k − rnk DsXnk ∆t Γn+1 k + Γnk 2 = 1 Pes 1 ∆s    r n+1 k+1/2 |DsXn+1k+1| + |DsXn+1k | Γn+1 k+1− Γn+1k ∆srn+1 k−1/2 |DsXn+1k | + |DsXn+1k−1| Γn+1 k − Γn+1k−1 ∆s + r n k+1/2 DsXnk+1 + DsXnk Γn k+1− Γnk ∆srn k−1/2 DsXnk + DsXnk−1 Γn k− Γnk−1 ∆s    (3.18) Where k is “half-integer” index, i.e. k = 0.5, 1.5, . . . , M − 1/2. According to the

new interfacial marker position Xn+1

k obtained in the previous step and the boundary

condition (Γn+1

k − Γn+1k−1)/∆s = 0, where k = 0.5 or M + 1/2, on the endpoints

of the interface, Eq.(3.18) results in a symmetric tri-diagonal linear system which can be solved easily. More importantly, the total mass of surfactant is conserved numerically; that is,

X k Γn+1 k |DsXn+1k |rn+1k−1/2∆s = X k Γn k|DsXnk|rnk−1/2∆s (3.19)

The summation in Eq.(3.19) is the mid-point rule discretization for the integral in Eq.(2.25).The thing that the total mass is conserved can be derived by taking the summation of both sides of Eq.(3.18) and using the boundary condition.

Note that, when the artificial velocity is imposed on our equations like Eq.(3.16), Eq.(3.18) also have to impose the artificial velocity as follows.

Γn+1 k − Γnk ∆t rn+1 k DsXn+1k + rnk DsXnk 2 + rn+1 k DsXn+1k − rnk DsXnk ∆t Γn+1 k + Γnk 2 −   r n+1 k+1/2UAn+1k+1/2  Γn+1 k+1 + Γn+1k  /2 − rn+1 k−1/2UAn+1k−1/2  Γn+1 k + Γn+1k−1  /2 2∆s +r n k+1/2UAnk+1/2  Γn k+1+ Γnk  /2 − rn k−1/2UAnk−1/2  Γn k + Γnk−1  /2 2∆s    = 1 Pes 1 ∆s    r n+1 k+1/2 |DsXn+1k+1| + |DsXn+1k | Γn+1 k+1− Γn+1k ∆srn+1 k−1/2 |DsXn+1k | + |DsXn+1k−1| Γn+1 k − Γn+1k−1 ∆s + r n k+1/2 DsXnk+1 + DsXnk Γn k+1− Γnk ∆srn k−1/2 DsXnk + DsXnk−1 Γn k − Γnk−1 ∆s    .(3.20)

(17)

### 4.1 Oscillating drop

First, it is important to check how much the drop volume leaks or bulks up since it affects the computational accuracy. We test a simple case that a spheroid drop in the static liquid will be shrunk and turn into a sphere drop. Consider a computational domain Ω = [0, 1] × [−1, 1] in the axisymmetrical cylindrical coordinates, and put a spheroid drop with the semi-major axis whose length is a = 0.75 and with the semi-minor axis whose length is b = 1/3 into our domain. Then, given the computational mesh h = ∆r = ∆z, ∆t = h/8, Reynolds number Re = 25, and Capillary number Ca = 0.04 for our fluid.

0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 r−axis z−axis Time = 0 , (V a − Ve)/Ve = 0 % 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 r−axis z−axis Time = 0.25 , (V a − Ve)/Ve = 0.12587 % 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 r−axis z−axis Time = 0.5 , (V a − Ve)/Ve = 0.26542 % 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 r−axis z−axis Time = 0.75 , (V a − Ve)/Ve = 0.36258 % 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 r−axis z−axis Time = 1 , (Va − Ve)/Ve = 0.35525 % 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 r−axis z−axis Time = 1.25 , (Va − Ve)/Ve = 0.40839 % 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 r−axis z−axis Time = 1.5 , (Va − Ve)/Ve = 0.50401 % 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 r−axis z−axis Time = 1.75 , (Va − Ve)/Ve = 0.6027 %

(18)

0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 r−axis z−axis Time = 2 , (Va − Ve)/Ve = 0.68024 % 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 r−axis z−axis Time = 2.25 , (Va − Ve)/Ve = 0.75274 % 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 r−axis z−axis Time = 2.5 , (Va − Ve)/Ve = 0.83533 % 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 r−axis z−axis Time = 2.75 , (Va − Ve)/Ve = 0.92554 %

Figure 4.1: In the fluid velocity field, we observe the time evolution of a spheroid drop in a static liquid with Re = 25, Ca = 0.04 and ∆t = h/8.

In figure 4.1, the oscillation of the clean drop depends on the fluid velocity field and the incompressibility of the drop. The spheroid drop with a = 0.75 and b = 1/3 oscillates from the time T = 0 to T = 2, and its steady state is the sphere drop with the radius r  0.4368. After T = 2, the amplitude of vibration is smaller than the one before T = 2, but the drop volume leaks with the time evolution. The reason of this phenomenon is that there is the spurious velocities in our computational process. The spurious velocities influence the computation of the new marker positions on the interface .Thus, we only do our best to reduce the relative error of the drop volume. The relative error of the drop volume is smaller than 1% at T = 2.75 in figure 4.1. So, This is an acceptable error range for us.

### 4.2 Convergence test of fluid velocity and surfactant concentration

After above subsection, we impose the surfactant on the interface to observe the drop behavior. Throughout this paper, we choose β = 1 for the surface tension equation Eq.(3.2). In order to make our research convinced, we should carry out the convergence study of the immersed boundary method. We use the different mesh h = ∆r = ∆z = 1/32, 1/64, 1/128, 1/256 to perform our computation relatively. The Lagrangian mesh is cho-sen as ∆s ≈ h and the time step size is ∆t = h/8. The solutions are computed up to time T = 1. Since it is not easy to obtain the analytical solution from our problem, we choose

the finest mesh as our reference solution to compute the L∞ error between the reference

solution and the solution which is solved from the coarser grid. Here, we choose Re = 1,

(19)

Table 1: The mesh refinement analysis of the velocity ur, uz, and the surfactant

concen-tration Γ.

h kur− URre fk∞ Rate kuz− UZ re fk∞ Rate kΓ − Γre fk∞ Rate

1/32 8.5478E − 03 - 1.1470E − 02 - 1.2303E − 02

-1/64 2.4752E − 03 1.788 4.6027E − 03 1.317 2.4629E − 03 2.321

1/128 7.5121E − 04 1.720 1.8433E − 03 1.320 4.9265E − 04 2.322

Table 1 shows that the error of kur− URre fk∞, kuz− UZ re fk∞, and kΓ − Γre fk∞converge

to zero where URre f, UZ re f, and Γre f are r-axial, z-axial reference velocities, and reference

surfactant concentration relatively. The convergent rates for ur, uz, and Γ are about 1.7,

1.3, and 2.3 relatively, but our method is a second order scheme. Note that since the fluid velocities are defined at the center of the uniform grid in r − z coordinates, the refined mesh will not coincide with the same grid locations. Thus, we approximate the coarser grid by averaging the velocities obtained from the reference grid near the coarser grid. We attribute this is part of the reason why the rate of convergence behaves less than second-order. Similarly, the surfactant concentration is defined at “half-integer” grid. When we average the reference grid, we should use the weighted average since the concentration is the mass per unit surface area.

### 4.3 Clean vs. contaminated interface

There is different behavior for clean interface and contaminated interface, so we will study it in this subsection. We compare a clean drop with a contaminated drop in an

extensional flow, ur = −0.5Gr and uz = Gz, where G is the principal strain rate [5]. By

the comparison, we observe that the surface tension reduces when the surfactant is on the interface. Here, we use the mesh h = ∆x = ∆y = 1/64, the Lagrangian grid with size ∆s ≈ h, and the time step size ∆t = h/8.

4.3.1 The steady state of the clean interface in the extensional flow

In general, a drop will be stretched along the direction of fluid flow if it is put into the extensional flow. But there is surface tension on the interface, it may restrict a tendency of stretching drop. Therefore, before the comparison the clean interface with contaminated interface, we first observe a property for the clean drop. We fix a small Reynolds number, Re = 1, such that our Navier-Stokes equations are similar to the Stokes equations. Then, we attempt to change Capillary number and get the numerical result that the drop achieves a steady state.

(20)

0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis , z−right−point = 0.4 r−axis time =0 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis , z−right−point = 0.43685 r−axis time =0.25 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis , z−right−point = 0.45546 r−axis time =0.5 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis , z−right−point = 0.46548 r−axis time =0.75 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis , z−right−point = 0.47121 r−axis time =1 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis , z−right−point = 0.47455 r−axis time =1.25 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis , z−right−point = 0.47645 r−axis time =1.5 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis , z−right−point = 0.47748 r−axis time =1.75

Figure 4.2: The drop will form a steady state when the Capillary number is smaller than the critical value. Here, the dotted line is the initial interface and Ca = 0.2.

If Capillary number is 0.2 such as figure 4.2, then the drop with radius r = 0.4 is stretched slowly at about T = 0.5 ,and the tip points of the interface on the z-axis are about at the positions (r, z) ≈ (0, ±0.47). In figure 4.3, Capillary number is 0.27, and we find that this adjustment can not decelerate the stretching drop after T = 2. The reason is that Capillary number has a critical value in this case and it is between 0.2 and 0.27. Hence, the stretching drop will achieve a steady state quickly as the Capillary number is smaller than the critical value. But the stretching drop will not stop when the Capillary

(21)

to the viscosity µ, the strain rate G, the drop’s radius r, and the clean surface tension

σc. Hence, the larger strain rate and smaller surface tension can strengthen the stretching

phenomenon. 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis , z−right−point = 0.4 r−axis time =0 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis , z−right−point = 0.4949 r−axis time =0.5 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis , z−right−point = 0.54999 r−axis time =1 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis , z−right−point = 0.59339 r−axis time =1.5 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis , z−right−point = 0.63507 r−axis time =2 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis , z−right−point = 0.68175 r−axis time =2.5 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis , z−right−point = 0.74252 r−axis time =3 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis , z−right−point = 0.83688 r−axis time =3.5

Figure 4.3: If the Capillary number is larger than the critical value, then the drop continue to be stretched after T = 2. Here, the dotted line is the initial interface and Ca = 0.27.

4.3.2 The effect for the controlled markers on the interface

The initial markers on the interface are defined at the uniform grid, but the markers will not be distributed uniformly after the drop is put into the extensional flow with Re = 1 and Ca = 0.25. Most of interfacial markers concentrate in the vicinity of the tip points

(22)

of the drop as the fluid flows along the z-axial direction. For accurate computation, we impose an artificial velocity on Eq.(2.19) and Eq.(2.29) to redistribute interfacial marker positions uniformly. 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis r−axis time =0, (Va − Ve)/Ve = 0 % no surfactant Pes = 0.125 Pes = 12.5 Pes = 1250 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis r−axis time =0.5, (Va − Ve)/Ve = 0.028247 % no surfactant Pes = 0.125 Pes = 12.5 Pes = 1250 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis r−axis time =1, (Va − Ve)/Ve = 0.047636 % no surfactant Pes = 0.125 Pes = 12.5 Pes = 1250 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis r−axis time =1.5, (Va − Ve)/Ve = 0.10333 % no surfactant Pes = 0.125 Pes = 12.5 Pes = 1250 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis r−axis time =2, (Va − Ve)/Ve = 0.23924 % no surfactant Pes = 0.125 Pes = 12.5 Pes = 1250

Figure 4.4: In an extensional flow, compare the clean drop with the contaminated drops whose Peclet numbers are different. Note that the markers on interface are not distributed uniformly. 0 0.5 1 1.5 0 0.5 1 1.5 2 2.5 3 3.5 arc−length surfactant concentration time = 0 Pes = 0.125 Pes = 12.5 Pes = 1250 0 0.5 1 1.5 0 0.5 1 1.5 2 2.5 3 3.5 arc−length surfactant concentration time = 0.5 Pes = 0.125 Pes = 12.5 Pes = 1250 0 0.5 1 1.5 0 0.5 1 1.5 2 2.5 3 3.5 arc−length surfactant concentration time = 1 Pes = 0.125 Pes = 12.5 Pes = 1250 0 0.5 1 1.5 0 0.5 1 1.5 2 2.5 3 3.5 arc−length surfactant concentration time = 1.5 Pes = 0.125 Pes = 12.5 Pes = 1250 0 0.5 1 1.5 0 0.5 1 1.5 2 2.5 3 3.5 arc−length surfactant concentration time = 2 Pes = 0.125 Pes = 12.5 Pes = 1250

Figure 4.5: The time evolution of the surfactant dynamics in a extensional flow. Note that the markers on interface are not distributed uniformly.

The time and the relative volume error of drop with Pes = 12.5 is on the top of figure

4.4 and figure 4.6. In figure 4.4, the relative volume error of the drop is 0.23924% at T = 2, and it is larger than 0.024377% which is the relative volume error of the drop in figure 4.6. So, redistributing markers on the interface uniformly can reduce the volume error. Both in figure 4.4 and figure 4.6, the clean drop reaches a steady state as the Capillary number is smaller than the critical value.

(23)

0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis r−axis time =0, (Va − Ve)/Ve = 0 % no surfactant Pes = 0.125 Pes = 12.5 Pes = 1250 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis r−axis time =0.5, (Va − Ve)/Ve = 0.014915 % no surfactant Pes = 0.125 Pes = 12.5 Pes = 1250 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis r−axis time =1, (Va − Ve)/Ve = 0.038418 % no surfactant Pes = 0.125 Pes = 12.5 Pes = 1250 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis r−axis time =1.5, (Va − Ve)/Ve = 0.036789 % no surfactant Pes = 0.125 Pes = 12.5 Pes = 1250 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis r−axis time =2, (Va − Ve)/Ve = 0.024377 % no surfactant Pes = 0.125 Pes = 12.5 Pes = 1250

Figure 4.6: In an extensional flow, compare the clean drop with the contaminated drops whose Peclet numbers are different. Note that the markers on interface are distributed uniformly. 0 0.5 1 1.5 0 0.5 1 1.5 2 2.5 3 3.5 arc−length surfactant concentration time = 0 Pes = 0.125 Pes = 12.5 Pes = 1250 0 0.5 1 1.5 0 0.5 1 1.5 2 2.5 3 3.5 arc−length surfactant concentration time = 0.5 Pes = 0.125 Pes = 12.5 Pes = 1250 0 0.5 1 1.5 0 0.5 1 1.5 2 2.5 3 3.5 arc−length surfactant concentration time = 1 Pes = 0.125 Pes = 12.5 Pes = 1250 0 0.5 1 1.5 0 0.5 1 1.5 2 2.5 3 3.5 arc−length surfactant concentration time = 1.5 Pes = 0.125 Pes = 12.5 Pes = 1250 0 0.5 1 1.5 0 0.5 1 1.5 2 2.5 3 3.5 arc−length surfactant concentration time = 2 Pes = 0.125 Pes = 12.5 Pes = 1250

Figure 4.7: The time evolution of the surfactant dynamics in a extensional flow. Note that the markers on interface are distributed uniformly.

We put a drop into a uniaxial extensional flow such that the drop is stretched by the flow. In order to find the differences from the clean drop and the contaminated drop, we impose the surfactant on the clean interface and its surfactant concentration is Γ = 0.5.

There are three Peclet numbers, Pes = 0.125, Pes = 12.5, and Pes = 1250, in figure

4.6 and figure 4.7. The concentration curves has different shapes for different Peclet numbers. Whether we observe figure 4.4 or figure 4.6, we find that the interface move more slowly as the Peclet number is larger. The reason is that the high Peclet number reduce the effect of the surfactant diffusion. According to the direction of the fluid flow, the surfactant concentration in the vicinity of the tip points of the drop is higher than others. The evidence is in figure 4.5 and figure 4.7.

(24)

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0.392 0.3921 0.3922 0.3923 0.3924 0.3925 0.3926 0.3927 0.3928 0.3929 0.393 time

total mass of surfactant

Figure 4.8: This is the total mass of surfactant for Pes = 0.125, Pes = 12.5, and Pes =

1250.

The diffusion term in Eq.(2.29) influences that the surfactant diffuses from the high concentration to the low concentration. Therefore, the surfactant is carried by the fluid flowing, but the diffusion effect is very small if Peclet number is very large. In figure 4.7, it is more obvious that the surfactant concentrate in the vicinity of the interfacial endpoints as Peclet number is larger. Let’s attempt to explain why the interface move more slowly as the Peclet number is larger. It is known for us that the surfactant reduce the surface tension on the interface. Smaller Peclet number makes the surfactant distributed uniformly, so the surface tension is reduced averagely. But larger Peclet number makes the surfactant concentrate in the vicinity of the interfacial endpoints. Most of surface tension on the interface only reduce a little, but the higher concentration near the interfacial endpoints makes the tip points of the drop become sharper. By the comparison, the interface move more slowly as the Peclet number is larger. In figure 4.8, no matter which Peclet number we choose, the total mass of surfactant is conserved very well by our method.

### 4.4 Forming the dumbbell-shaped drop in another fluid velocity field

In order to ensure our scheme and concept, we want to observe whether there is a similar numerical result for the drop in another velocity field. Consider the velocity field

, ur = −0.5Gπr cos(πz) and uz = G sin(πz), which satisfies Eq.(2.17) and the boundary

condition ur = 0, ∂uz/∂r = 0 on r = 0 in the axisymmetric cylindrical coordinates. Put a

drop into the velocity field with Re = 1 and Ca = 0.5. Then we impose the surfactant on the interface. In the following result, we have imposed the artificial velocity to uniformly redistribute the marker positions.

(25)

0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis r−axis time = 0, (Va − Ve)/Ve = 0 % 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis r−axis time = 3, (Va − Ve)/Ve = 0.38462 % 0 0.5 1 1.5 2 2.5 3 1.5706 1.5706 1.5706 1.5706 1.5706 1.5706 1.5706 1.5706 1.5706 1.5706 time

total mass of surfactant

Figure 4.9: The left two figures are the process that the drop is stretched from T = 0 to T = 3. The total mass of surfactant is on the right hand side.

Observing the left two subfigures in figure 4.9, the drop with Γ = 0.5 and Pes = 50

is stretched along z = −1 and z = 1 in the velocity field. At T = 3, the surfactant concentration is higher near the stretched parts than other parts of the interface. The reason is that the surfactant is carried along the direction of the fluid flowing. We can compare the drops with different Peclet numbers in figure 4.10. The stretched part of the

drop with Pes = 50 is sharper than the one with Pes = 0.5 since there is more surfactant

concentrating on the tip parts. In addition, the clean drop is more smooth than others.

0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis r−axis time = 0 no surfactant Pes = 0.5 Pes = 50 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis r−axis time = 1 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis r−axis time = 2 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis r−axis time = 3

(26)

0 1 2 3 0 1 2 3 4 5 6 arc−length surfactant concentration time = 0 Pes = 0.5 Pes = 50 0 1 2 3 0 1 2 3 4 5 6 arc−length surfactant concentration time = 1 Pes = 0.5 Pes = 50 0 1 2 3 0 1 2 3 4 5 6 arc−length surfactant concentration time = 2 Pes = 0.5 Pes = 50 0 1 2 3 0 1 2 3 4 5 6 arc−length surfactant concentration time = 3 Pes = 0.5 Pes = 50

Figure 4.11: The time evolution of the surfactant dynamics in this velocity field.

For the curve of surfactant concentration with Pes = 0.5 in figure 4.11, it is more

smooth and more straight than the one with Pes = 50 since larger Peclet number reduces

the diffusion term in Eq.(2.29). On the right hand side of figure 4.9, the total mass of surfactant is conserved. Hence, this result is similar to the one about the drop in the extensional flow.

### 4.5 Moving contact line on the slip boundary

In this case, we put a clean drop onto a solid, and the initial contact angle is a right angle. We restrict the computational domain Ω = [0, 1] × [0, 1] in r − z coordinates. In axisymmetric geometries, most of boundary conditions are the same as the previous cases

except the slip boundary conditions uz = 0 and ur = (∂ur/∂z) on z = 0, where  is a

coefficient about the slip length. Instead of the solver for the software package Fishpack, we solve Eq.(2.15) by the IMSL subroutine “DLSLXD” since the boundary conditions on

z = 0 have changed. According to Young’s condition σ cos θssl = σsg, we get the force

Fc = σc(cos θs− cos θd) at moving contact line along r-axis where θsis the equilibrium

contact angle and θdis the dynamic contact angle.

4.5.1 Convergence test of fluid velocity and surfactant concentration

First, we should carry out the convergence study of the immersed boundary method. We use the different mesh h = ∆r = ∆z = 1/32, 1/64, 1/128, 1/256 to perform our computation relatively. The Lagrangian mesh is chosen as ∆s ≈ h/2 and the time step size is ∆t = h/8. The solutions are computed up to time T = 1. Since it is not easy to obtain the

(27)

analytical solution from our problem, we choose the finest mesh as our reference solution

to compute the L∞ error between the reference solution and the solution which is solved

from the coarser grid. Here, we choose Re = 1, Ca = 0.5, Pes = 50, θs = 60◦, and  = h

for our computation.

Table 2: The mesh refinement analysis of the velocity ur, uz, and the surfactant

concen-tration Γ.

h kur− URre fk∞ Rate kuz− UZ re fk∞ Rate kΓ − Γre fk∞ Rate

1/32 1.0782E − 03 - 2.9943E − 03 - 2.7585E − 03

-1/64 4.3824E − 04 1.299 1.3635E − 03 1.135 1.0808E − 03 1.352

1/128 1.5433E − 04 1.506 5.0570E − 04 1.431 3.7413E − 04 1.531

Table 2 shows that the error of kur− URre fk∞, kuz− UZ re fk∞, and kΓ − Γre fk∞converge

to zero where URre f, UZ re f, and Γre f are r-axial, z-axial reference velocities, and reference

surfactant concentration relatively. The convergent rates for ur, uz, and Γ are about 1.5,

1.4, and 1.5 relatively, but our method is a second order scheme. The reason is that it produces some errors as the data location which we want are interpolated by the reference data.

4.5.2 Clean vs. contaminated interface with θs = 60◦

In this subsection, we want to compare the clean interface with contaminated

inter-face as the equilibrium contact angle is 60◦. Since the surfactant is insoluble and only

distributed on the liquid-to-liquid interface, σsl and σsgare not affected by the surfactant.

Therefore, the force Fc at the moving contact line must be different from the one on the

clean interface. That is Fc = σsg− σsl− σ cos θd where σsg− σsl = σccos θs.

In order to achieve the equilibrium contact angle quickly, we set  = h = 1/128,

∆t = h/8, Re = 100, and Weber number We = ReCa = 0.1. We set Pes = 500 such that

there is the obvious fluctuation on the surfactant concentration curve. In figure 4.12, we find that the tip points of interfaces on z-axis move slowly after about T = 4. The tip point of contaminated interface on z-axis is farther from the original point than The tip point

of clean interface on z-axis. It is imaginable since the term σ cos θd in the force on the

moving contact line is smaller than σccos θs. Hence, θd must be so small (θd < θs) that it

can achieve the equilibrium and σccos θs = σ cos θd.

The fluctuation of concentration curve is presented in figure 4.13. We can compare this figure with figure 4.12. Note that the right endpoint of the curve is defined at the point which is closest to the moving contact line. In the beginning, the surfactant concentrate in the vicinity of the moving contact line since the slipping drop produce the similar velocity field as the extensional flow out along r-axis. After T = 4, the concentration curve becomes straight gradually since the drop becomes a steady state gradually.

(28)

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z−axis r−axis time = 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z−axis r−axis time = 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z−axis r−axis time = 2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z−axis r−axis time = 4 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z−axis r−axis time = 6 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z−axis r−axis time = 8

Figure 4.12: The drop is on the slip boundary. Clean drop (solid line) and contaminated drop (‘− · − · −·’) are in this figure.

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.2 0.22 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 Maker arc−length surfactant concentration time = 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.2 0.22 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 Maker arc−length surfactant concentration time = 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.2 0.22 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 Maker arc−length surfactant concentration time = 2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.2 0.22 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 Maker arc−length surfactant concentration time = 4 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.2 0.22 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 Maker arc−length surfactant concentration time = 6 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.2 0.22 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 Maker arc−length surfactant concentration time = 8

(29)

0 1 2 3 4 5 6 7 8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8

Relative Error of Volume : %

Time 0 1 2 3 4 5 6 7 8 0.4707 0.4708 0.4709 0.471 0.4711 0.4712 0.4713 0.4714 0.4715 0.4716 0.4717 Surfactant Mass Time 0 1 2 3 4 5 6 7 8 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8

Cosine Value of the Dynamic Contact Angle

Time 0 1 2 3 4 5 6 7 8 40 45 50 55 60 65 70 75 80 85 90

Dynamic Contact Angle

Time

Figure 4.14: This figure presents the relative volume error, the total mass of surfactant,

the variation of cos θd, and the variation of θd with time.(solid line : clean drop, ‘− · − · −·’

: contaminated drop)

In figure 4.14, the relative volume error is smaller than 1% from T = 0 to T = 8. The total mass of surfactant is conserved by our method. For the clean drop, the contact angle

is about 60◦after T = 4, but the contact angle is about 45for the contaminated drop since

the surfactant has changed the original equilibrium contact angle. In fact, our Reynolds number is so large that it influences the smoothness for the curves about concentration and contact angle.

4.5.3 Clean vs. contaminated interface with θs = 120◦

(30)

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z−axis r−axis time = 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z−axis r−axis time = 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z−axis r−axis time = 3 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z−axis r−axis time = 6 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z−axis r−axis time = 9 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z−axis r−axis time = 12

Figure 4.15: The drop is on the slip boundary. Clean drop (solid line) and contaminated drop (‘− · − · −·’) are in this figure.

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.2 0.22 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 Maker arc−length surfactant concentration time = 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.2 0.22 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 Maker arc−length surfactant concentration time = 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.2 0.22 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 Maker arc−length surfactant concentration time = 3 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.2 0.22 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 Maker arc−length surfactant concentration time = 6 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.2 0.22 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 Maker arc−length surfactant concentration time = 9 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.2 0.22 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 Maker arc−length surfactant concentration time = 12

(31)

0 2 4 6 8 10 12 −1 0 1 2 3 4 5 6

Relative Error of Volume : %

Time 0 1 2 3 4 5 6 7 8 0.4707 0.4708 0.4709 0.471 0.4711 0.4712 0.4713 0.4714 0.4715 0.4716 0.4717 Surfactant Mass Time 0 2 4 6 8 10 12 −0.8 −0.7 −0.6 −0.5 −0.4 −0.3 −0.2 −0.1 0 0.1

Cosine Value of the Dynamic Contact Angle

Time 0 2 4 6 8 10 12 80 90 100 110 120 130 140

Dynamic Contact Angle

Time

Figure 4.17: This figure presents the relative volume error, the total mass of surfactant,

the variation of cos θd, and the variation of θd with time.(solid line : clean drop, ‘− · − · −·’

: contaminated drop)

In figure 4.15, we find that the tip points of interfaces on z-axis move slowly after about T = 6. The tip point of contaminated interface on z-axis is closer to the original

point than The tip point of clean interface on z-axis. In this case, the term σ cos θd in the

force on the moving contact line must be equal to σccos θs to achieve the steady state.

Hence, the new equilibrium contact angle is larger than original contact angle after the interface is contaminated.

The fluctuation of concentration curve is presented in figure 4.16. We can compare this figure with figure 4.15. Note that the right endpoint of the curve is defined at the point which is closest to the moving contact line. In the beginning, the surfactant concentrate in the vicinity of the moving contact line since the slipping drop produce the similar

(32)

velocity field as the extensional flow out along z-axis. After T = 6, the concentration curve becomes straight gradually since the drop becomes a steady state gradually.

In figure 4.17, the relative volume error is larger than 1% after T = 2, and it is over 5% at T = 12. This result is a little bad for us. The reason may be that in order to achieve the equilibrium contact angle quickly, larger Reynolds number influences the accuracy. The total mass of surfactant is conserved by our method. For the clean drop, the contact

angle is about 120◦after T = 4, but the contact angle is about 135for the contaminated

drop since the surfactant has changed the original equilibrium contact angle. In fact, our Reynolds number is so large that it influences the smoothness for the curves about concentration and contact angle.

4.5.4 Imposing initial velocity on our moving contact line case

In this subsection, we impose initial velocity, ur = Grz2and uz= −2Gz3/3, on our case.

This velocity satisfies the slip boundary condition, axisymmetric boundary condition, and ∇ · u = 0. In order to reduce the effect of the velocity field to make contact line moving,

we choose G = 0.5 and the equilibrium contact angle θs = 0◦. Given initial contact angle

120◦, Re = 10, We = 0.1 and Pe

s= 625 for the drop and the velocity field. Compare the

clean drop with the contaminated drop as follows.

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z−axis r−axis time = 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z−axis r−axis time = 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z−axis r−axis time = 3 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z−axis r−axis time = 6 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z−axis r−axis time = 9 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z−axis r−axis time = 12

Figure 4.18: The drop is on the slip boundary. Clean drop (solid line) and contaminated drop (‘− · − · −·’) are in this figure.

(33)

In figure 4.18, the clean drop is very close to the contaminated drop since the shear rate is smaller than other previous cases. But the contact angle of the contaminated interface is smaller than the one of clean interface since the surfactant reduces the contact angle more quickly than the clean one. We can observe the phenomenon more clearly in figure 4.19. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z−axis r−axis time = 12

Figure 4.19: It is an enlargement in figure 4.18 at T = 12.

In figure 4.20, the surfactant concentration at moving contact line is higher than the other parts near moving contact line at first. The reason is that the moving contact line slips along the direction of positive r−axis, and the fluid flows along the same direction. Therefore, the surfactant concentrate in the vicinity of moving contact line. After T = 9, the surfactant diffuses from moving contact line to other parts since the contact angle is going to achieve the equilibrium contact angle.

In figure 4.21, the relative volume error is about 1.5% at T = 12. The total mass of

surfactant is conserved by our method. For the clean drop, the contact angle is about 35◦

at T = 12, but the contact angle is about 0for the contaminated drop since the surfactant

reduces the surface tension such that the contact angle at moving contact line achieves the equilibrium contact angle more quickly.

In fact, if we choose the larger shear rate, then the effect of the fluid velocity field is greater than the effect of the slip boundary. Hence, it is difficult to see that the contact angle at moving contact line achieves the equilibrium contact angle.

(34)

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.2 0.22 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 arc−length surfactant concentration time = 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.2 0.22 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 arc−length surfactant concentration time = 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.2 0.22 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 arc−length surfactant concentration time = 3 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.2 0.22 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 arc−length surfactant concentration time = 6 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.2 0.22 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 arc−length surfactant concentration time = 9 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.2 0.22 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 arc−length surfactant concentration time = 12

Figure 4.20: The time evolution of surfactant concentration depended on figure 4.18.

0 2 4 6 8 10 12 −1.5

−1 −0.5

0 Relative Error of Volume : %

Time 0 2 4 6 8 10 12 0.1767 0.1767 0.1767 0.1767 0.1767 0.1767 0.1767 0.1767 0.1767 0.1767 0.1767 Surfactant Mass Time 0 2 4 6 8 10 12 −0.5 0 0.5 1

Cosine Value of the Dynamic Contact Angle

Time 0 2 4 6 8 10 12 0 20 40 60 80 100 120

Dynamic Contact Angle

Time

Figure 4.21: This figure presents the relative volume error, the total mass of surfactant,

the variation of cos θd, and the variation of θd with time.(solid line : clean drop, ‘− · − · −·’

+7

Animal or vegetable fats and oils and their fractiors, boiled, oxidised, dehydrated, sulphurised, blown, polymerised by heat in vacuum or in inert gas or otherwise chemically

Milk and cream, in powder, granule or other solid form, of a fat content, by weight, exceeding 1.5%, not containing added sugar or other sweetening matter.

I want to start by saying there are two kinds of Olympic Games: Summer Olympic Games and Winter Olympic Games. But the Winter Olympic Games are much smaller than the

Estimated resident population by age and sex in statistical local areas, New South Wales, June 1990 (No. Canberra, Australian Capital

First Taiwan Geometry Symposium, NCTS South () The Isoperimetric Problem in the Heisenberg group Hn November 20, 2010 13 / 44.. The Euclidean Isoperimetric Problem... The proof

[r]

Salas, Hille, Etgen Calculus: One and Several Variables Copyright 2007 © John Wiley &amp; Sons, Inc.. All

6 《中論·觀因緣品》，《佛藏要籍選刊》第 9 冊，上海古籍出版社 1994 年版，第 1

The first row shows the eyespot with white inner ring, black middle ring, and yellow outer ring in Bicyclus anynana.. The second row provides the eyespot with black inner ring

• ‘ content teachers need to support support the learning of those parts of language knowledge that students are missing and that may be preventing them mastering the

• helps teachers collect learning evidence to provide timely feedback &amp; refine teaching strategies.. AaL • engages students in reflecting on &amp; monitoring their progress

Robinson Crusoe is an Englishman from the 1) t_______ of York in the seventeenth century, the youngest son of a merchant of German origin. This trip is financially successful,

fostering independent application of reading strategies Strategy 7: Provide opportunities for students to track, reflect on, and share their learning progress (destination). •

Strategy 3: Offer descriptive feedback during the learning process (enabling strategy). Where the

How does drama help to develop English language skills.. In Forms 2-6, students develop their self-expression by participating in a wide range of activities

volume suppressed mass: (TeV) 2 /M P ∼ 10 −4 eV → mm range can be experimentally tested for any number of extra dimensions - Light U(1) gauge bosons: no derivative couplings. =&gt;

• Formation of massive primordial stars as origin of objects in the early universe. • Supernova explosions might be visible to the most

(Another example of close harmony is the four-bar unaccompanied vocal introduction to “Paperback Writer”, a somewhat later Beatles song.) Overall, Lennon’s and McCartney’s

There are existing learning resources that cater for different learning abilities, styles and interests. Teachers can easily create differentiated learning resources/tasks for CLD and

Microphone and 600 ohm line conduits shall be mechanically and electrically connected to receptacle boxes and electrically grounded to the audio system ground point.. Lines in

• Thresholded image gradients are sampled over 16x16 array of locations in scale space. • Create array of

Digital PCR works by partitioning a sample into many individual real-time PCR reactions, some portion of these reactions contain the target molecules(positive) while others do