Two-Phase Flow: Level Set and Surface Tension

println("Julia ", VERSION)

Julia 1.12.6

Overview

A two-phase flow carries two immiscible fluids separated by an interface — an air bubble in water, an oil droplet in a carrier, a rising gas slug. Rasmah models it with a level-set method: a signed-distance field φ marks which phase a point is in (φ < 0 inside the droplet, φ > 0 outside), and the continuum surface force (CSF) turns the interfacial surface tension into a body force f_σ = −σ κ ∇φ δ(φ) concentrated on the interface.

Here σ is the surface-tension coefficient, κ = ∇·(∇φ/|∇φ|) the interface curvature, and δ a smoothed Dirac delta localising the force to the interface.

The Laplace pressure jump

The classic test for any surface-tension model is a static droplet: with no flow, the two momentum terms −∇p and f_σ must balance, so the pressure jumps across the interface by exactly

\[\Delta p = p_{\text{in}} - p_{\text{out}} = \frac{\sigma}{R}\]

for a 2D circle of radius R (Laplace–Young; the 3D sphere is 2σ/R). This is an exact prediction, so it’s the perfect verification.

Building the droplet problem

droplet_levelset(xc, R) gives the analytic signed distance of a circle, and two_phase_operator assembles the (u, p) operator from the two fluids’ densities/viscosities and the surface tension. The droplet is stationary, so u = 0 on the boundary and the convective term vanishes:

using Rasmah
using LinearAlgebra

function square(nx, ny)
    nv = (nx + 1) * (ny + 1)
    verts = Matrix{Float64}(undef, 3, nv)
    for j in 0:ny, i in 0:nx
        idx = i + j * (nx + 1) + 1
        verts[1, idx] = i / nx; verts[2, idx] = j / ny
    end
    faces = Matrix{Int}(undef, 3, 2nx * ny)
    k = 0
    for j in 0:(ny - 1), i in 0:(nx - 1)
        a = i + j * (nx + 1) + 1; b = a + 1; c = a + nx + 1; d = c + 1
        k += 1; faces[:, k] = [a, b, c]
        k += 1; faces[:, k] = [b, d, c]
    end
    TriangleMesh(verts, faces)
end

m  = square(24, 24)
U  = FESpace(m, ReferenceFE(2, 2); valuetype = VectorValue{2})   # velocity
P  = FESpace(m, ReferenceFE(2, 1))                               # pressure

R, σ = 0.2, 1.0
op = two_phase_operator(U, P, 1.0, 1.0, 0.1, 0.1, σ, droplet_levelset([0.5, 0.5], R))

FEOperator{MultiFieldFESpace{Tuple{FESpace{TriangleMesh{Matrix{Float64}, Matrix{Int64}, Nothing}, ReferenceFE, VectorValue{2}, Float64}, FESpace{TriangleMesh{Matrix{Float64}, Matrix{Int64}, Nothing}, ReferenceFE, ScalarValue, Float64}}, BlockStyle, Tuple{Matrix{Int64}, Matrix{Int64}}}, MultiFieldFESpace{Tuple{FESpace{TriangleMesh{Matrix{Float64}, Matrix{Int64}, Nothing}, ReferenceFE, VectorValue{2}, Float64}, FESpace{TriangleMesh{Matrix{Float64}, Matrix{Int64}, Nothing}, ReferenceFE, ScalarValue, Float64}}, BlockStyle, Tuple{Matrix{Int64}, Matrix{Int64}}}, Vector{Rasmah.ResidualTerm}, Nothing}(MultiFieldFESpace{Tuple{FESpace{TriangleMesh{Matrix{Float64}, Matrix{Int64}, Nothing}, ReferenceFE, VectorValue{2}, Float64}, FESpace{TriangleMesh{Matrix{Float64}, Matrix{Int64}, Nothing}, ReferenceFE, ScalarValue, Float64}}, BlockStyle, Tuple{Matrix{Int64}, Matrix{Int64}}}((FESpace{TriangleMesh{Matrix{Float64}, Matrix{Int64}, Nothing}, ReferenceFE, VectorValue{2}, Float64}(TriangleMesh{Matrix{Float64}, Matrix{Int64}, Nothing}([0.0 0.041666666666666664 … 0.9583333333333334 1.0; 0.0 0.0 … 1.0 1.0; 0.0 5.0e-324 … 1.83987082e-315 0.0], [1 2 … 599 600; 2 27 … 600 625; 26 26 … 624 624], nothing), ReferenceFE(2, :simplex, 2), VectorValue{2}, [1 2 … 599 600; 2 27 … 600 625; … ; 627 628 … 2397 2399; 628 630 … 2399 2401], [1 3 … 1197 1199; 2 4 … 1198 1200; … ; 1255 1259 … 4797 4801; 1256 1260 … 4798 4802], [0.0 0.041666666666666664 … 1.0 0.9791666666666667; 0.0 0.0 … 0.9791666666666667 1.0], 2401), FESpace{TriangleMesh{Matrix{Float64}, Matrix{Int64}, Nothing}, ReferenceFE, ScalarValue, Float64}(TriangleMesh{Matrix{Float64}, Matrix{Int64}, Nothing}([0.0 0.041666666666666664 … 0.9583333333333334 1.0; 0.0 0.0 … 1.0 1.0; 0.0 5.0e-324 … 1.83987082e-315 0.0], [1 2 … 599 600; 2 27 … 600 625; 26 26 … 624 624], nothing), ReferenceFE(2, :simplex, 1), ScalarValue, [1 2 … 599 600; 2 27 … 600 625; 26 26 … 624 624], [1 2 … 599 600; 2 27 … 600 625; 26 26 … 624 624], [0.0 0.041666666666666664 … 0.9583333333333334 1.0; 0.0 0.0 … 1.0 1.0], 625)), BlockStyle(), [0, 4802, 5427], Int64[], 5427, ([1 3 … 1197 1199; 2 4 … 1198 1200; … ; 1255 1259 … 4797 4801; 1256 1260 … 4798 4802], [4803 4804 … 5401 5402; 4804 4829 … 5402 5427; 4828 4828 … 5426 5426])), MultiFieldFESpace{Tuple{FESpace{TriangleMesh{Matrix{Float64}, Matrix{Int64}, Nothing}, ReferenceFE, VectorValue{2}, Float64}, FESpace{TriangleMesh{Matrix{Float64}, Matrix{Int64}, Nothing}, ReferenceFE, ScalarValue, Float64}}, BlockStyle, Tuple{Matrix{Int64}, Matrix{Int64}}}((FESpace{TriangleMesh{Matrix{Float64}, Matrix{Int64}, Nothing}, ReferenceFE, VectorValue{2}, Float64}(TriangleMesh{Matrix{Float64}, Matrix{Int64}, Nothing}([0.0 0.041666666666666664 … 0.9583333333333334 1.0; 0.0 0.0 … 1.0 1.0; 0.0 5.0e-324 … 1.83987082e-315 0.0], [1 2 … 599 600; 2 27 … 600 625; 26 26 … 624 624], nothing), ReferenceFE(2, :simplex, 2), VectorValue{2}, [1 2 … 599 600; 2 27 … 600 625; … ; 627 628 … 2397 2399; 628 630 … 2399 2401], [1 3 … 1197 1199; 2 4 … 1198 1200; … ; 1255 1259 … 4797 4801; 1256 1260 … 4798 4802], [0.0 0.041666666666666664 … 1.0 0.9791666666666667; 0.0 0.0 … 0.9791666666666667 1.0], 2401), FESpace{TriangleMesh{Matrix{Float64}, Matrix{Int64}, Nothing}, ReferenceFE, ScalarValue, Float64}(TriangleMesh{Matrix{Float64}, Matrix{Int64}, Nothing}([0.0 0.041666666666666664 … 0.9583333333333334 1.0; 0.0 0.0 … 1.0 1.0; 0.0 5.0e-324 … 1.83987082e-315 0.0], [1 2 … 599 600; 2 27 … 600 625; 26 26 … 624 624], nothing), ReferenceFE(2, :simplex, 1), ScalarValue, [1 2 … 599 600; 2 27 … 600 625; 26 26 … 624 624], [1 2 … 599 600; 2 27 … 600 625; 26 26 … 624 624], [0.0 0.041666666666666664 … 0.9583333333333334 1.0; 0.0 0.0 … 1.0 1.0], 625)), BlockStyle(), [0, 4802, 5427], Int64[], 5427, ([1 3 … 1197 1199; 2 4 … 1198 1200; … ; 1255 1259 … 4797 4801; 1256 1260 … 4798 4802], [4803 4804 … 5401 5402; 4804 4829 … 5402 5427; 4828 4828 … 5426 5426])), Rasmah.ResidualTerm[Rasmah.ResidualTerm{Rasmah.var"#two_phase_residual##0#two_phase_residual##1"{Float64, Float64, Float64, Float64, Float64, Rasmah.var"#droplet_levelset##0#droplet_levelset##1"{Vector{Float64}, Float64}, Float64}}(1, Rasmah.var"#two_phase_residual##0#two_phase_residual##1"{Float64, Float64, Float64, Float64, Float64, Rasmah.var"#droplet_levelset##0#droplet_levelset##1"{Vector{Float64}, Float64}, Float64}(1.0, 1.0, 0.1, 0.1, 1.0, Rasmah.var"#droplet_levelset##0#droplet_levelset##1"{Vector{Float64}, Float64}([0.5, 0.5], 0.2), 0.01)), Rasmah.ResidualTerm{Rasmah.var"#two_phase_residual##2#two_phase_residual##3"}(2, Rasmah.var"#two_phase_residual##2#two_phase_residual##3"())], nothing)

Fix u = 0 on the boundary and pin one pressure degree of freedom (to set the pressure level), then Newton-solve the static droplet:

walls = vcat(dofs(U, BoundingBox([0.0, 0.0], [1.0, 0.0])),
             dofs(U, BoundingBox([0.0, 0.0], [0.0, 1.0])),
             dofs(U, BoundingBox([1.0, 0.0], [1.0, 1.0])),
             dofs(U, BoundingBox([0.0, 1.0], [1.0, 1.0])))
nU = ndofs(U)
fixed = vcat(walls, nU .+ [1])
u = solve(op, DofBC(fixed; g = zeros(length(fixed))), zeros(nU + ndofs(P)); method = :newton)
p = u[(nU + 1):end]

625-element Vector{Float64}:
  0.0
  0.0001582111883969454
  0.0005184658151742731
  0.0005770999035259059
  0.0002698684336812754
 -0.0003542207954214097
 -0.001176533452733736
 -0.002027550250854811
 -0.0027068981586691136
 -0.0030146936901211435
  ⋮
 -0.020167932580407427
 -0.019768929934626114
 -0.01949709338579158
 -0.019309719109108076
 -0.019219804603866707
 -0.019248548783790474
 -0.01939124104021701
 -0.019597955647172042
 -0.01965183392984351

Verifying the Laplace jump

The pressure dofs are the mesh vertices (P1), so we read p at the vertex nearest the droplet centre (inside) and one far from it (outside):

center = argmin([sum(abs2, m.vertices[1:2, i] .- [0.5, 0.5]) for i in 1:size(m.vertices, 2)])
corner = argmin([sum(abs2, m.vertices[1:2, i] .- [0.05, 0.05]) for i in 1:size(m.vertices, 2)])
println("p_in − p_out = ", round(p[center] - p[corner], digits = 4))
println("σ / R        = ", σ / R)

p_in − p_out = 4.9237
σ / R        = 5.0

The computed jump converges to σ/R as the mesh refines and the interface width shrinks. You can visualise the pressure field with vtk(m; field = p).

Next steps

Navier–Stokes — the single-phase limit (σ = 0) this model generalises.
Stokes flow — the creeping-flow saddle system under the two-phase operator.