mdrop.jl

# # Magnetic liquid droplet in magnetic field

# The behaviour of droplets under the action of magnetic fields of different configurations is an important issue in many domains, including microfluidics (Seeman et al. 2012), mechanics of tissues (Douezan et al. 2011; Frasca et al. 2014), studies of dynamic self-assembly (Timonen et al. 2013) and many others. After the successful synthesis of magnetic fluids in the late sixties by (Rosensweig 1985) an exciting story about droplets of magnetic fluid began. At first, the elongation of the droplets in an external field was observed and explained by Arhipenko et al. (1978). Later on, different droplet configurations were experimentally observed in the high-frequency rotating field, which was subject to our numerical study in [^1].

# Here I present the essential ingredients to reproduce the results of the paper with the new tools SurfaceTopology, LaplaceBIE and ElTopo (optional).

using LinearAlgebra
using GeometryTypes
using SurfaceTopology
using LaplaceBIE
using AbstractPlotting, GLMakie
## using ElTopo

# The calculation requires normal vector and vertex areas methods. Latter one gives the area of 1/3 of the vertex ring, so the sum is the area of the droplet.

function normals(vertices,topology)
    n = Point{3,Float64}[]
    for v in 1:length(vertices)
        s = Point(0,0,0)
        for (v1,v2) in EdgeRing(v,topology)
            s += cross(vertices[v1],vertices[v2])
        end
        normal = s ./ norm(s)
        push!(n,normal)
    end
    return n
end

function vertexareas(points,topology)
    vareas = zeros(Float64,length(points))
    for face in Faces(topology)
        v1,v2,v3 = face
        area = norm(cross(points[v3]-points[v1],points[v2]-points[v1])) /2
        vareas[v1] += area/3
        vareas[v2] += area/3
        vareas[v3] += area/3
    end
    return vareas
end

# These are tools to keep track if the calculation is sensible. Since we are working in an incompressible fluid, the volume needs to be constant, and the energy of the droplet decreases until equilibrium is reached.

function surfacevolume(points,topology)
    normal0 = [0,0,1]
    s = 0
    for face in Faces(topology)
        y1 = points[face[1]]
        y2 = points[face[2]]
        y3 = points[face[3]]

        normaly = cross(y2-y1,y3-y1)
        normaly /= norm(normaly)

        area = norm(cross(y2-y1,y3-y1))/2
        areaproj = dot(normaly,normal0)*area
        volume = dot(y1 + y2 + y3,normal0)/3*areaproj

        s += volume
    end
    return s
end

function energy(points,normals,faces,psi,mup,gammap,H0)
    vareas = vertexareas(points,faces)
    Area = sum(vareas)

    s = 0
    for xkey in 1:length(points)
        s += psi[xkey]*dot(H0,normals[xkey]) * vareas[xkey]
    end

    Es = gammap * Area
    Em = 1/8/pi * (1 - mup) * s

    return Es+Em
end

# For calculating the equilibrium of the droplet, we use a curvatureless algorithm for a vicious droplet developed by Zinchenko 1997. He found a way to calculate velocity generated by a surface tension without explicitly calculating the curvature, which makes it easier to implement and from some tests also more stable. The following method accepts surface defined by points, normals and faces (or topology) and surface force as well as surface tension γ and viscosity η (which only affects the scaling of time).

function stokesvelocity(points,normals,faces,forcen,etaP,gammap)
    vareas = vertexareas(points,faces)
    velocityn = zeros(Float64,length(points))

    for xkey in 1:length(points)

        x = points[xkey]
        nx = normals[xkey]
        fx = forcen[xkey]
                
        s = 0
        for ykey in 1:length(points)
            if ykey==xkey
                continue
            end

            y = points[ykey]
            ny = normals[ykey]
            fy = forcen[ykey]

            ### Need to check a missing 2
            s += vareas[ykey]*1 ./8/pi/etaP* dot(y-x,nx+ny)/norm(y-x)^3*(1-3*dot(y-x,nx)*dot(y-x,ny)/norm(y-x)^2) * gammap
            s += vareas[ykey]*1 ./8/pi/etaP* ( dot(nx,ny)/norm(x-y) + dot(nx,x -y)*dot(ny,x-y)/norm(x-y)^3 )*(fy - fx)
        end
        velocityn[xkey] = s
    end
    return velocityn
end

# Now having methods defined, we can include a sphere and proceed with calculation.

include("sphere.jl")
msh = unitsphere(2)
vertices, faces = msh.vertices, msh.faces

## Now let's do something fun. Visualize the process with Makie in real time.
x = Node(msh)
y = lift(x->x,x)

scene = Scene(show_axis=false)

wireframe!(scene,y,linewidth = 3f0)
mesh!(scene,y, color = :white, shading = false)

display(scene)

## Initial parameters
H0 = [4.,0.,0.]
etap = 1.
gammap = 1.
μ = 10.

t = 0.
Δt = 0.1
N = 100
volume0 = surfacevolume(vertices,faces)

record(scene, "mdrop.gif", 1:N) do i # for i in 1:N
    n = normals(vertices,faces) 

    psi = surfacepotential(vertices,n,faces,μ,H0)
    P∇ψ = tangentderivatives(vertices,n,faces,psi)
    Hn = normalderivatives(vertices,n,faces,P∇ψ,μ,H0)

    E = energy(vertices,n,faces,psi,μ,gammap,H0)
    rV = surfacevolume(vertices,faces)/volume0
    @show E,rV

    Ht = [norm(j) for j in P∇ψ]

    # The force generated (M⋅∇)H, which includes a jump of magnetization at the surface and force coming from the whole bulk.
    tensorn = μ*(μ-1)/8/pi * Hn.^2 + (μ-1)/8/pi * Ht.^2

    vn = stokesvelocity(vertices,n,faces,tensorn,etap,gammap)

    vertices .+= n .* vn * Δt

    msh = HomogenousMesh(vertices,faces) 
    
    ### ElTopo stabilization
    ## par = SurfTrack(allow_vertex_movement=true)
    ## msh = stabilize(msh,par)

    push!(x,msh)
    AbstractPlotting.force_update!()

    global vertices, faces = msh.vertices, msh.faces
    global t += Δt
end

# ![](mdrop.gif)

# [^1]: Erdmanis, J. & Kitenbergs, G. & Perzynski, R. & Cebers, A. (2017) Magnetic micro-droplet in rotating field: numerical simulation and comparison with experiment