7.5. Stokes’ law for a falling droplet with insoluble surfactants
In Section 5.3.2, a transient Stokes law is solved. We will now slightly modify this problem and exchange the rigid spherical object by a droplet. For simplicity, we still assume that the capillary number is small, i.e. the droplet does not deform when falling. Thereby, we do not have to consider moving meshes. The necessary modifications are first in the mesh class, so that also a droplet mesh is generated:
# Modification (besides renaming some interfaces): Add a droplet mesh
self.line(pSnorth,pSsouth,name="droplet_axisymm")
self.plane_surface("droplet_axisymm", "droplet_outside",name="droplet") # droplet domain
Then, in the constructor __init__ of the problem class, we add e.g. the viscosity for the droplet as well:
# Modifications: Also define viscosity of the droplet
self.droplet_radius = 0.25 * milli * meter # radius of the spherical object
self.droplet_density = 1200 * kilogram / meter ** 3 # density of the sphere
self.outer_density = 1000 * kilogram / meter ** 3 # density of the liquid
self.droplet_viscosity = 5 * milli * pascal * second # viscosity
self.outer_viscosity = 1 * milli * pascal * second # viscosity
self.surface_tension=50*milli*newton/meter # Also define a surface tension (although it does not matter on a static mesh)
In the define_problem() method, we remove the no-slip boundary condition and add the equations of the droplet:
# Do not add a no slip
# eqs += DirichletBC(velocity_x=0, velocity_y=0) @ "droplet_outside" # and no-slip on the object
# droplet equations
inertia_correction = self.droplet_density * vector([0, 1]) * partial_t(U)
deqs=MeshFileOutput()
deqs+=NavierStokesEquations(dynamic_viscosity=self.droplet_viscosity,mass_density=self.droplet_density,bulkforce=inertia_correction)
deqs+=DirichletBC(velocity_x=0)@"droplet_axisymm"
deqs+=NavierStokesFreeSurface(surface_tension=self.surface_tension,static_interface=True)@"droplet_outside"
deqs+=ConnectVelocityAtInterface()@"droplet_outside"
deqs += DragContribution(U) @ "droplet_outside" # Also add the tractions from the inside to the drag
self.add_equations(deqs@"droplet")
if __name__ == "__main__":
with FallingDropletProblem() as problem:
problem.run(0.5*second,startstep=0.05*second,outstep=True) # solve and output
When we have static meshes (i.e. no equations for the mesh positions are added), we must add static_interface=True to the NavierStokesFreeSurface. With that, the action of the Lagrange multiplier of the kinematic boundary condition (6.5) will be added to the velocity, not the mesh motion. Since the latter is not allowed to move, only an adjustment of the velocity can guarantee the kinematic boundary condition to hold. Thereby, the kinematic boundary condition is effectively replaced by a zero normal flow condition, i.e. (4.15).
The results of this modifications are shown in Fig. 7.5. The value of the surface tension does not matter here, since on static meshes, the droplet will remain perfectly spherical, i.e. corresponding to a zero capillary number. However, surface tension gradients can still drive Marangoni flow, as we will see next.
To obtain surface tension gradients, we add an insoluble surfactant to the droplet-outside interface. The corresponding transport equation for a surfactant interface concentration \(\Gamma\) with surface diffusivity \(D_S\) reads
In pyoomph, it is again trivial to implement this:
from falling_droplet import *
from pyoomph.expressions.phys_consts import *
class SurfactantTransportEquation(InterfaceEquations):
def __init__(self):
super(SurfactantTransportEquation, self).__init__()
self.D=1e-9*meter**2/second # diffusivity
def define_fields(self):
self.define_scalar_field("Gamma","C2",testscale=scale_factor("temporal")/scale_factor("Gamma"))
def define_residuals(self):
u=var("velocity") # velocity at the interface
G,Gtest=var_and_test("Gamma")
self.add_residual(weak(partial_t(G)+div(u*G),Gtest))
self.add_residual(weak(self.D*grad(G), grad(Gtest)))
Again, if we solve it on an interface, i.e. a manifold with co-dimension, div() and grad() will expand to their surface counterparts as discussed in Section 4.3.2.
Warning
This transport equation is only valid without mass transfer. When mass transfer is considered, the normal interface motion would not coincide with the normal interface velocity. Thereby, this equation must be changed.
The only thing we have to do is adding this equation to the interface, setting a suitable initial condition and scaling and modifying the surface tension so that it depends on \(\Gamma\). The simplest surface tension relation with surfactants is just \(\sigma(\Gamma)=\sigma_0-RT\Gamma\) with the gas constant \(R\) and temperature \(T\). Instead of modifying the problem class directly, we just add the additional_equations and modify the surface_tension in the run script of the simulation:
if __name__ == "__main__":
with FallingDropletProblem() as problem:
Gamma0=1*micro*mol/meter**2
problem.set_scaling(Gamma=Gamma0)
add_ieqs=SurfactantTransportEquation()+InitialCondition(Gamma=Gamma0)+MeshFileOutput()
problem.additional_equations+=add_ieqs@"droplet/droplet_outside"
T=20*celsius
problem.surface_tension=50*milli*newton/meter-gas_constant*T*var("Gamma")
problem.run(0.5*second,startstep=0.05*second,outstep=True) # solve and output
The surfactants get advected to the top of the droplet and hamper the flow in the droplet, cf. Fig. 7.5.
Fig. 7.5 Droplet falling due due gravitiy without (left) and with surfactants (right).