二、 Level Set Method
2.5 Navier-Stoke Equations
2.5.2 Marker and Cell Method
Refer to [4], [9] and [12]. We use marker and cell (MAC) method and let the computation domain Ω be squared with staggered grid in Cartesian coordinates. This specially defined grid decomposes the computational domain into cells with velocities and forces defined on the cell faces and scalars defined at cell centers such as pressure, density and viscosity. We take the mesh size h = ∆x = ∆y and let 𝑁𝑥, 𝑁𝑦 be the
number of mesh in the x and y directions respectively. In general, we take the same number of grid points in each direction, and denote that 𝑁 = 𝑁𝑥 = 𝑁𝑦. MAC grid is as the following diagram comes from [9],
Figure 5 : The computational domain Ω using staggered grid with mesh h.
First Part : Density and Viscosity Are Constant 2.5.3 Introduce Navier-Stokes Equations
The Navier-Stokes equations are
∂w
∂t + w ∙ ∇ w + ∇𝑝 = 1
𝑅𝑒∆w + 𝑓
∇ ∙ w = 0 (7) where w = 𝑣(𝑡,𝑥,𝑦)u(t,x, y) = the velocity in the 𝑥 − direction
the velocity in the 𝑦 − direction
and p : p(t,x,y) denotes the pressure. f(t,x,y) is the external force defined on the computational domain, which contains the forces resulted from the following, a. Surface tension. The form is τκ(ϕ)∇ϕδ(ϕ), where τ is the surface
tension coefficient, κ is curvature,
and κ ϕ = ϕy2ϕxx−2ϕ(ϕ xϕyϕxy+ϕx2ϕyy
x2+ϕy2)3/2 ,
δ is a one-dimension Dirac delta function and ϕ is our signed distance function.
b. Gravity.
c. Frictional force of moving contact line.
2.5.4 Governing equations
Therefore we need to solve the coupling system :𝑢n+1 − 𝑢n
Similar arguments for v can be done.
Discretizing (8), we obtain
𝑢i+1 diffusion term. Therefore, we choose
∆t = τ ∙ min {𝑅𝑒
[𝑢t]n+1 =3𝑢n +1−4𝑢n+𝑢n −1
2∆t + 𝑂(∆t2), and 𝑢𝑢x+ 𝑣𝑢y n+1 = 𝑢𝑢x+ 𝑣𝑢y |n−1n
= 2[𝑢𝑢x + 𝑣𝑢y]n − [𝑢𝑢x+ 𝑣𝑢y]n−1 + 𝑂(∆t2).
Therefore we need to solve the coupling system :
3𝑢n +1−4𝑢n+𝑢n −1
We use second-order projection method to solve Navier-Stokes equations.
Step1. Using Helmholtz-type solver to solve
3w ∗−4w n+w n −1
2∆t + w ∙ ∇ w |n−1n + ∇𝑝n =Re1 ∆w ∗ (9) Where w ∗ is the intermediate velocity field between w n and w n+1. Here w ∗ is solved implicitly, the constraint of ∆t will be diminished.
Step2. From Hodge decomposition, since w ∗ is a vector field, there are a potential function φn+1 and a divergence-free vector field w n+1 such that
w n +1−w ∗
2∆t/3 = −∇φn+1 (10) And if we take the divergence operator into (10), we obtain
∆φn+1 =2∆t3 ∇ ∙w ∗ (11) new schemes and the coupling system term by term, we obtain
∇𝑝n+1 = ∇𝑝n + ∇φn+1−Re1 ∆w ∗ (13)
Second Part : Density and Viscosity Are Not Constant
We solve Navier-Stokes equations by two steps as follows.
Step1. Prediction:
First, we solve it by Crank-Nicholson method in time
∇ ∙ 𝜇n+1∇𝐮∗ − 2ρn+12𝐮∆t∗ = 2∇𝑝n−12− 2ρn+12𝐮∆tn − ∇ ∙ (𝜇∇𝐮)n + (3𝐅n − 𝐅n−1) (15)
and (15) becomes Poisson type equation and it can be solved by linear system for 𝐮∗. Step2. Projection:
ρn+12𝐮n+1− 𝐮∗
∆t = −∇𝜓n+1 and then take divergence we have
∇ ∙ ∇𝜓n+1 ρn+12
= ∇ ∙ 𝐮∗
∆t After obtaining ∇𝜓n+1 we have
𝐮n+1 = 𝐮∗− ∆t ρn+12
∇𝜓n+1
∇𝑝n+12 = ∇𝑝n−12+ ∇𝜓 + ∇ ∙ (𝜇n+1(𝐮n+1− 𝐮∗))
2.6 Level set Equation
In order to define the evolution of our implicit function ϕ, we use the simple level set equation ϕt + u ∙▽ϕ = 0, where ϕ is a signed distance function, and u is the velocity. By solving the Navier-Stokes equations, we can obtain the velocity u of each point on computational domain at this time step. Therefore, according to the equation ϕt+ u ∙▽ϕ = 0, we can get ϕ at next time step. Here, we use HJ WENO to discretize the spatial terms to fifth-order accuracy and TVD Runge-Kutta to get third-order accuracy in time.
2.7 Summary
Refer to figure 6, our procedure is as follows.
1. Construct the level set function on computational domain where ϕ 𝑥 = 0 stands for the interface.
2. In order to adjust the level set function to a signed distance function, we must re-initialize the level set function. It’s important that re-initialization doesn’t change the shape of interface.
3. Define Heaviside function and use it to define the corresponding regularized density function and the regularized viscosity function.
4. Solve the Navier-Stokes equations to obtain the velocity of each point on computational domain at this time step.
5. Using the velocity obtained from Navier-Stokes equations and combining it with the level set equation, we can get the level set function at next time step.
6. Remark that the level set function obtained from process 5 is not a signed distance function, therefore we should return to process 2.
Figure 6
3 Moving Contact Line Problems
3.1 Introduction
The moving contact line problem is very practical and interesting, because it is in our daily life. For example, the capillarity, lubricated pipeline transport, fibers, and the bead on a leaf. We know that a bead on the floor and a drop of mercury on the floor whose shapes are different clearly, this phenomenon is mainly caused by surface tension.
Figure 7
Refer to [2]. The surface tension occurs on interface, therefore as the above diagram, let ςsg, ςsl and ς be the surface tension coefficients between solid-gas, solid-liquid, and gas-liquid phases, respectively. Because the different surface tensions, there are three forces at points A as Figure 6. These forces would attain to balance, at this time, ςsg, ςsl, ς and the contact angle satisfy the well-known Young’s equation
ςsg = ςsl + ς ∙ cosθ0 (16) here θ0 is the static equilibrium contact angle.
When the shape of interface changes, the contact angle θ will differ from θ0 until an equilibrium is established. Let θ𝑑 represent the dynamic equilibrium contact angle at the contact line, the friction force Fs at the contact line satisfies the following equation
Fs = ςsg − ςsl− ς ∙ cosθ𝑑 (17) By (16) and (17), we can obtain that Fs = ς ∙ (cosθ0− cosθ𝑑).
It’s clearly there are two challenges needed to be overcome as following if we use level set methods to solve moving contact line problems.
1. How to calculate the dynamic equilibrium contact angle using the level set method?
2. How to distribute the friction force at the contact line using the level set method?
3.2 Calculate Contact Angle
In order to obtain the friction force at the contact line, from equation (17), we need to calculate the dynamic equilibrium contact angle at the contact line using level set method.
First, we know that interface which is near the boundary approximates tangent line of the interface. It’s reasonable if the mesh size is small enough.
Here, we assume that the level set function is a signed distance function. As the following diagram,
Figure 8
where θ is the contact angle, and the signed distance function is defined at cell center.
What we know are y, z and mesh size ∆x and unknown is h, by the property of can’t ascertain that cosθ is positive or negative. In other words, we should determine that contact angle is acute or obtuse.
Refer to figure 7, we define the part of signed distance function which represents liquid is negative, the other part which represents gas is positive. Remark that signed distance function is defined at cell centers. Let ϕ 𝑖, 𝑗 be our signed distance function, if ϕ 𝑖, 1 > 0 and ϕ 𝑖 + 1,1 < 0 for some 𝑖, it shows that interface is between the two points. For this 𝑖, if ϕ 𝑖, 2 > ϕ 𝑖, 1 , we can say that contact angle is acute, otherwise, the contact angle is obtuse. As the following diagrams,
Figure 9 : ϕ 𝑖, 2 > 0, ϕ 𝑖, 1 > 0, ϕ 𝑖 + 1,1 < 0 and ϕ 𝑖, 2 > ϕ 𝑖, 1 , θ is acute.
Figure 10 : ϕ 𝑖, 2 < 0, ϕ 𝑖, 1 > 0, ϕ 𝑖 + 1,1 < 0 and ϕ 𝑖, 2 < ϕ 𝑖, 1 , θ is obtuse.
Figure 11 : ϕ 𝑖, 2 > 0, ϕ 𝑖, 1 > 0, ϕ 𝑖 + 1,1 < 0 and ϕ 𝑖, 2 < ϕ 𝑖, 1 , θ is obtuse.
By the same way, we can calculate the dynamic equilibrium contact angle using level set method in the other case. Therefore, we have Fs .
3.3 Distribute the Friction Force at the Contact Line
We have calculated the friction force at the contact line, but the force terms are defined at cell faces. Therefore we need to distribute the friction force at the contact line using discrete delta function as (5)
δε ϕ = 2ε1 1 + cos πϕ
ε if |ϕ| < 𝜀 0 otherwise Here, taking ε = 2.5∆x, ∆x is the mesh size.
Now, we should calculate distance between the source of force and the point of MAC grid where the force terms are defined at. As the following diagram
Figure 12
According to our MAC grid, force is defined at cell faces like point A and B. Signed distance function is defined at cell center like point a and b, therefore y and z we have known. Now, we need to calculate Ao and Bo .
cb = 𝐳
𝐲 + 𝐳∙ ∆x dc = do ∙ cotθ = 1
2∆x ∙ cotθ dB = dc + cb − Bb = 1
2∆x ∙ cotθ + 𝐳
𝐲 + 𝐳∙ ∆x −1 2∆x dA = ∆x − dB = 3
2∆x −1
2∆x ∙ cotθ − 𝐳
𝐲 + 𝐳∙ ∆x Ao = dA 2+ do 2
Bo = dB 2+ do 2
By the same way, we can calculate distance between the source of force and the point of MAC grid where the force terms are defined at. To combine discrete delta function (5), it’s easy to distribute the friction force at the contact line.
Here, we let the incidence of the friction force at the contact line be a circle with radius 2.5∆x . Refer to the following diagram,
Figure 13
By section 3.2, we can obtain the friction force at the contact line Fs . By section 3.3, we can obtain the lengths of F , … , FsF1 respectively. sF9
We know that the incidence of discrete delta function (5) is a circle, but we let the force outside of the boundary be zero. In order to maintain the integral of delta function is 1, we let the discrete delta function inside of the boundary multiply by 2.
Taking the force F1 which is resulted from the friction force at the contact line as an example,
F1 = Fs × 2 × δε F sF1 where ε = 2.5∆x .
4 Numerical Results
4.1 A Rising Bubble
In this easy case, we try to simulate a rising bubble in the water under the influence of gravity. The fluid is at a standstill initially. Our computational domain Ω = 0,1 × [0,2], using 128× 256 grid points with ∆t = 0.25∆x. Viscosity for the fluid inside the bubble is equal to μ = 0.001873. Viscosity for the fluid outside the bubble is equal to μ = 0.00089. We take the density inside the bubble to be 0.8 and the density outside the bubble to be 1. The surface tension is 0.05. The initial position of the bubble is a circle. The circle is centered at (0.5, 0.65) with radius 0.15. In Figure 4.1, we plot the evolution of the bubble at time T = 0,100128,200128,300128,400128,500128,600128,700128 .
Figure 4.1 : The results for the above examples.
4.2 Merging Two Bubbles with the Same Density
Here, we compute the interaction of two bubbles of the same density under the influence of gravity. Our computational domain Ω = 0,1 × [0,2], using 128× 256 grid points with ∆t = 0.25∆x. The fluid is at a standstill initially. Viscosity for the fluid inside the two bubbles is equal to μ = 0.001873. Viscosity for the fluid outside the two bubbles is equal to μ = 0.00089. We take the density inside the two bubbles to be 0.8 and the density outside the bubbles to be 1. The surface tension is 0.05. The initial position of the two bubbles is two circles. The lower one is centered at (0.5, 0.35) with radius 0.1. The upper one is centered at (0.5, 0.65) with radius 0.15. In Figure 4.2, we plot the evolution of the two bubbles at time T = 0,250512,500512,750512,1000512,
1250
512 ,1500512 ,1750512 .
Figure 4.2 : The results for the above examples.
4.3 Rayleigh-Taylor Instability ( I )
Here, we compute the Rayleigh-Taylor instability under the influence of gravity. Our computational domain Ω = 0,1 × [0,2], using 128 × 256 grid points with ∆t = 0.025∆x. The fluid is at a standstill initially. Viscosity for the upper fluid is equal to μ = 0.05. Viscosity for the lower fluid is equal to μ = 0.01. We take the density of the upper fluid to be 10 and the density of the lower fluid to be 1. The surface tension is 0.05. Initially, the two kinds of fluids are separated by the interface. The interface is a part of circle. The circle is centered at (0.5, 2.3) with radius 1.2. In Figure 4.3, we plot the evolution at time T = 0,2569 ,25618 ,25627 ,25636 ,25645 ,25654 ,25663 .
Figure 4.3 : The results for the above examples.
4.4 Rayleigh-Taylor Instability ( II )
In this case, the two kinds of fluids are separated by the interface initially. The interface is a part of circle. The circle is centered at (0.5, -0.3) with radius 1.5. All other parameters are as case 4.3. In Figure 4.4, we plot the evolution at time T = 0,2569 ,25618 ,25627 ,25636 ,25645 ,25654 ,25663 .
Figure 4.4 : The results for the above examples.
4.5 Emit ( I )
Here, we try to simulate the air bubble emits from water surface under the influence of gravity. Our computational domain Ω = 0,1 × [0,2], using 128× 256 grid points with ∆t = 0.025∆x. The fluid is at a standstill initially. Viscosity for the fluid A is equal to μ = 0.001. Viscosity for the fluid B is equal to μ = 0.5. We take the density of the fluid A to be 1 and the density of the fluid B to be 10. The surface tension is 0.05. The bubble is a circle which is centered at (0.5, 0.6) with radius 0.1. In Figure 4.5, we plot the evolution of the bubble at time T = 0,120512,276512,288512,300512,312512,
324
512,360512 .
Figure 4.5 : The results for the above examples.
4.6 Emit ( II )
Most parameters are as 4.5. But we take the density of the fluid A to be 1 and the density of the fluid B to be 100 here. In Figure 4.6, we plot the evolution of the bubble at time T = 0,120512,240512,246512,252512,258512,264512,360512 .
Figure 4.6 : The results for the above examples.
4.7 Drip ( I )
Here, we try to simulate the water drop drips into the water under the influence of gravity. Our computational domain Ω = 0,1 × [0,2], using 128× 256 grid points with ∆t = 0.025∆x. The fluid is at a standstill initially. Viscosity for the fluid A is equal to μ = 0.089. Viscosity for the fluid B is equal to μ = 0.1873. We take the density of the fluid A to be 10 and the density of the fluid B to be 8. The surface tension is 0.1. The bubble is a circle which is centered at (0.5, 1.2) with radius 0.1. In Figure 4.7, we plot the evolution of the bubble at time T = 0,330512,340512,350512,360512,
370
512,380512,500512 .
Figure 4.7 : The results for the above examples.
4.8 Drip ( II )
Most parameters are as 4.7. But we take the density of the fluid A to be 10 and the density of the fluid B to be 5 here. Viscosity for the fluid A is equal to μ = 0.2.
Viscosity for the fluid B is equal to μ = 0.09. The surface tension is 0.05. In Figure 4.8, we plot the evolution of the bubble at time T = 0,100512,170512,200512,240512,270512,
280
512,320512 .
Figure 4.8 : The results for the above examples.
4.9 Moving Contact Line (acute)
In this case, we compute the behavior of a bubble on the floor under the influence of the friction force at the contact line without gravity. Our computational domain Ω = 0,1 × [0,1], using 128× 128 grid points with ∆t = 0.025∆x. We set the static equilibrium contact angle θ0 = 90. The fluid is at a standstill initially. The interface is a part of circle. The circle is centered at (0.5, -0.15) with radius 0.3, therefore the initial contact angle θ = 60. The viscosity is 0.001 in both fluids, the interfacial tension is 0.03. In Figure 4.9, we plot the evolution of the bubble at time T = 0,25,45,65,85, 2 .
Figure 4.9 : The results for the above examples.
4.10 Moving Contact Line (obtuse)
As 4.9, we compute the behavior of a bubble on the floor under the influence of the friction force at the contact line without gravity. The fluid is at a standstill initially. To set the static equilibrium contact angle θ0 = 90. The interface is a part of circle. Here, we set the circle is centered at (0.5, 0.1) with radius 0.2, therefore the initial contact angle θ = 120. All other parameters are as 4.9. In Figure 4.10, we plot the evolution of the bubble at time T = 0,25,45,65,85, 2 .
Figure 4.10 : The results for the above examples.
5 Concluding Remarks
We have introduced the numerical scheme of how to use level set methods to describe interface and displayed many numerical experiments of two phase flows. These results show that the efficiency of level set methods describe interface with large topological change.
But, there are many details need to be improved and there are many problems what we can challenge. For example, the re-initialization of level set function is essential if we use level set methods. Because it uses iteration, we must cost much time.
Therefore, it’s a good problem that how to increase our work efficiency.
The problem of surfactant is also very practical, many authors have researched this problem carefully. But most of them are based on front tracking where the free interfaces are explicitly tracked. It’s a challenge that how to solve the problem of surfactant using level set methods [10].
We need to solve the equation defined on the interface Γ 𝑓t + u ∙ ∇𝑓 − 𝐧 ∙ (∇u𝐧)𝑓 = 𝐷𝑠∇𝑠2𝑓
Where f is the consistency of surfactant, n is the normal vector to Γ directed towards outside, 𝐷𝑠 is the surfactant diffusivity, ∇𝑠= (𝐼 − 𝐧⨂𝐧)∇ is the surface gradient, u is velocity.
Since the level set method doesn’t like the front-tracking method which can control the information on interface completely. Therefore we should combine level set method with the other numerical methods or extend the equation to a band near the interface instead of restraining on the interface.
It’s a challenge to extend this method to three-dimension. But in order to describe physical phenomenon accurately, it’s a meaningful work we must do. As it should be, to maintain the efficiency and accuracy when extending this method into three-dimension is very important.
References
[1] Y. C. Chang, T.Y. Hou, B. Merriman, S. Osher, “A Level Set Formulation of Eulerian Interface Capturing Methods for Incompressible Fluid Flows.”Journal of Computational Physics 124, 449-464 (1996).
[2] Hauxiong Huang, Dong Liang, and Brian Wetton, “Computation of a Moving Drop/ Bubble on a Solid Surface Using a Front-Tracking Method.”August 23, 2004.
[3] Michael Renardy, Yuriko Renardy, and Jie Li, “Numerical Simulation of Moving Contact Line Problems Using a Volume-of-Fluid Method. ” Journal of Computational Physics 171, 243-263 (2001).
[4] S. Osher, R. Fedkiw, “Level Set Methods and Dynamic Implicit Surfaces.”
Springer, 2003.
[5] A. Harten, B. Engquist, S. Osher, and S. Chakravarthy, “Uniformly High-order Accurate Essentially Non-Oscillatory Schemes III. ” J. Comput. Phys. 71, 231-303 (1987).
[6] C. W. Shu, and S. Osher, “ Efficient Implementation of Essentially Non Oscillatory Shock Capturing Schemes.”J. Comput. Phys. 77, 439-471 (1988).
[7] C. W. Shu, and S. Osher, “ Efficient Implementation of Essentially Non Oscillatory Shock Capturing Schemes II (two).”J. Comput. Phys. 83, 32-78 (1989).
[8] S. Osher, and J. Sethian, “Fronts Propagating with Curvature Dependent Speed:
Algorithms Based on Hamilton-Jacobi Formulations.” J. Comput. Phys. 79, 12-49 (1988).
[9] H. C. Tseng, “Numerical Methods and Applications for Immersed Interface Problems.”Master Thesis.
[10]J. J. Xu, Z. L. Li, J. Lowengrub, H. Zhao, “A Level-set Method for Interfacial Flows with Surfactant.”J. Comput. Phys. 212, 590-616 (2006).
[11]D. P. Peng, B. Merriman, S. Osher, H. K. Zhao, and M. J. Kang, “A PDE-Based Fast Local Level Set Method.” J. Comput. Phys. 155, 410-438 (1999).
[12]M. C. Lai, Y. H. Tseng, H. X. Huang, “An immersed boundary method for interfacial flows with insoluble surfactant.”submitted.