7.8. Beads on a string: a free surface driven by a polymer stress
The Rayleigh-Plateau filament of Section 7.6 was Newtonian, and it did what a Newtonian filament does: the neck thins ever faster and pinches off. If the liquid is a polymer solution instead, the outcome changes completely. The thread between two forming drops is a uniaxial extensional flow, which is exactly the flow in which the Oldroyd-B chains of Section 6.3.5 stretch without bound, and the axial stress they build up arrests the pinch-off. What is left is a beads-on-a-string structure: nearly spherical drops joined by threads that keep thinning exponentially instead of breaking [8].
This section combines the two ingredients – the moving mesh with reconstruction-based remeshing from the previous section, and the log-conformation viscoelasticity of Section 6.3.5 – to reproduce the corresponding verification model of the COMSOL Polymer Flow Module.
Lengths are scaled by the unperturbed filament radius \(R_0\), stresses by \(\sigma/R_0\) and time by the inertio-capillary time \(\tau=\sqrt{\rho R_0^3/\sigma}\), which leaves three dimensionless numbers,
i.e. the Ohnesorge number, the Deborah number and the solvent fraction of the total viscosity \(\eta_0=\eta_\mathrm{s}+\eta_\mathrm{p}\). The dimensionless solvent and polymer viscosities are then \(\beta\,\mathrm{Oh}\) and \((1-\beta)\mathrm{Oh}\), and the dimensionless relaxation time is numerically \(\mathrm{De}\). The reference case is
def __init__(self):
super().__init__()
self.Oh = 3.16 # Ohnesorge number
self.De = 94.9 # Deborah number, i.e. the dimensionless relaxation time
self.beta = 0.25 # solvent viscosity ratio
self.perturbation = 0.05 # amplitude of the initial radius perturbation
self.wavenumber = 0.5 # r(z,0) = 1 + perturbation*cos(wavenumber*z)
self.L = 8 * pi # domain length, i.e. two wavelengths
self.elements_per_radius = 3 # resolution across the *local* radius
self.max_elements_per_radius = 5 # upper bound on the resolution, i.e. lower bound on h
self.thinnest_thread = 0.02 # thinnest thread the mesh is allowed to resolve
self.interface_spline_points = 81
# Run control, all overridable from the command line with e.g. -P tend=100
self.tend, self.outstep, self.maxstep = 50.0, 5.0, 1.0
on a domain of two wavelengths, \(0\le z\le 8\pi\), starting from \(r(z,0)=1+0.05\cos(z/2)\) at rest. Note that these are plain attributes of the problem class, so each of them can be overridden on the command line with -P, which is used for the second parameter set at the end of this section.
The mesh class is the one of the previous section almost verbatim – gmsh size fields that follow the local radius on the interface and the distance to the interface on the axis, and a spline through get_boundary_coordinates() whenever we remesh. The only change is that the initial interface is not a straight line but the perturbed cosine profile. That size control is not a refinement here but a necessity: the thread ends up twenty-five times thinner than the beads that feed it, and a mesh which is uniform on the scale of a bead has no elements across the thread at all.
The equations are where the two chapters meet:
def define_problem(self):
self.add_mesh(FilamentMesh())
self.set_coordinate_system("axisymmetric")
# The Navier-Stokes viscosity is the *solvent* one; the polymer contributes its stress
# through the viscoelastic equations
eqs = NavierStokesEquations(mass_density=1, dynamic_viscosity=self.beta * self.Oh)
eqs += ViscoelasticEquations(model=OldroydB(), relaxation_time=self.De,
polymer_viscosity=(1 - self.beta) * self.Oh)
eqs += HyperelasticSmoothedMesh()
eqs += MeshFileOutput()
eqs += RemeshWhen(RemeshingOptions())
eqs += AxisymmetryBC() @ "axisymm"
# The reference imposes periodicity at the two ends. The initial perturbation has a maximum
# at both of them, so the solution is mirror-symmetric there and a symmetry plane - free
# slip for the flow, a mesh that may only slide radially - imposes the same thing.
eqs += DirichletBC(mesh_y=True, velocity_y=0) @ ["top", "bottom"]
# No traction condition of its own is needed for the polymer: the momentum residual is
# assembled in weak form, so the free surface sees the total stress, solvent plus polymer.
eqs += NavierStokesFreeSurface(surface_tension=1) @ "interface"
eqs += ExtremumObservables(min_r=var("mesh_x")) @ "interface"
eqs += TextFileOutput() @ "interface"
eqs += AssignZetaCoordinatesByEulerianCoordinate("x") @ "bottom"
eqs += AssignZetaCoordinatesByEulerianCoordinate("x") @ "top"
eqs += AssignZetaCoordinatesByEulerianCoordinate("y") @ "axisymm"
eqs += AssignZetaCoordinatesByArclength(sort_along_axis="y+") @ "interface"
self.add_equations(eqs @ "liquid")
Again, the viscosity handed to the NavierStokesEquations is the solvent viscosity only. The polymer enters exclusively through the ViscoelasticEquations, which add both the evolution equation for \(\log\mathbf{C}\) and the polymer stress in the momentum equation.
The two ends are treated as symmetry planes, while the reference uses a periodic flow condition to mimic an infinite filament. The initial perturbation has a maximum at both \(z=0\) and \(z=8\pi\), so the periodic solution is mirror-symmetric there anyway.
The run itself just steps through the output times, letting the adaptive time stepping choose the steps in between, and records the minimum radius through the ExtremumObservables as before.
RemeshWhen rebuilds the mesh whenever the elements have deformed too much, which over these 300 time units happens six times.
Fig. 7.8 Filament profiles at \(t=0,20,30,100\) and \(300\), to be compared with Fig. 1 of the reference model. The instability grows slowly, collapses into the beads-on-string structure between \(t\approx20\) and \(t\approx30\), and from then on merely drains the threads into the beads, which become spherical.
The quantity the reference reports is the minimum radius. It drops by a factor of four in the first thirty time units and then crosses over to a straight line in \(\log r\). That is the elasto-capillary regime: surface tension pulls on the thread, the polymer stress balances it, and the balance can only hold if the thread thins at the rate at which the polymer relaxes,
with the factor three coming from the uniaxial kinematics [21].
Fig. 7.9 The minimum radius, to be compared with Fig. 2 of the reference model, whose curve is the dashed line. Both leave the fast stage at \(\log_{10}r\approx-0.65\) around \(t=30\) and follow the asymptotic slope from there.
Fitting the computed curve for \(t>100\) gives a decay rate of \(3.65\cdot10^{-3}\) against the predicted \(1/(3\lambda)=3.51\cdot10^{-3}\), i.e. the thread thins \(4\,\%\) faster than the slender asymptote – unsurprisingly, since that asymptote neglects both the solvent viscosity and the drainage into the beads. The minimum radius ends at \(\log_{10}r=-1.085\) at \(t=300\), against \(-1.06\) read off the reference figure. This run has about 28 000 degrees of freedom and takes some ten minutes.
Lowering the viscosity and the elasticity to \(\mathrm{Oh}=0.4\), \(\mathrm{De}=0.8\) changes the picture: elasticity no longer holds the thread together on its own and inertia has time to act on it, so each thread grows a satellite drop before it breaks. Since every parameter is an attribute of the problem class, this second case needs no second script:
python beads_on_string.py -P Oh=0.4 De=0.8 tend=19.75 outstep=0.25 maxstep=0.25 \
elements_per_radius=5 max_elements_per_radius=8 thinnest_thread=0.005
The finer mesh is not optional – the reference notes the same requirement – and it is what makes this the expensive one of the two: the size field keeps refining as the thread thins, so the problem grows from 67 000 to nearly 200 000 degrees of freedom on the way to \(t=19.75\), and the run takes about a quarter of an hour.
Fig. 7.10 Compare Fig. 3 of the reference model: the elasticity-dominated case at \(t=300\) (left) keeps a smooth thread, while the low-Ohnesorge case (right) grows a satellite drop on each of the two threads.