using Plots

, Swalbe

21.4 s

Droplet coalescence

A fairly simple showcase is the dynamics of two coalescing sessile droplets. Of course we are not the first ones to take a look at this problem and in fact the literature concerning coalescing droplets is vast. For example references [1,2,3,4,5] study the coalescence dynamics of freely suspended droplets. While references [6,7] focus on coalescence of sessile droplets on partially wetting substrates with contact angles below $π / 2$ . We try to model the later, sessile droplets on a partially wetting substrate. Considering also the limits of the model we push towards higher contact angles to document when the small slope hypothesis fails for this problem.

The simulation is performed within the quasi two-dimensional (fast prototyping) system. Everything, besides the visualization with Plots.jl, is contained in the Swalbe.jl package. However for convenience we define two more functions. The first one contains the simulation, so we can run it with different parameters just by switching input arguments. The second one is for data analysis.

Experiment and analysis functions

The fluid dynamics simulation run_drop_coal()
Computation of the bridge height brdige_height()

One minor addition, for convenience we use a recipe that works quite similar to matplotlibs markevery function.

md"# Droplet coalescence

A fairly simple showcase is the dynamics of two coalescing sessile droplets.

Of course we are not the first ones to take a look at this problem and in fact the literature concerning coalescing droplets is vast.

For example references [[1](https://www.cambridge.org/core/services/aop-cambridge-core/content/view/0DE6817B486607109F5A32018951B5C0/S002211209900662Xa.pdf/coalescence_of_liquid_drops.pdf),[2](https://www.cambridge.org/core/services/aop-cambridge-core/content/view/CB56EE5EDDAFF2E9503CFEB0D2110662/S0022112003004646a.pdf/div-class-title-inviscid-coalescence-of-drops-div.pdf),[3](https://journals.aps.org/prl/pdf/10.1103/PhysRevLett.95.164503?casa_token=kdvf_DF7sK0AAAAA%3AN1Kw_LdxZMAfmGjHF5Cj39M7tugTp1frIl7HWV32tQw19igd43ha4Sk18H6fpJf2ZHruWHmQWK48),[4](https://arxiv.org/pdf/1406.5842.pdf),[5](https://journals.aps.org/prl/pdf/10.1103/PhysRevLett.106.114501?casa_token=rMfqSWjkqFUAAAAA%3AHWLZ6mGKdkI5nZj2PUNaHMDsnL_33bqD36iRgScHUpDQM-sIavUNOkwRobjobvm8fOQ9CStgJB_6)] study the coalescence dynamics of freely suspended droplets.

While references [[6](https://aip.scitation.org/doi/abs/10.1063/1.4824108?showFTTab=true&journalCode=phf),[7](https://journals.aps.org/prl/pdf/10.1103/PhysRevLett.111.144502?casa_token=I4xo5sbhnq4AAAAA%3A89YP7q4bHs2kJ36AKKgWzhrXRdCdNAUgcbQdXkEWFNHD_29yORY51rqk1na3HLToTi3xuBW4KWyp)] focus on coalescence of sessile droplets on partially wetting substrates with contact angles below $\pi/2$.

We try to model the later, sessile droplets on a partially wetting substrate.

Considering also the limits of the model we push towards higher contact angles to document when the small slope hypothesis fails for this problem.

The simulation is performed within the quasi two-dimensional (fast prototyping) system.

Everything, besides the visualization with [**Plots.jl**](https://github.com/JuliaPlots/Plots.jl), is contained in the [**Swalbe.jl**](https://github.com/Zitzeronion/Swalbe.jl) package.

However for convenience we define two more functions.

The first one contains the simulation, so we can run it with different parameters just by switching input arguments.

The second one is for data analysis.

## Experiment and analysis functions

- The fluid dynamics simulation `run_drop_coal()`

- Computation of the bridge height `brdige_height()`

One minor addition, for convenience we use a `recipe` that works quite similar to *matplotlibs* `markevery` function.

7.9 ms

@recipe function f(::Type{Val{:samplemarkers}}, x, y, z; step = 10)

n = length(y)

sx, sy = x[1:step:n], y[1:step:n]

# add an empty series with the correct type for legend markers

@series begin

seriestype := :path

markershape --> :auto

x := [Inf]

y := [Inf]

end

# add a series for the line

@series begin

primary := false # no legend entry

markershape := :none # ensure no markers

seriestype := :path

seriescolor := get(plotattributes, :seriescolor, :auto)

x := x

y := y

end

# return a series for the sampled markers

primary := false

seriestype := :scatter

markershape --> :auto

x := sx

y := sy

end

335 ms

run_drop_coal

run_drop_coal()

Simulation for droplet coalescence in quasi two dimensions.

"""

run_drop_coal()

Simulation for droplet coalescence in quasi two dimensions.

"""

function run_drop_coal(

sys::Swalbe.SysConst_1D;

r₁=115,

r₂=115,

θ₀=1/9,

verbos=true,

dump = 100,

fluid=zeros(sys.Tmax÷dump, sys.L)

)

println("Simulating droplet coalecense")

state = Swalbe.Sys(sys)

drop_cent = (sys.L/3, 2*sys.L/3)

state.height .= Swalbe.two_droplets(sys, r₁=r₁, r₂=r₂, θ₁=θ₀, θ₂=θ₀, center=drop_cent)

Swalbe.equilibrium!(state)

println("Starting the lattice Boltzmann time loop")

for t in 1:sys.Tmax

if t % sys.tdump == 0

mass = 0.0

mass = sum(state.height)

if verbos

println("Time step $t bridge height is $(round(state.height[Int(sys.L/2)], digits=3))")

end

Swalbe.filmpressure!(state, sys)

Swalbe.h∇p!(state)

Swalbe.slippage!(state, sys)

state.F .= -state.h∇p .- state.slip

Swalbe.equilibrium!(state)

Swalbe.BGKandStream!(state)

Swalbe.moments!(state)

96.4 ms

bridge_height

bridge_height()

Function that tracks the height or thickness at a single well defined lattice point.

"""

bridge_height()

Function that tracks the height or thickness at a single well defined lattice point.

"""

function bridge_height(data; t1=100, t_dur=size(data)[1], pos=size(data)[2]÷2)

bridge_h = []

for i in 1:t_dur

push!(bridge_h, data[i, pos])

end

return bridge_h

end

3.2 ms

false

# Need this for plotting purposes

begin

os = true

if Sys.iswindows()

os = false

elseif Sys.islinux()

os = true

end

84.0 μs

Initial condition

The plot below describes the initial condition of our numerical experiment. Two droplets that are barely touching with a contact angle θ of 20° (or π/9).

Due to surface tension (γ) we expect that these two droplets will coalesce into a single one. That is because a single droplet does have less surface area than two smaller ones and therefore it is an energetically favorable state.

md"## Initial condition

The plot below describes the initial condition of our numerical experiment.

Two droplets that are barely touching with a contact angle θ of 20° (or π/9).

Due to surface tension (γ) we expect that these two droplets will coalesce into a single one.

That is because a single droplet does have less surface area than two smaller ones and therefore it is an energetically favorable state."

396 μs

begin

# The initial configuration of the numerical experiment

rad = 500

h = Swalbe.two_droplets(Swalbe.SysConst_1D(L=1024), r₁=rad, r₂=rad)

x_ = 1:1024

p0 = plot(x_ , h, w=3, aspect_ratio=7, label="Initial conf.", xlabel="x [Δx]", ylabel="h(x)")

ylims!(0,40)

end

11.6 s

Experiment

We start the numerical experiment or simulation by simply calling the function that we defined before: run_drop_coal().

The data will be saved to an array that contains n-time steps and the whole height field. Something that is only possible with numerical simulations is to change parameters at will. We do this by changing the surface tension γ of our test liquid. Of course changing the contact angle θ or the kinematic viscosity ν would as well be possible.

md"## Experiment

We start the numerical experiment or simulation by simply calling the function that we defined before: `run_drop_coal()`.

The data will be saved to an array that contains n-time steps and the whole height field.

Something that is only possible with numerical simulations is to change parameters at will.

We do this by changing the surface tension γ of our test liquid.

Of course changing the contact angle θ or the kinematic viscosity ν would as well be possible."

405 μs

begin

# Sphere radius

sphere_rad = 500

# Different surface tension values

γs = [0.00008, 0.0003, 0.0006, 0.0008]

# Array that contains the simulation data

data_merge = zeros(20000, 1024, length(γs))

# Loop different surface tension values

for γ in enumerate(γs)

# System parameter, δ=50 can still considered small to medium slippage

sys = Swalbe.SysConst_1D(L=1024, Tmax=2000000, δ=50.0, γ=γ[2])

# The experiment

data_merge[:,:,γ[1]] = run_drop_coal(sys, r₁=sphere_rad, r₂=sphere_rad)

end

204 s

Data

Below is an example of the data that was generated during the simulation. We saved all Tmax time steps of all L lattice points in data_merge[:,:,γ] where we can specify the desired γ value in the last column of the array. Displayed is the data of the run with γ=8×10⁻⁵, or simply the first run.

A short word of caution. The lattice Boltzmann method is not perfect. On the plus side we have a relatively simple algorithm and the scales well. One of the downsides however is the fact that this method is derived from the Boltzmann equation. Therefore trying to make sense of these surface tension values is not too helpful. Interpreting results in terms of SI units is often not the best way. Personally I like to work with dimensionless numbers and relevant length and time scales.

Instead of adding SI units to the plots we normalize with characteristic quantities of the experiment.

md"### Data

Below is an example of the data that was generated during the simulation.

We saved all `Tmax` time steps of all `L` lattice points in `data_merge[:,:,γ]` where we can specify the desired `γ` value in the last column of the array.

Displayed is the data of the run with γ=8×10⁻⁵, or simply the first run.

A short word of caution.

[The lattice Boltzmann method](https://link.springer.com/book/10.1007/978-3-319-44649-3) is not perfect.

On the plus side we have a relatively simple algorithm and the scales well.

One of the downsides however is the fact that this method is derived from the [Boltzmann equation](https://en.wikipedia.org/wiki/Boltzmann_equation).

Therefore trying to make sense of these surface tension values is not too helpful.

Interpreting results in terms of **SI** units is often not the best way.

Personally I like to work with dimensionless numbers and relevant length and time scales.

Instead of adding **SI** units to the plots we normalize with characteristic quantities of the experiment."

634 μs

20000×1024 Matrix{Float64}:
 0.05       0.05       0.05       0.05       …  0.05       0.05       0.05
 0.05       0.05       0.05       0.05          0.05       0.05       0.05
 0.05       0.05       0.05       0.05          0.05       0.05       0.05
 0.05       0.05       0.05       0.05          0.05       0.05       0.05
 0.05       0.05       0.05       0.05          0.05       0.05       0.05
 0.05       0.05       0.05       0.05       …  0.05       0.05       0.05
 0.05       0.05       0.05       0.05          0.05       0.05       0.05
 ⋮                                           ⋱                        
 0.0499808  0.0499808  0.0499808  0.0499808     0.0499808  0.0499808  0.0499808
 0.0499808  0.0499808  0.0499808  0.0499808  …  0.0499808  0.0499808  0.0499808
 0.0499808  0.0499808  0.0499808  0.0499808     0.0499808  0.0499808  0.0499808
 0.0499808  0.0499808  0.0499808  0.0499808     0.0499808  0.0499808  0.0499808
 0.0499808  0.0499808  0.0499808  0.0499808     0.0499808  0.0499808  0.0499808
 0.0499808  0.0499808  0.0499808  0.0499808     0.0499808  0.0499808  0.0499808

# Data containing the γ=0.00008 simulation

data_merge[:,:,1]

68.6 ms

Since a plot says more than a few thousand columns, here is a plot of three different time steps. The initial kink between the two droplets is quickly smoothed out and the bridge (h₀) connecting the two droplets is growing.

md"Since a plot says more than a few thousand columns, here is a plot of three different time steps.

The initial kink between the two droplets is quickly smoothed out and the bridge (h₀) connecting the two droplets is growing."

271 μs

begin

# Time evolution for a single surface tension

lf = 8

p01 = plot(p0)

plot!(x_ ,data_merge[1000,:,1], w=3, alpha=0.6, label="t=100000Δt")

plot!(x_ ,data_merge[5000,:,1], w=3, alpha=0.6, label="t=500000Δt")

plot!(x_ ,data_merge[15000,:,1], w=3, alpha=0.6, label="t=1500000Δt",

xlims=(150,1024),

# legend=:bottomright,

legendfontsize = lf, # legend font size

tickfont = (14), # tick font and size

guidefont = (15), # label font and size

grid=:none)

end

2.3 ms

Under the bridge

The next point on the list is to measure the evolution of the bridge height h₀. We expect that the height will grow according to some power law, in fact given theoretical assumptions, see refs. [6,7], we hope to observe $h_{0} (t) \propto t^{α}$ with $α = 2 / 3$ .

Luckily this is a fairly easy task, since by definition the touching point of the droplets is at the center of the substrate.

md"### Under the bridge

The next point on the list is to measure the evolution of the bridge height h₀.

We expect that the height will grow according to some power law, in fact given theoretical assumptions, see refs. [6,7], we hope to observe $h_0(t) \propto t^{\alpha}$ with $\alpha = 2/3$.

Luckily this is a fairly easy task, since by definition the touching point of the droplets is at the center of the substrate."

356 μs

begin

bridges = zeros(20000, length(γs))

for i in 1:length(γs)

bridges[:, i] = bridge_height(data_merge[:,:,i])

end

490 ms

Now there are two data sets

Film thicknesses: $h (t, x)$
bridge heights: $h_{0} (t, γ)$

Especially the bridge heights should likely collapse to a single curve upon rescaling with a characteristic time scale $τ$ which we define as

$τ = \frac{μ R}{γ},$

where $R$ is the radius of the sphere from which the spherical cap has been cut, $γ$ is the surface tension and $μ$ is the liquid's viscosity. This time scale can be understood as a capillary time scale, see reference [8].

md"Now there are two data sets

- Film thicknesses: $h(t,x)$

- bridge heights: $h_0(t,γ)$

Especially the bridge heights should likely collapse to a single curve upon rescaling with a characteristic time scale $\tau$ which we define as

$$\tau = \frac{\mu R}{\gamma},$$

where $R$ is the radius of the sphere from which the spherical cap has been cut, $\gamma$ is the surface tension and $\mu$ is the liquid's viscosity.

This time scale can be understood as a capillary time scale, see reference [[8](https://journals.aps.org/pre/abstract/10.1103/PhysRevE.100.033313)]."

713 μs

τ(γ, R, μ)

Computes a characteristic time scale.

"""

τ(γ, R, μ)

Computes a characteristic time scale.

"""

function τ(γ; R=500, μ=1/6)

time_c = 0.0

time_c = R*μ/γ

return time_c

end

2.2 ms

Results

In the following we show the results collected. First being the evolution of the bridge height h₀(t). This information may not be too insightful, because the data is given in terms of lattice Boltzmann units (l.b.u.). To overcome this issue we have two natural scales. On the one hand a time scale defined by τ and on the other hand a length scale given by the spheres radii R.

τ ... time
R ... length

In the plot below we see that pure data spans four different curves, based on their respective surface tensions.

md"## Results

In the following we show the results collected.

First being the evolution of the bridge height h₀(t).

This information may not be too insightful, because the data is given in terms of lattice Boltzmann units (l.b.u.).

To overcome this issue we have two natural scales.

On the one hand a time scale defined by **τ** and on the other hand a length scale given by the spheres radii **R**.

- τ ... time

- R ... length

In the plot below we see that *pure* data spans four different curves, based on their respective surface tensions.

720 μs

begin

labels_gamma = ["γ=0.00008" "γ=0.0003" "γ=0.0006" "γ=0.0008"]

γ_h = [0.00008, 0.0003, 0.0006, 0.0008]

time_lbm = 100:100:2000000

p1 = plot(time_lbm,

bridges,

xlabel="t [Δt]",

ylabel="h₀",

label=labels_gamma,

xticks=([0:1000000:2000000;], ["0", "1×10⁶", "2×10⁶"]),

xlims=(0,2100000),

ylims=(0, 40),

w = 3, # line width

st = :samplemarkers, # some recipy stuff

step = 1800, # density of markers

marker = (8, :auto, 0.6), # marker size

legend=:bottomright,

legendfontsize = lf, # legend font size

tickfont = (14), # tick font and size

guidefont = (15), # label font and size

grid=:none

)

end

7.1 ms

In the box below the lattice Boltzmann time steps are normalized and nondimensionalized with τ.

md"In the box below the lattice Boltzmann time steps are normalized and nondimensionalized with τ."

240 μs