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 F_s .

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 + 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 ²+ do ²

Bo = dB ²+ do ²

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 F_s . By section 3.3, we can obtain the lengths of F , … , F_sF1 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 = F_s × 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,¹⁰⁰₁₂₈,²⁰⁰₁₂₈,³⁰⁰₁₂₈,⁴⁰⁰₁₂₈,⁵⁰⁰₁₂₈,⁶⁰⁰₁₂₈,⁷⁰⁰₁₂₈ .

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,²⁵⁰₅₁₂,⁵⁰⁰₅₁₂,⁷⁵⁰₅₁₂,¹⁰⁰⁰₅₁₂,

1250

512 ,¹⁵⁰⁰₅₁₂ ,¹⁷⁵⁰₅₁₂ .

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,₂₅₆⁹ ,₂₅₆¹⁸ ,₂₅₆²⁷ ,₂₅₆³⁶ ,₂₅₆⁴⁵ ,₂₅₆⁵⁴ ,₂₅₆⁶³ .

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,₂₅₆⁹ ,₂₅₆¹⁸ ,₂₅₆²⁷ ,₂₅₆³⁶ ,₂₅₆⁴⁵ ,₂₅₆⁵⁴ ,₂₅₆⁶³ .

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,¹²⁰₅₁₂,²⁷⁶₅₁₂,²⁸⁸₅₁₂,³⁰⁰₅₁₂,³¹²₅₁₂,

324

512,³⁶⁰₅₁₂ .

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,¹²⁰₅₁₂,²⁴⁰₅₁₂,²⁴⁶₅₁₂,²⁵²₅₁₂,²⁵⁸₅₁₂,²⁶⁴₅₁₂,³⁶⁰₅₁₂ .

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,³³⁰₅₁₂,³⁴⁰₅₁₂,³⁵⁰₅₁₂,³⁶⁰₅₁₂,

370

512,³⁸⁰₅₁₂,⁵⁰⁰₅₁₂ .

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,¹⁰⁰₅₁₂,¹⁷⁰₅₁₂,²⁰⁰₅₁₂,²⁴⁰₅₁₂,²⁷⁰₅₁₂,

280

512,³²⁰₅₁₂ .

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 θ₀ = 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,²₅,⁴₅,⁶₅,⁸₅, 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 θ₀ = 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,²₅,⁴₅,⁶₅,⁸₅, 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𝐧)𝑓 = 𝐷_𝑠∇_𝑠²𝑓

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.