二、 The governing equations
2.3 Controlling the interfacial markers uniformly … 6
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
If we allow surfactant diffusion along the interface, we obtain the surfactant transport-diffusion equation as
r 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,
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 2π
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 Xnk = 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 are as follows.
1. Compute the surface tension and unit tangent vector on the interface as σnk = σc 1 + ln 1 − βΓnk
, (3.2)
τnk = DsXnk
DsXnk , (3.3)
both of which hold for sk−1/2 = (k − 1/2)∆s. Then we define the interfacial force as Fnk = 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˜∆˜un+1+ 1
ReCafn, (3.7)
˜u = ub, on ∂Ω, (3.8)
˜∇h· un+1= 0, (3.9)
˜∇2hφn+1= 3 ˜∇h· ˜u 2∆t , ∂φ
∂n = 0, on ∂Ω (3.10)
˜∇hφn+1= 3
˜un+1− un+1
2∆t , (3.11)
pn+1= pn+ φn+1− 1
Re˜∇h· ˜un+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∂r∂r,∂z∂)· and the Laplace operator
∆˜h˜un+1= ∂2un+1r
∂r2 + 1 r
∂un+1r
∂r + ∂2un+1r
∂z2 − un+1r
r2 , ∂2un+1z
∂r2 + 1 r
∂un+1z
∂r + ∂2un+1z
∂z2
!
˜∇2hφ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+1k = 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+1k = k∆s 2π
XM k0=0
DsUn+1k · τn+1k ∆s − Xk k0=0
DsUn+1k · τn+1k ∆s. (3.15) From the above equations, we have the modified marker positions,
Xn+1k = Xnk + ∆tUn+1k + ∆tUAn+1k τn+1k . (3.16) 5. Update surfactant concentration distribution Γn+1k . 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)
Now we discretize the above equation by the Crank-Nicholson scheme in a sym-metric way as
Γn+1k − Γnk new interfacial marker position Xn+1k obtained in the previous step and the boundary condition (Γn+1k − Γ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+1k |DsXn+1k |rn+1k−1/2∆s =X
k
Γnk|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+1k − Γnk
4 Numerical results
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
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, Ca = 0.25, and Pes = 12.5.
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.
0 0.2 0.4 0.6 0.8 1
z−axis , z−right−point = 0.4
r−axis
z−axis , z−right−point = 0.43685
r−axis
z−axis , z−right−point = 0.45546
r−axis
z−axis , z−right−point = 0.46548
r−axis
z−axis , z−right−point = 0.47121
r−axis
z−axis , z−right−point = 0.47455
r−axis
z−axis , z−right−point = 0.47645
r−axis
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 number is larger than the critical value. The Capillary number Ca ≈ µGr/σc is relative
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.
z−axis , z−right−point = 0.4
r−axis
z−axis , z−right−point = 0.4949
r−axis
z−axis , z−right−point = 0.54999
r−axis
z−axis , z−right−point = 0.59339
r−axis
z−axis , z−right−point = 0.63507
r−axis
z−axis , z−right−point = 0.68175
r−axis
z−axis , z−right−point = 0.74252
r−axis
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
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.
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.
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.
0 0.2 0.4 0.6 0.8 1
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.
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.
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.
0 0.2 0.4 0.6 0.8 1
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
Figure 4.10: The difference of the drop with different Peclet numbers
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 θs+σsl = σ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
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.
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.