At the beginning of each time step, e.g., step n, the variables Xnk = X(αk, n∆t), γk+1/2n = γ(αk+1/2, n∆t), un = u(x, n∆t), and pn−1/2 = p(x, (n − 1/2)∆t) are all given. The details of the numerical time integration are as follows.
1. Compute the surface tension and unit tangent on the interface.
For a bubble:
σk+1/2n = σc¡
1 + ln¡
1 − ηγk+1/2n ¢¢
(4.46) τnk+1/2 = DαXnk
|DαXnk| (4.47)
both of which hold for αk+1/2 = (k +1/2)∆α. Then the interfacial force at the fluid/fluid markers Xk, k = 1, . . . , M are given by
Fnk = Dα(σnk+1/2τnk+1/2). (4.48) Note that Fn0 = FnM due to the periodicity of F .
For a drop:
σkn = σc(1 + ln (1 − ηγkn)) (4.49) τnk = DsXnk
|DsXnk| (4.50)
both of which hold for αk = k∆α, and γkn is simply approximated by the average of γk−1/2n and γk+1/2n . Then we compute the interfacial force at the fluid/fluid markers Xk, k = 1, . . . , M − 1 by
Fnk = Ds(σknτnk) . (4.51) The unbalanced Young force at the contact lines (k = 0 and k = M) can be computed by
Fn0 = (σs2− σs1− σn0 cos θ0n) e1 (4.52) FnM = − (σs2− σs1− σMn cos θMn ) e1 (4.53) where e1 = (1, 0) and cos θkn= −τnk· e1.
2. Distribute the force from the markers to the fluid by fn(x) =
M −1X
k=1
Fnkδh(x − Xnk) ∆α + C X
k=0,M
Fnkδh(x − Xnk) , (4.54)
where the smooth version of Dirac delta function in [36] is used. Here, C = 0.5 is used for a bubble enclosed in a fluid while C = 1 is used in the moving contact line problem.
3. Solve the Navier-Stokes equations. This can be done by the following second-order accurate projection method [2], where the nonlinear term is approximated by the Adams-Bashforth scheme and the viscous term is approximated by the Crank-Nicholson scheme. For convenience, we use a simple notation ψ|b[n]a[n−1] to represent bψn− aψn−1.
∇h· (µ∇hu∗) − 3Re
2∆tu∗ = ReHn (4.55)
where
Hn = un−1− 4un
2∆t + ((u · ∇h) u) |2[n][n−1]+ ∇hpn
− 1 Re
¡∇h·¡
µ∇huT¢¢
|2[n][n−1]− 1
ReCaf |2[n][n−1] (4.56) BCs for a shear flow:
u∗ = ry, on x = a and x = b, u∗ = rc, on y = c,
u∗ = rd, on y = d, v∗ = 0, on ∂Ω
(4.57)
BCs for a moving contact line problem:
u∗ = β∂u∗
∂y , v∗ = 0, on y = c, u∗ = 0, otherwise
(4.58)
∇2hφn+1= ∇h· u∗
2∆t/3 , ∂φn+1
∂n = 0, on ∂Ω (4.59)
un+1 = u∗− 2∆t
3 ∇hφn+1, (4.60) pn+1 = pn+ φn+1− 1
Re∇h· u∗. (4.61) Here, ∇h is the standard centered difference operator on the staggered grid. One can see that the above Navier-Stokes solver involves solving two Helmholtz equations for velocity u∗ = (u∗, v∗) and one Poisson equation for pressure. These elliptic equations are solved by using the geometric multigrid solver.
4. Interpolate the new velocity on the fluid lattice points onto the marker where (UA)n+1k is computed by a mid-point integration rule of Eq. (4.29).
It is worth noting that at the positions of contact lines k = 0 and k = M, we only update the x coordinates by the velocity along the solid wall, and keep their y coordinates fixed as Y0n+1 = YMn+1 = c.
This is consistent with the unbalanced Young forces which are applied only along the solid wall (x direction) as shown in Eqs. (4.52)-(4.53). As mentioned before, there is no artificial velocity applied at the moving contact lines; that is, (UA)n+10 = (UA)n+1M = 0.
5. Update surfactant concentration distribution γkn+1 by the numerical algorithm (4.37). Basically, we can rearrange (4.37) to a simpler form
à 1
2). Note that the unknown γn+1 appears only on the left-hand side of the above equation, we simply have a linear system Aγ = b. Since the essential data Xn+1and (UA)n+1are obtained in previous steps, we can get Sαn+1 involving to the tri-diagonal matrix A and the right-hand side of the system. The resultant symmetric tri-diagonal linear system can be solved easily by the Thomas algorithm.
Chapter 5
Numerical results
In this chapter, we first verify the convergence of our numerical schemes for the concerning problems, a bubble in a shear flow and a drop adheres to a solid substrate. If the surface tension is a constant, the motion of the interface is affected only by the curvature force in the normal direction (Capillary effect only), a corresponding benchmark will be shown after convergence test. When the interface is contaminated by surfactants, the surface tension can be altered and regarded as a function of surfactant which is derived by a linear or nonlinear Langmuir equation, therefore the maragoni effect appears. In the consequent, the effect of the surfactant for a bubble in a shear flow will be presented. Also, several moving contact line phenomena which are derived by the unbalanced Young’s forces will attract one’s attention.
Finally, we consider a two-phase flow with different fluid properties (both of the viscosity and density are piecewise constant functions). In this run, not only the surface tension evolve the interface but an external force (gravity) will affect the evolution of entire fluid field.
5.1 Convergent test
5.1.1 For a bubble in a shear flow
First of all, we carry out the parameter setting in the convergence study of the present method. Here, we perform different computations with varying Cartesian meshes h = ∆x = ∆y = 0.04, 0.02, 0.01, 0.005. 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 the analytical solution is not available in these simulations, we choose the results obtained from the finest mesh as our reference solution
and compute the L2 error between the reference solution and the solution obtained from the coarser grid. Table 5.2 shows the mesh refinement analysis of the velocity u, v, and the surfactant concentration γ. One can see that
h ku − urefk2 rate kv − vrefk2 rate kγ − γrefk2 rate 0.04 4.9739E-03 - 4.1656E-03 - 1.4551E-02 -0.02 2.1476E-03 1.21 1.8169E-03 1.20 6.3542E-03 1.20 0.01 6.9859E-04 1.62 6.2180E-04 1.55 2.2329E-03 1.51 Table 5.1: The mesh refinement analysis of the velocity u, v, and the surfac-tant concentration γ.
the error decreases substantially when the mesh is refined, and the rate of convergence is about 1.5. Notice that, the fluid variables are defined at the staggered grid and the surfactant concentration is defined at ”half-integer”
indices of the grid, so when we refine the mesh, the numerical solutions will not coincide with the same grid locations. In these runs, we simply use a linear interpolation to compute the solutions at the desired locations. We attribute this is part of the reason why the rate of convergence behaves less than second-order. However, we believe that the key factor reducing the rate of convergence comes from using the discrete delta function, a function with compact support and only C1 smoothness.
5.1.2 For the moving contact line problem
Since the motion of the interface in the moving contact line problem involves not only fluid-fluid surface tension but surface tension between fluid and solid near the contact line points. There are additional treatments in the numer-ical processes, that is why we have two verifications for the convergence of the numerical schemes.
We consider a computational domain Ω = [−1, 1] × [0, 1] where a half of circular drop with radius 0.5 is initially attached on the bottom of the domain, and both left and right contact angles are π/2 initially. The initial fluid velocity are all set to be zero. The surface tension is affected by the sur-factants following the equation of state given by Eq. (2.38), here, σc= 1 and η = 0.3. Other parameters are initial surfactant concentration γ(α, 0) = 1, the Reynolds number Re = 10, the capillary number Ca = 0.1, the surface Peclet number P es = 20, and the slip length β = h/4, where h is the mesh width. Here, we perform four different computations with varying Carte-sian mesh h = ∆x = ∆y = 1/32, 1/64, 1/128, 1/256. The Lagrangian mesh
h ku − urefk2 rate kv − vrefk2 rate kγ − γrefk2 rate 1/16 5.8079e-03 - 3.3148e-03 - 3.1818e-02 -1/32 2.9639e-03 0.97 1.9179e-03 0.79 1.7977e-02 0.82 1/64 1.4773e-03 1.00 1.0805e-03 0.82 9.8698e-03 0.86 1/128 5.9179e-04 1.32 4.6628e-04 1.21 4.2087e-03 1.23 Table 5.2: The mesh refinement analysis of the velocity u, v, and the surfac-tant concentration γ.
h kX − Xrefk∞ rate | cos(θ) − cos(θref)| rate |A − Aref|/Aref rate
1/16 4.7646e-02 - 1.0285e-01 - 1.1476e-01
-1/32 2.4225e-02 0.98 6.3948e-02 0.69 5.6724e-02 1.02 1/64 1.1426e-02 1.08 3.6915e-02 0.79 2.7473e-02 1.05 1/128 4.1700e-03 1.45 2.1708e-02 0.77 1.3145e-02 1.05 Table 5.3: The mesh refinement analysis of interface positions, the contact angles, and the area of drop.
is chosen as ∆α ≈ h and the time mesh size is ∆t = h/10. The solutions are computed up to time T = 6.25. Again, the analytical solution is not available in these simulations, we choose the results obtained from the finest mesh as our reference solution and compute the L2 error between the ref-erence solution and the solution obtained from the coarser grid. Table 5.2 shows the mesh refinement analysis of the velocity u, v, and the surfactant concentration γ. One can see that the error decreases substantially when the mesh is refined, and the rate of convergence is about first-order. Notice that, the fluid variables are defined at the staggered grid and the surfactant con-centration is defined at ”half-integer” indices of grid, so when we refine the mesh, the numerical solutions will not coincide with the same grid locations.
In these runs, we simply use a linear interpolation to compute the solutions at the desired locations. Table 5.3 shows the L∞ errors of the interface po-sitions, the cosine value of the contact angles, and the area loss of the drop for different meshes. One can see the rate of convergence is about first-order which is consistent with the results shown in Table 5.2.
During the drop deformation, the Lagrangian markers will gradually sweep into contact lines and cause clustered distribution near those points.
Fig. 5.1 shows the comparison of the local stretching factor |Xα| = Sα at T = 6.25 with (solid line) or without (dashed line) the implementation of
0 1 0.05
0.1 0.15 0.2 0.25 0.3 0.35 0.4
Sα comparison at T = 1.5625
Figure 5.1: Comparison of stretching factor |Xα| with (UA = 0, dashed line) and (UA6= 0, solid line), where h = 1/128.
equi-distributed technique of Lagrangian markers. Here, we only show the case with surfactant η = 0.3. One can see this technique does preserve the markers more uniformly.
In addition, it is well-known that the fluid leakage often appears in the simulation of immersed boundary method. In [37], Peskin and Printz pro-posed an improved volume (area in 2D) conservation scheme for the immersed boundary method by constructing a discrete divergence operator based on the interpolation scheme. Here, however, the area loss is not that significant, thus no modification is applied. One should mention that with the implemen-tation of present dynamical control of Lagrangian markers, one can reduce the area loss significantly than without implementing it.