In this section we will show our numerical results of simulating two-phase flows. Here we thank [22] who provided fluid code of solving Navier-Stokes equation with different densities ρ1 and ρ2 and viscosities µ1 and µ2.
Example.1
Our first example is to test single bubble simulation. Here we give our setting.
Inside interface:ρ1 = 1,µ1 = 0.02 and outside interface:ρ2 = 40, µ2 = 0.005 and let σ = 0.01. We let initial interface to be a circle which center at (0, −0.2) with radius r = 0.15. We take ∆t = 14∆x.
Figure 19: Simulation of example.1
Note that the inner area variation is less than 0.1% and if area variation is positive then it means area increases. On the other hand, if area variation is negative then area decreases. We contour φ = 0 in this and later examples.
Example.2
In this example, we simulate the motion of two bubbles in two phases. We use two circles to represent interface of two bubbles. If the initial distance between two bubbles is far enough, then the two bubbles will not merge each other. Note that [8] had details about how far that two bubbles will not merge.
Figure 20: Simulation of example.2
We use level function to represent the interface of two circles and we first choose the case which do not merge
φ1(x, y) = (p
x2+ (y − 0.1)2− 0.15) × (p
x2 + (y + 0.3)2− 0.15) Note that φ1is not signed distance function and we need apply reinitialization process at the beginning.
Here we give our setting. Inside interface:ρ1 = 1, µ1 = 0.02 and outside interface:ρ2 = 40, µ2 = 0.00025, σ = 0.005 and ∆t = 14∆x. Note that since the deformation of bubble shape is much heavier and so the area variation increase more. But the area variation still less than 1% which is also a good result.
Example.3
In this example, our setting is just as example.1 except the level function and we let the distance between two bubbles to be smaller. So the level function is
φ2(x, y) = (p
x2+ y2− 0.15) × (p
x2 + (y + 0.3)2− 0.14)
We find that two bubbles merge with each other and then break up into two parts. This example tells us a very important information which is that this coupled VOF and level set method also can be applied to the changing of interface topology like merging and breaking and the inner volume variation is still less than 1%.
Example.4
In this example we apply steady shear flow to simulation one bubble in two-phase flow. The circle is centered at (0, 0) with radius 0.1 in the compu-tational domain [−0.5, 0.5] × [−0.5, 0.5] and we set the boundary condition of velocity field as ub = (2y, 0). We set the gravity g = 0, σ = 0.1, ρ1 = ρ2 = 1 and µ1 = 0.1, µ2 = 1.Note that the area variation is less than 10−3% since we let the time step ∆t = 401∆x.
Figure 21: Simulation of example.3
Figure 22: Simulation of example.4
5 Insoluble Surfactant on Interface
5.1 Introduction
Since we have introduced the VOF method and Local Level Set Method in above chapter, we now couple these two methods to solve insoluble surfactant on interface. Note that we have introduced how to solve heat equation on circle which is static(i.e. interface do not move) and in this chapter interface will be moved by velocity field. So the surfactant problem is more complicate than heat equation on circle since the interface is no more static.
We refer to [21] which uses level set method to handle the interface proper-ties and solves surfactant equation by local level set method. We also couple VOF into level set method to conserve inner area of bubble. Since we use local level set method as [21], so we can expect that the mass of surfactant will not be conserved in our numerical scheme and this is the drawback of local level set method. However the local level set is still powerful since it just defined on Cartesian grid and the computational domain is just a band of interface.
5.2 Governing Equations with Insoluble Surfactant
We want to solve Navier-Stokes equation with two-phase flow and surfactant on the interface, so the Navier-Stokes equation becomes
ρ(∂u
∂t + u · ∇u) = −∇p + ρg − ∇ · µ[∇u + (∇u)T] + σκn − ∇sσ (61)
∇ · u = 0
where (σκn − ∇sσ) is called Laplace-Young equation and since the effect of surface tension is just on the interface, we approximate σκn by σκ∇H where H is a heaviside function which is defined as above chapter. Note that σκn is the capillary force and ∇sσ is Marangoni force where ∇s is surface gradient which is ∇s = (I − n ⊗ n)∇ that we have mentioned above.
Since the surfactant will depress the surface tension coefficient σ, so the surfactant and σ have a relation.We refer [21] to write down the relation :
σ(Γ) = σ0− RT Γ (62)
where Γ is surfactant concentration, σ0 is surface tension coefficient for clean interface which means without surfactant, R is the ideal gas constant and T
is absolute temperature. Then we refer [22] to non-dimensionalize Eq.(62) as the form :
σ(Γ) = σ0(1 − βΓ) (63)
where β satisfies 0 ≤ β ≤ 1.We can know that from Eq.(63), if surfactant concentration Γ is lager, then the surface tension coefficient σ is smaller.
Since the surfactant will diffuse and convect along the interface, so we can have the diffusion-convection equation of surfactant as follow:
Γt+ u · ∇Γ − n · (∇un)Γ = 1
P e∆sΓ (64)
where P e is Peclet number.
Next we describe how to calculate the mass of surfactant. Since the surfactant is just defined on interface which we use φ to describe. The surface integral of Γ on the interface ∂Ω is
Z
∂Ω
Γ(x)ds = Z
Γ(x)δ(φ)|∇φ|dx (65)
where δ is a 1Dδ-function which we refer [18] as
δ(x) = So the total mass of surfactant M is
M =
When solving Eq.(64) on moving interface relative to two-phase flow, the process is complicate since there are many equations needed to solve and we couple VOF and level set to handle the properties of interface, so we describe the outline of whole method for one time step as below :
Step1 Solve the Navier-Stokes equation to get the velocity u and note that our velocity u is defined on the cell face.(e.g.solve Eq.(61))
Step2 Solve the advection equation of VOF by finite volume method and reconstruct the interface by PLIC method.
Step3 Move the interface of level function by the velocity u.
Step4 Reinitialize φ to be signed distance function within a band by the VOF method.
Step5 Extend the surfactant Γ off interface by normal extension method.
Step6 Evaluate surfactant equation within the band of interface,that is to solve Eq.(64) using local level set method.
Step7 Evolve Eq.(63) and solve ∇sσ. Note that we also define σ in a band of interface and solve ∇sσ by Cartesian coordinate. Finally we return to Step1.
The process of whole method is as above and there are seven steps to handle, so the priority of solving which equation must be decided.
5.4 Algorithm of The Method
We have given the details from Step1 to Step6 in above chapters and we just write the algorithm of Step7 since the surfactant equation is a little different as comparing with heat equation.
The surfactant equation can be rewritten as
Γt+ u · ∇Γ − n · (∇un)Γ = 1
P e(∆Γ − ∂2Γ
∂n2 − κ∂Γ
∂n) (66)
and again we use semi-implicit Crank-Nicholson method to solve Eq.(64) which we have mentioned above. Note that Crank-Nicholson is just used in diffusion part which is ∆sΓ and we treat the convection part with explicit discretization, then the discretization form is :
Γn+1− Γn
We refer [21] to write the form of Eq.(66) and [21] suggested that in the discretization for spatial of Eq.(66), central difference is used in all term except this term u·∇Γ which is discretized by upwind-WENO-3 scheme which we have introduced. Note that we encounter the data structure problem again as in above chapter.
Note that σx and σy are discretized by central difference.
5.5 Numerical Results
The effect of surfactant is to decrease the surface tension coefficient and so the interface will deform much heavier than clean interface.
In this section we refer the example of [22] and we apply the steady shear flow with boundary condition ub = (2y, 0) in the computational domain [−0.5, 0.5] × [−0.5, 0.5].The initial bubble is centered at (0, 0) with radius r = 0.15.
Example.1
We give our setting as follows.Initial surface tension coefficient σ0 = 0.5, gravity g = 0, ρ1 = ρ2 = 1, µ1 = 0.1, µ2 = 1 and Peclet number P e = 1.
Note that our time step is as the example of above chapter, ∆t = 401 ∆x and we consider three cases that are β = 0(clean), β = 0.25 and β = 0.5(See figure23).
As our expectation, if β is larger then the interface will deforms heavier since we know that
σ(Γ) = σ0(1 − βΓ)
which means that concentration affects the surface tension coefficient.
Note that our inner area is conserved within 0.01% but the surfactant concentration which computed by Eq.(64) will not be conserved as we predict before.(See figure24)
Figure 23: Time evolution of a bubble in shear flow with β = 0(black), β = 0.25(blue) and β = 0.5(red)
6 Conclusion and Future Work
In this article we introduce some methods to simulate two-phase flow with or without surfactant. We also introduce level set method to solve heat equation on a circle and we can extend the level set method to 3D directly;
that is to solve heat equation on sphere φ(x, y, z) = px2+ y2+ z2 − 1. We can also couple level set method and IIM(immerse interface method) to solve heat equation on the 2D domain with jump condition on the interface. The application of level set is very widely and it is really a valuable method.
The VOF method used in this article is to conserve the inner area and we get nice result in our numerical simulation. By coupling VOF and level set method, we are not only able to reconstruct the interface accurately but also to conserve the volume of inner area. Moreover, the method is also adapted to changing of interface topology such as two bubbles merging or one bubble braking into two parts.
Figure 24: (a) Total surfactant in each time step for β = 0.25. (b) Total area of drop in each time step for β = 0.25. (c) Total surfactant in each time step for β = 0.5. (d) Total area of drop in each time step for β = 0.5.
In our simulation of two-phase flow with surfactant on interfaces, the inner area is still conserved as we expect but the mass of surfactant is not conserved which does not make sense in our real life. So our future work is to develop another method to let the mass of surfactant to be conserved, that is we use other method to handle Eq.(64) and get feasible result in mass of surfactant.
References
[1] S. Osher, R. Fedkiw (2003). Level Set Methods and Dynamic Implicit Surfaces, Springer
[2] J. A. Sethian (1996). Level Set Methods and Fast Marching Methods, Cambrige
[3] C. W. Hirt, B. D. Nichols (1981). Volume of Fluid Method for the Dy-namics of Free Boundaries, J. Comput. Phys 39, 201-225
[4] C.-W. Shu (1989). Efficient Implementation of Essentially Non-oscillatory Shoch-Capturing Schemes II, J. Comput. Phys 83, 32-78 [5] S. Osher, C.-W. Shu (1991). High-Order Essentially Nonoscillatory
Schemes for Hamilton-Jacobi Equations, SIAM J. on Num. Analysis, Vol.28,No.4 907-922
[6] M. Sussman, P. Smereka, S. Osher (1994). A Level Set Approach for Computing Solutions to Incompressible Two-Phase Flow, J. Comput.
Phys 114, 146-159
[7] H. K. Zhao, T. Chan, B. Merriman, S. Osher (1996). A Variation Level Set Approach to Multiphase Motion, J. Comput. Phys 127, 179-195 [8] L. Chen, S. V. Garimella, J. A. Reizes, E. Leonardi (1997). Motion of
Interacting Gas Bubbles in A Viscous Liquid Including Wall Effects and Evaporation, Numer. Heat Trans,Part A, 31, 629-654
[9] S. Chen, B. Merriman, S. Osher, P. Smereka (1997). A Simple Level Set Method for Solving Stefan Problems, J. Comput. Phys 135, 8-29
[10] R. Scardovelli,S. Zaleski (1999). Direct Numerical Simulation of Free-Surface and Interfacial Flow, Annu. Rev. Fluid Mech. 31, 567-603 [11] D. Gueyffier, J. Li, A. Nadim, R. Scardovelli, S. Zaleski (1998).
Volume-of-Fluid Interface Tracking with Smoothed Surface Stress Methods for Three-Dimensional Flows, J. Comput. Phys 152, 423-456
[12] D. Peng, B. Merriman, S. Osher, H. K. Zhao, M. Kang (1999). A PDE-Based Fast Local Level Set Method, J. Comput. Phys 155, 410-438 [13] G.-S. Jiang, D. Peng (2000). Weighted ENO Schemes for
Hamilton-Jacobi Equations, SIAM J. Sci. Comput. Vol.21, No.6 2126-2143
[14] M. Sussman, E. G. Puckett (2000) A Coupled Level Set and Volume-of-Fluid Method for Computing 3D and Axisymmetric Incompressible Two-Phase Flows, J. Comput. Phys 162, 301-337
[15] M. Beralmio, L.-T. Cheng, S. Osher (2001) Variational Problems and Partial Differential Equations on Implicit Surfaces, J. Comput. Phys 174, 759-780
[16] G. Son (2001) A Numerical Method for Bubble Motion with Phase Change, Numer. Heat Trans,Part B, 39, 509-523
[17] G. Son, N. Hur (2002) A Coupled Level Set and Volume-of-Fluid Method for the Buoyancy-Driven Motion of Fluid Particles, Numer.
Heat Trans,Part B, 42, 523-542
[18] J.-J. Xu, H.-K. Zhao (2002) An Eulerian Formulation for Solving Par-tial DifferenPar-tial Equations Along A Moving Interface, J. of Sci. Comp.
Vol.19, 573-594
[19] J. B. Greer (2006) An Improvement of A Recent Eulerian Method for Solving PDEs on General Geometries, J. of Sci. Comp. Vol.29, 321-352 [20] H. Tseng (2006) Numerical methods and applications for immersed
in-terface problems, Master Thesis.
[21] J.-J. Xu, Z. Li, J. Lowengrub, H.-K. Zhao (2006) A level-set method for interfacial flows with surfactant, J. Comput. Phys 212,590-616
[22] M.-C. Lai, Y.-H. Tseng, H. Huang , An immersed boundary method for interfacial flows with insoluble surfactant, J.Comput.Phys., in press 2008