Commit 3902aa49 authored by s.zitz's avatar s.zitz
Browse files

WIP: Sliding droplet interacting with defect

parent b1a13753
Pipeline #23277 failed with stage
in 51 minutes and 43 seconds
#=
First attempt to generate a set of data for the sliding of a droplet on patterned substrates.
I am going to measure four different pattern geometries with four different contact angle contrasts.
The four geometries are:
- A square
- A equilateral triangle
- An ellipse
- A circle
All have approximately the same area.
=#
using DrWatson
@quickactivate :JuThinFilm
using DataFrames, CSV, CodecZlib, Mmap, Images
# Define some constants
lx = 512
ly = 256
maxtime = 50000
θₑ = 1/9 # Contact angle on the substrate, expect on the patch
δ₀ = 4.0 # Slip length, larger better for geometry
T = Float64
# The center of the lattice
center = (lx÷2, ly÷2)
# Things to be measured
measureframe = DataFrame()
dropletframe = DataFrame()
contactlineframe = DataFrame()
# Geometric parameters for the patterns
boxside = 62
triangleside = 94.3
circlerad = 35
largehalfax, smallhalfax = 25, 49
# Some structs that contain mostly default values
parameters = JuThinFilm.params{T}(lx=lx, ly=ly, mt=maxtime)
collisionparameters = JuThinFilm.collisionparams{T}(δ=δ₀)
pressureparameters = JuThinFilm.pressurestats{T}(θ=θₑ)
dists = default_distribution(T, parameters.lx, parameters.ly)
lat = JuThinFilm.lattice{T}()
# Loop over three different droplet volumes
for R₀ in [55.8 69.75 83.7]
# Loop over the four different patch shapes
for shape in ["box" "triangle" "ellipse" "circle"]
# Loop over five different contact angle contrasts, last one is less wettable than the surrounding substrate
for δᵪ in [1/36 1/18 1/12 1/9 -1/18]
vol = π/3*R₀^3*(2+cospi(θₑ)*(1-cospi(θₑ))^2)
println("Sliding experiment with $shape pattern and wetting contrast $(rad2deg(π*δᵪ)) and initial volume $(round(vol,digits=3))")
# Generate the initial droplet
mom = singledroplet(T, lx, ly, R₀, θₑ, center)
diffsurf = []
energy = []
contactline = []
positions = [0, 0, 0]
MaxVel = []
θₙ = fill(θₑ,lx,ly)
angle_switch = true
# And the defect with the correct shape and contrast
# Forces needed for the simulation
Fpressure = default_force(T, lx, ly)
Fsum = default_force(T, lx, ly)
Fslip = default_force(T, lx, ly)
Fbody = default_force(T, lx, ly)
# Get the first equilibrium here
equilibrium!(dists.feq, mom.height, mom.velocity.x, mom.velocity.y, mom.pressure, lat.ci, lat.wi, parameters)
copy!(dists.ftemp,dists.feq)
# Start of the main time loop
for t in 1:parameters.mt
# Only change the substrate after the droplet has started sliding!
if findmax(mom.height, dims=1)[2][:1][1][1] > lx÷2+150 && t > 1000 && angle_switch
if shape == "box"
θₙ, P = boxpattern(lx, ly, pressureparameters, δₐ=-δᵪ, side=boxside)
elseif shape == "triangle"
θₙ, P = trianglepattern(lx, ly, pressureparameters, δₐ=-δᵪ, side=triangleside)
elseif shape == "ellipse"
θₙ, P = ellipsepattern(lx, ly, pressureparameters, δₐ=-δᵪ, a=largehalfax, b=smallhalfax)
elseif shape == "circle"
θₙ, P = ellipsepattern(lx, ly, pressureparameters, δₐ=-δᵪ, a=circlerad, b=circlerad)
end
# Changes the substrate pattern only once and not ever time step!
angle_switch = false
end
# Performe the measurements
if t%10 == 0
# Measurements of surface area and three phase contact line lenght
JuThinFilm.measuresurfaceandenergy!(diffsurf, energy, mom, θₙ)
JuThinFilm.measurecontactlinelengthandposition!(contactline, positions, t, mom)
JuThinFilm.measuremaxvelocity(MaxVel, mom)
end
if t%5000 == 0
# IO for progression
mass = round(sum(mom.height),digits=3)
println("time step $t and mass $mass")
end
# Compute forces
capillarypressure!(mom.pressure, padarray(mom.height, Pad(:circular, 1, 1)), θₙ, pressureparameters)
h∇p!(Fpressure.x, Fpressure.y, mom.height, padarray(mom.pressure, Pad(:circular, 1, 1)))
slippage!(Fslip.x, Fslip.y, mom.height, mom.velocity.x, mom.velocity.y, collisionparameters)
inclinedplane!(Fbody.x, Fbody.y, mom.height, t, 1e-4, true, false, tonset=100, tsmooth=1000)
sumforces!(Fsum, Fpressure, Fslip, Fbody)
# Get the equilibrium and then performe a collision.
equilibrium!(dists.feq, mom.height, mom.velocity.x, mom.velocity.y, mom.pressure, lat.ci, lat.wi, parameters)
BGKandStream!(dists, Fsum, collisionparameters, kind="WFM")
macroscopicmoments!(mom, dists)
end
# Save measurements and the last snapshot of the height field
angle_diff = rad2deg(δᵪ*π)
# All measurements saved here
measureframe["L_shape_$(shape)_contrast_$(angle_diff)_vol_$(R₀)"] = contactline # ./ contactline[1] but use it afterwards not already!
measureframe["S_shape_$(shape)_contrast_$(angle_diff)_vol_$(R₀)"] = diffsurf
measureframe["E_shape_$(shape)_contrast_$(angle_diff)_vol_$(R₀)"] = energy
contactlineframe["TPC_position_shape_$(shape)_contrast_$(angle_diff)_vol_$(R₀)"] = positions
measureframe["MaxVel_shape_$(shape)_contrast_$(angle_diff)_vol_$(R₀)"] = MaxVel
# Finial fluid configuration saved here
dropletframe["height_shape_$(shape)_contrast_$(angle_diff)_vol_$(R₀)"] = vec(mom.height)
dropletframe["velx_shape_$(shape)_contrast_$(angle_diff)_vol_$(R₀)"] = vec(mom.velocity.x)
dropletframe["vely_shape_$(shape)_contrast_$(angle_diff)_vol_$(R₀)"] = vec(mom.velocity.y)
dropletframe["pres_shape_$(shape)_contrast_$(angle_diff)_vol_$(R₀)"] = vec(mom.pressure)
end
end
end
savepath = "data/sims/dynamic_different_shape_angle/"
results_measures = "measurements_0.csv"
results_drop = "final_hdvars_0.csv"
results_TPC = "contactline_0.csv"
CSV.write(savepath * results_measures, measureframe)
CSV.write(savepath * results_drop, dropletframe)
CSV.write(savepath * results_TPC, contactlineframe)
zip1 = `gzip data/sims/dynamic_different_shape_angle/measurements_0.csv`
zip2 = `gzip data/sims/dynamic_different_shape_angle/final_hdvars_0.csv`
zip3 = `gzip data/sims/dynamic_different_shape_angle/contactline_0.csv`
run(zip1)
run(zip2)
run(zip3)
println("Simulations done results are saved in $savepath")
# TODO: Add a send email flag, to be notified if the simulation has finished
\ No newline at end of file
......@@ -6,29 +6,29 @@ What I want to see is that the droplet, deformes and then starts sliding.
using DrWatson
@quickactivate :JuThinFilm
using JuThinFilm, Plots
using Plots, Images
GR.inline("png") # Otherwise the animation will break...
palette(:tab10) # Hopefully the matplotlib colorscheme
# Some parameters to define the simulation
lx = 100
ly = 100
maxtime = 10000
θ₀ = 1/3
θₑ = 1/6
lx = 512
ly = 256
maxtime = 5000
θ₀ = 1/9
θₑ = 1/9
δ₀ = 4.0
h₀ = 25
h₀ = 55.8
T = Float64
center = (Int(lx/2), Int(ly/2))
mom = singledroplet(T, lx, ly, h₀, θ₀, center)
# Some parameter defining at the start
parameters = params{T}(lx=lx, ly=ly, mt=maxtime)
collisionparameters = collisionparams{T}(δ=δ₀, τ=1.0)
pressureparameters = pressurestats{T}(θ=θₑ)
dists = distributions{T}(para=parameters)
lat = lattice{T}()
parameters = JuThinFilm.params{T}(lx=lx, ly=ly, mt=maxtime)
collisionparameters = JuThinFilm.collisionparams{T}(δ=δ₀, τ=1.0)
pressureparameters = JuThinFilm.pressurestats{T}(θ=θₑ)
dists = JuThinFilm.distributions{T}(para=parameters)
lat = JuThinFilm.lattice{T}()
Fpressure = default_force(T, lx, ly)
Fsum = default_force(T, lx, ly)
......@@ -37,37 +37,41 @@ Fbody = default_force(T, lx, ly)
maxu = []
maxv = []
# Get the first equilibrium here
equilibrium!(dists, mom, lat, parameters, T)
copy!(dists.fout,dists.feq)
equilibrium!(dists.feq, mom.height, mom.velocity.x, mom.velocity.y, mom.pressure, lat.ci, lat.wi, parameters)
copy!(dists.ftemp,dists.feq)
# Now move to the time loop
anim = @animate for t in 1:parameters.mt
if t%500 == 0
if t%5000 == 0
println("time step $t")
end
push!(maxu, maximum(mom.velocity.x))
push!(maxv, maximum(mom.velocity.y))
l = @layout [a b; c{0.2h}; d{0.2h}]
u = extrema(mom.velocity.x) # Yeah the velocity is not taken from the center of the droplet.
v = extrema(mom.velocity.y)
uₘ = maximum(abs.(u))
vₘ = maximum(abs.(v))
push!(maxu,max(vₘ, uₘ))
l = @layout [a b; c{0.2h}]
p1 = surface(mom.height, c=:viridis)
p2 = contourf!(mom.height, aspect_ratio=1, c=:viridis)
p2 = contourf(mom.height, aspect_ratio=1, c=:viridis)
p3 = plot(maxu, xlabel = "Δt", ylabel = "u", lw=3, label="max(u)")
p4 = plot(maxv, xlabel = "Δt", ylabel = "v", lw=3, label="max(v)")
p = plot(p1, p2, p3, p4, layout = l)
plot!(p[1], zlims=(0, 14))
# p4 = plot(maxv, xlabel = "Δt", ylabel = "v", lw=3, label="max(v)")
p = plot(p1, p2, p3, layout = l)
plot!(p[1], zlims=(0, 10))
plot!(p[2])
plot!(p[3])
copy!(dists.ftemp,dists.fout)
# Compute forces
JuThinFilm.kernel_pressure_2!(mom.pressure, mom.height, pressureparameters, T)
filmpressure!(Fpressure, mom, T)
slippage!(Fslip, mom, collisionparameters, T)
constantforce!(Fbody, mom, t, 1e-4, true, false, T)
sumforces!(Fsum, Fpressure, Fslip, Fbody)
capillarypressure!(mom.pressure, padarray(mom.height, Pad(:circular, 1, 1)), pressureparameters)
h∇p!(Fpressure.x, Fpressure.y, mom.height, padarray(mom.pressure, Pad(:circular, 1, 1)))
slippage!(Fslip.x, Fslip.y, mom.height, mom.velocity.x, mom.velocity.y, collisionparameters)
# This set of parameters does work for the sliding droplet :)
inclinedplane!(Fbody.x, Fbody.y, mom.height, t, 1e-4, true, false, tonset=100, tsmooth=1000)
sumforces!(Fsum, Fpressure, Fslip, Fbody) #
# Get the equilibrium and then performe a collision.
equilibrium!(dists, mom, lat, parameters, T)
JuThinFilm.collidestreamBGK_periodicGUO!(dists, Fsum, lat, collisionparameters, T)
JuThinFilm.macroscopicmoments!(mom, dists, Fsum)
equilibrium!(dists.feq, mom.height, mom.velocity.x, mom.velocity.y, mom.pressure, lat.ci, lat.wi, parameters)
BGKandStream!(dists, Fsum, collisionparameters, kind="WFM")
macroscopicmoments!(mom, dists)
end
gif(anim, "animations/sliding_drop.gif", fps = 60)
\ No newline at end of file
gif(anim, "animations/sliding_drop.gif", fps = 500)
\ No newline at end of file
......@@ -355,7 +355,7 @@ function trianglepattern(lx, ly, stats::JuThinFilm.pressurestats; center=(lx÷2,
θₙ = ones(T, lx, ly)
height = sqrt(3)/2*side
xcoords = [center[1]-side/2 center[1]+side/2 center[1]]
ycoords = [center[2]-height/3 center[2]-height/3 center[2]+2*height/3]
ycoords = [center[2]+height/3 center[2]+height/3 center[2]-2*height/3]
# Check if the vertices are on the grid
for i in 1:length(xcoords)
if xcoords[i] < 1 || xcoords[i] > lx
......
......@@ -190,6 +190,23 @@ function measurecontactlinelength!(perimeter, mom)
push!(perimeter, sol)
end
function measurecontactlinelengthandposition!(perimeter, position, t, mom)
len, wid = size(mom.height)
height = copy(mom.height)
height[height .<= 0.06] .= 0
sol = 0
for i in 2:len-1
for j in 2:wid-1
# Basically if one or more of my neighbors is dry I have to be contactline!
if (height[i+1, j] == 0 && height[i, j] > 0) || (height[i-1, j] == 0 && height[i, j] > 0) ||(height[i, j+1] == 0 && height[i, j] > 0) ||(height[i, j-1] == 0 && height[i, j] > 0)
append!(position, (i, j, t))
sol += 1
end
end
end
push!(perimeter, sol)
end
"""
measuremaxheight!(maxheight, mom)
......@@ -265,6 +282,15 @@ function measurecoxvoinov!(coxvoi, mom, col, pstats)
push!(coxvoi, coxvoinovtuple)
end
function measuremaxvelocity(vel, mom)
u = extrema(mom.velocity.x)
v = extrema(mom.velocity.y)
uₘ = maximum(abs.(u))
vₘ = maximum(abs.(v))
realmax = max(vₘ, uₘ)
push!(vel, realmax)
end
#=========================================== Small helpers ================================================#
function ∇h_periodic!(resx, resy, arg)
# works as well :)
......
......@@ -152,7 +152,7 @@ function dummyrun()
h∇p!(Fpressure.x, Fpressure.y, mom.height, padarray(mom.pressure, Pad(:circular, 1, 1)))
slippage!(Fslip.x, Fslip.y, mom.height, mom.velocity.x, mom.velocity.y, collisionparameters)
inclinedplane!(Fbody.x, Fbody.y, mom.height, t, 1e-4, true, false)
sumforces!(Fsum, Fpressure, Fslip) #, Fbody
sumforces!(Fsum, Fpressure, Fslip, Fbody) #, Fbody
# Get the equilibrium and then performe a collision.
equilibrium!(dists.feq, mom.height, mom.velocity.x, mom.velocity.y, mom.pressure, lat.ci, lat.wi, parameters)
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment