PaperlayerView exact revision on arXiv

Source-linked AI summary

A discrete geometric approach for simulating the dynamics of thin viscous threads

Basile Audoly, Nicolas Clauvelin, Pierre-Thomas Brun, Miklós Bergou, Eitan Grinspun, Max Wardetzky

arXiv:1202.4971v2physics.flu-dyncs.CGmath.DGnlin.PS

TL;DR

Thin viscous-thread dynamics requires a robust treatment of coupled stretching, bending, twisting, inertia, and large rotations, especially for unsteady flows. The paper develops a discrete Lagrangian model using reduced centerline/spin coordinates, discrete geometric twist, and variational internal stress. It establishes consistency with the smooth equations and validates the method against reference coiling solutions.

  • Problem

    Unsteady thin-thread dynamics involves coupled deformation modes, nonlinear rotations, and numerically stiff equations, while existing approaches often address narrower cases.

  • Method

    The paper builds a fully discrete Lagrangian viscous-thread model with reduced centerline/spin coordinates, discrete-geometric twist, and variationally derived internal viscous stress.

  • Results

    The method is validated against reference solutions for steady coiling, with simulations closely matching the exact-value approximation when the thread radius is small.

  • Takeaways & Limitations

    The model supports robust and efficient simulation of unsteady thin viscous jets with combined inertia, stretching, bending, twisting, large rotations, and surface tension.

  • Takeaways & Limitations

    The linear implicit scheme is easier to use than a fully nonlinear implicit scheme but does not preserve conservation laws associated with problem symmetries.

Abstract

from arXiv · show

We present a numerical model for the dynamics of thin viscous threads based on a discrete, Lagrangian formulation of the smooth equations. The model makes use of a condensed set of coordinates, called the centerline/spin representation: the kinematical constraints linking the centerline's tangent to the orientation of the material frame is used to eliminate two out of three degrees of freedom associated with rotations. Based on a description of twist inspired from discrete differential geometry and from variational principles, we build a full-fledged discrete viscous thread model, which includes in particular a discrete representation of the internal viscous stress. Consistency of the discrete model with the classical, smooth equations is established formally in the limit of a vanishing discretization length. The discrete models lends itself naturally to numerical implementation. Our numerical method is validated against reference solutions for steady coiling. The method makes it possible to simulate the unsteady behavior of thin viscous jets in a robust and efficient way, including the combined effects of inertia, stretching, bending, twisting, large rotations and surface tension.

1. Introduction

Thin viscous-thread dynamics combines stretching, bending, twisting, inertia, and large rotations, creating nonlinear, stiff equations that are difficult to solve analytically and numerically. The paper addresses this with a fully discrete, geometrically controlled model designed for efficient simulation of both steady and unsteady behavior.

  • Context: Unsteady thin-thread flows can produce complex patterns when relative motion suppresses steady coiling, motivating robust simulations beyond steady solutions.The fluid-mechanical sewing-machine experiment generates more than ten patterns by varying surface velocity and fall height.
  • Context: Stretching, bending, and twisting are coupled by finite-rotation nonlinearities, while fourth-order spatial derivatives make the governing equations numerically stiff.These properties eliminate analytical solutions and make controlled spatial discretization important.
  • Related work: Thin-filament simulations benefit from dimensionally reduced equations because general free-boundary flow methods are inefficient when thickness is small relative to longitudinal length scales.Existing approaches include marker-and-cell variants, GENSMAC, and implicit projection schemes for viscous flows.
  • Related work: Unlike models neglecting selected deformation modes, the proposed method targets combined twist, bending, stretching, inertia, and large rotations in non-steady dynamics.Previous work includes steady rotating threads and viscous strings without bending or twisting.
  • Model scope: The model assumes incompressible viscous fluid with circular cross-sections that remain disks, excluding tubes and general non-axisymmetric cross-sections.Surface tension is assumed to round cross-sections quickly compared with the flow timescale.
  • Proposed approach: A fully discrete formulation defines strain rates, viscous forces, bending, and twist geometrically, allowing stable simulations with mesh sizes comparable to the smallest curvature radius.Discrete twist uses discrete differential geometry and holonomy, while the centerline/spin representation removes two rotational degrees of freedom.
  • Proposed approach: The paper establishes equivalence between formulations based on Kirchhoff equations and Rayleigh potentials, supporting a natural variational discretization.The Kirchhoff formulation is intuitive for fluid mechanicians, whereas the Rayleigh-potential formulation facilitates discretization.

2. Mathematical toolbox

The mathematical toolbox introduces moving-frame calculus for thin-thread mechanics, including Darboux vectors, covariant derivatives, and tangent-normal projections. These tools distinguish frame-aware derivatives from ordinary spatial derivatives and support the later geometric formulation.

  • Projections: Vectors are decomposed into tangent and perpendicular components using a projection operator defined by P ⊥(q, a) = a − (q · a) q.The tangent direction is taken as the third frame vector.
  • Moving frames: An orthonormal moving frame is described through a Darboux vector, which represents infinitesimal rotations and specializes to angular velocity in time or twist-curvature in space.The same definition applies when the continuous parameter is time or arc length.
  • Covariant calculus: The covariant derivative measures derivatives in the frame moving with the triad, unlike the ordinary derivative.It is introduced for vector fields associated with an orthonormal frame.
  • Projections: Covariant differentiation preserves tangent and normal subspaces, so the corresponding projections commute with the covariant time derivative.This compatibility does not hold for the regular derivative with the projection operators.

3. Smooth setting: a Lagrangian description of viscous threads

The paper reformulates thin-viscous-thread mechanics in Lagrangian variables so the geometric quantities can be discretized naturally. The formulation follows material particles, distinguishes reference and actual arc lengths, and uses a centerline with an associated tangent.

  • Lagrangian formulation: The Lagrangian reformulation is introduced to enable a geometric discretization of twist and a natural discrete viscous-thread model.The usual equations are expressed in Eulerian variables, whereas the paper develops a Lagrangian alternative.
  • Kinematics: Because the thread can stretch, the formulation distinguishes reference arc length S from actual arc length s.This distinction is central to describing deformation between reference and current configurations.
  • Reference configuration: An infinite circular cylinder of constant radius a0 is used as a convenient reference configuration, although the equations do not depend on that choice or radius.The reference configuration need not coincide with the thread at any particular time.
  • Constitutive assumptions: The fluid is treated as incompressible, with the reference-to-current mapping chosen to preserve volume for the model equations.Volume preservation simplifies the formulation but is not strictly required by incompressibility.
  • Kinematics: The coordinate S follows fixed fluid particles, and the centerline is represented by x(S, t) with material tangent T(S, t).A prime denotes spatial differentiation and a dot denotes time differentiation.

3.3. Incompressibility: radius and related quantities

The model describes thin viscous threads using incompressibility, material-frame kinematics, and a reduced centerline/spin representation. Twist and centerline motion are coupled through a geometric identity that supplies the twisting strain rate for the constitutive law.

  • 3.3. Incompressibility: radius and related quantities: Incompressibility preserves the reference volume through A(S,t) = A0 ℓ(S,t) and I(S,t) = I0 ℓ2(S,t).The thread is modeled with locally cylindrical cross-sections and a uniform cylindrical reference configuration.
  • Material-frame kinematics: The material frame tracks cross-section rotation, whose tangential component defines kinematical twist rather than Frenet–Serret torsion.Kinematical twist can remain nonzero for a planar centerline, including a straight twisted configuration.
  • Material-frame kinematics: The Kirchhoff kinematical hypothesis keeps material cross-sections perpendicular to the centerline tangent and couples frame rotations to centerline motion.This condition is justified by the thin-thread, shearless-flow limit and removes relative sliding between cross-sections.
  • Darboux vectors: The angular velocity is ω(S,t) = t(S,t) × ˙t(S,t) + v(S,t)t(S,t), combining tangent motion with axial spin velocity.The axial spin velocity is the tangential component of the temporal Darboux vector.
  • Centerline/spin representation: The geometric identity couples centerline motion and twist, allowing x(S,t) and v(S,t) to parameterize the thread while providing et = ˙τ for the constitutive law.The twisting strain rate depends on both the rotational degree of freedom and the centerline motion; the discrete model follows this discretization strategy.
  • Centerline/spin representation: The centerline/spin representation eliminates two rotational degrees of freedom by retaining the centerline and incremental spin velocity for isotropic cross-sections.It avoids tracking the absolute direction of the transverse material vectors while respecting the compatibility constraint.

4. Equivalence with Kirchhoff equations for a thin viscous thread

The centerline/spin formulation identifies viscous tension and internal bending–twisting moments through dissipation potentials, then shows these stresses reproduce the classical Kirchhoff equations. This equivalence provides a natural basis for discretization and efficient implementation.

  • Constitutive laws: The scalar coefficient ns represents tension resisting stretching, while the vector m represents the internal moment from twisting and bending.These coefficients depend on the actual motion and are identified as viscous stresses.
  • Constitutive laws: The constitutive laws agree with those derived from three-dimensional Stokes equations and can be rewritten using Eulerian strain rates.The formulation also clarifies the relation between viscous threads and elastic rods through the Rayleigh-Taylor analogy.
  • Viscous forces and moments: The viscous stress produces a net centerline force combining stretching, twisting, and bending contributions.The force is obtained from Rayleigh dissipation potentials, with endpoint terms restored in the discrete model.
  • Viscous forces and moments: The corresponding twisting moment is obtained by projecting the derivative of the internal moment along the material tangent.This expression matches the twisting-moment balance in the Kirchhoff formulation.
  • The centerline/spin formulation is shown to be equivalent to the classical Eulerian equations for thin viscous threads.The equivalence is established by matching force and twisting-moment expressions derived from dissipation potentials and Kirchhoff equations.

5. Space discretization: the discrete viscous thread model

The discrete model extends the centerline/spin representation with polygonal centerlines, segment-based geometry, parallel transport, and variationally derived viscous forces. Its construction preserves compatibility while reducing rotational degrees of freedom and defining discrete twist geometrically.

  • The discrete formulation uses centerline/spin coordinates, parallel transport for twist, and discrete dissipation potentials to derive equations of motion.The centerline compatibility condition eliminates two of the three rotational degrees of freedom.
  • Discrete geometry: The centerline is represented by n + 2 vertices, with forces, masses, and dynamics integrated at each vertex.The discrete viscous force is designed to converge to the smooth force as n →∞.
  • Discrete geometry: Segment vectors, lengths, tangents, velocities, and axial strain rates are defined as integrated discrete counterparts of smooth quantities.Vertex and segment indexing distinguish quantities associated with positions from those associated with segments.
  • Material quantities: The discrete model conserves segment volume and mass while reconstructing radius and cross-sectional area through incompressibility.Segment subdivision is the exception when adaptive meshing is used.
  • Discrete twist: Parallel transport is the unique compatible rotation of minimal angle mapping one segment tangent to the next.It is represented by a rotation about the binormal through the turning angle and corresponds to twist-less configurations.

5.5. Discrete twist

Discrete twist is defined by decomposing the finite rotation between adjacent material frames into parallel transport and an axial rotation. This construction recovers the smooth twist decomposition in the vanishing-discretization limit.

  • The rotation between adjacent material frames is decomposed into parallel transport and an axial rotation about the tangent.The axial angle τi is the discrete angle of twist across a vertex.
  • Parallel transport supplies the centerline-induced rotation, while τi supplies the additional rotation needed to match the two material frames.The angle is uniquely defined modulo 2π.
  • The discrete twist definition is consistent with the smooth relation π = K + τi t in the smooth limit.The finite rotations converge toward the infinitesimal rotation, curvature, and twist components.

5.6. Rate of change of twisting strain

The discrete twisting strain rate is obtained from the material derivative of the discrete twist angle. Its expression includes a holonomy term that couples centerline motion to twisting through changing parallel transport.

  • The twisting strain rate is defined as the material derivative of the discrete angle of twist.It is a spatially integrated counterpart of the smooth twisting strain rate.
  • The discrete rate formula is expressed using centerline velocities and spin velocities in the centerline/spin representation.The construction rewrites the twisting strain rate as a function of uj and vj.
  • Geometric decomposition: The rotation between material frames is analyzed using parallel transport and polar angles associated with adjacent segment tangents.These angles are represented in Figure 4 to illustrate the decomposition used in the rate calculation.
  • The holonomy term captures changes in parallel transport caused by centerline motion.This term couples centerline motion with the twisting mode and is identified as geometrical in origin.

5.7. Rate of change of bending strain

The discrete model represents bending and twisting strain rates through centerline geometry, vertex-based tangents, and spin velocities. Its operators are designed to reproduce the corresponding smooth quantities for real motions while retaining a compact centerline/spin representation.

  • Strain-rate decomposition: Discrete strain rates decompose the strain-rate vector into twisting and bending components using tangent and perpendicular projections.The decomposition uses a vertex-based tangent and angle-dependent normalizing functions.
  • Normalization choices: The bending normalization is chosen as hb(ϕi) = 1, while the twisting choice follows from the parallel-transport definition of discrete twist.Any hb converging to one as ϕi approaches zero is equivalent in the smooth limit.
  • Velocity dependence: Virtual vertex and spin velocities allow discrete viscous forces to be obtained as gradients of the dissipation potential.The discrete strain-rate operators are linear in the virtual velocity argument.
  • Centerline/spin representation: The centerline/spin representation stores vertex positions as generalized coordinates and spin angular velocities in the generalized velocity, without explicitly tracking material-frame orientation.The twisting mode remains coupled to centerline motion through generalized velocities.
  • Kinematic operators: The discrete operators Vi and Wi reproduce tangent time derivatives and angular velocities when evaluated on real motions.These operators extend the corresponding smooth kinematic operators to segments.

5.9. Dissipation potentials

The model constructs viscous internal forces from discrete stretching, twisting, and bending dissipation potentials. These potentials use centerline/spin velocities and produce a sparse, configuration-dependent dissipation matrix that is formally consistent with the smooth model.

  • Generalized velocities: The generalized velocity collects vertex linear velocities and segment spin angular velocities, while virtual velocities define the dissipation potentials.The real and virtual velocity vectors share the centerline/spin ordering.
  • Dissipation potentials: Discrete viscous internal forces are derived from dissipation potentials representing stretching, twisting, and bending contributions.Stretching is summed over segments, whereas twisting and bending are summed over interior vertices because their strain rates are defined on different entities.
  • Discrete moduli: The discrete moduli depend on viscosity, cross-sectional area, segment lengths, and Voronoi-cell lengths, but not on velocities.Their relation Bi/Ci = 3/2 matches the smooth model and reflects incompressibility.
  • Smooth-limit consistency: As n approaches infinity, the discrete model formally converges to the smooth equations of motion and their internal viscous-force expression.The convergence follows from consistency of the discrete dissipation potential and is checked numerically later.
  • Internal forces: Viscous forces and twisting moments are obtained by differentiating the dissipation potential with respect to vertex and spin velocities.These quantities enter the discrete equations of motion alongside external loading and mass terms.

5.11. Surface tension and other forces

Surface tension is modeled through a cylindrical segment representation and a discrete capillary energy based on lateral area. Differentiating this energy yields capillary forces that reproduce endpoint contributions and converge to the smooth forces under refinement.

  • Cylindrical representation: The implementation represents each thread segment as a cylinder and derives surface tension from its lateral area.This cylindrical representation is simpler than the truncated-cone formulation used previously.
  • Capillary energy: Discrete capillary forces are obtained as the negative gradient of a capillary energy with respect to vertex positions.The capillary energy depends on vertex positions but not on the twist degree of freedom.
  • Approximation: The thin-thread approximation assumes slowly varying radius and neglects longitudinal curvature relative to azimuthal curvature.The approximation is not suitable for Rayleigh-Taylor instability analysis when the critical wavelength is comparable to the radius.
  • Force effect: Capillary forces tend to shorten and compact the thread by bringing endpoints together and flattening curved centerline regions.The force interpretation follows from an energy proportional to lateral surface area.
  • Endpoint and interior forces: At interior vertices, contributions from adjacent segments nearly cancel, whereas terminal vertices receive only one contribution.This captures the endpoint Dirac contributions of the smooth model, and the discrete forces converge as n approaches infinity.

6. Time discretization, numerical implementation

The implementation assembles sparse dissipation matrices and advances the constrained thread dynamics with a linear implicit time-stepping scheme. The resulting update requires solving a symmetric positive-definite linear system, while the mesh and constraint dimensions may change during simulation.

  • Sparse assembly: Local strain-rate operators depend only on neighboring vertices or segments, producing sparse vectors and matrices that support efficient storage and manipulation.The dissipation matrix is quadratic in virtual velocity and assembled from these sparse operators.
  • Band structure: The dissipation matrices for stretching, twisting, and bending are symmetric and band-diagonal because their underlying linear forms are sparse.Their band structure is illustrated in Figure 6.
  • Changing discretization: The number of vertices and the dissipation-matrix dimension may vary between time steps as vertices are created, discarded, or affected by changing constraints.Fluid properties are stored on segments and may also vary with time when adaptation or coupled processes are used.
  • Kinematic constraints: Kinematical constraints are imposed by dispatching independent degrees of freedom through a matrix B and prescribed constrained velocities through B′.The number of independent unknowns can change when constraints are created or destroyed, such as at first ground contact.
  • Time stepping: A linear implicit scheme evaluates viscous forces implicitly in velocity and explicitly in position, requiring only a linear solver.The authors report that this choice combines stability with ease of implementation.
  • Linear solve: The velocity update system is symmetric positive definite for every positive time increment, enabling efficient and robust solvers.The system incorporates constraints, mass, dissipation, external loading, and the previous velocity.

7. Interaction of the thread with other bodies

The model handles thread interaction with containers and obstacles through kinematical constraints, with refined insertion and collision procedures improving numerical smoothness and reproducibility.

  • 7.1. Interaction with the container: Container injection prescribes the outlet velocity and clamps the thread by constraining the first two vertices and blocking their connecting segment’s rotation.The ejection velocity is Uc = Qc/Ac relative to the container, while the container’s motion contributes to the exiting-fluid velocity.
  • 7.1.1. Simple container model: The simple container model periodically frees vertices that pass the opening and inserts new segments assigned prescribed length, volume, mass, and surface tension.Two vertices remain inside the container, and newly created segments use ℓc, Acℓc, ρAcℓc, and γ.
  • 7.1.1. Simple container model: The simple insertion scheme produces small-amplitude, high-frequency oscillations because the effective fall height jumps whenever a new vertex is added.The discontinuity is controlled by the spacing between vertices inside the container and can impair acceleration convergence and reproducibility.
  • 7.1.2. Refined container model: The refined container model keeps the first two vertices fixed relative to the container and continuously assigns incoming volume ϵQc to the second segment before splitting it at the target volume.This makes fall height vary smoothly and yields smoothly convergent acceleration; all examples in Sections 8 and 9 use this model.
  • 7.2.1. ‘Capture and continue’ mode: The basic capture-and-continue collision scheme causes delayed momentum transfer, irregular obstacle penetration, and spurious acceleration fluctuations.The delayed constraint activation and penetration depth are both tied to the discrete time step.
  • 7.2.2. ‘Time roll-back’ mode: Time roll-back recomputes colliding steps with contact enforced at the step endpoint, suppressing momentum-transfer delay and unwanted landing rugosity.When many collisions occur within one step, linearization becomes inaccurate; adaptive shortening limits the method to at most one collision per step but increases simulation time.

8. Validation in a steady coiling geometry

Steady-coiling simulations are compared with numerical-continuation reference solutions while varying fall height and discretization. The model reproduces the reference coiling branches and validates the effects of bending, stretching, gravity, inertia, collisions, and surface tension.

  • 8. Validation setup: Validation uses steady coiling under gravity, comparing the discrete simulation with N. Ribe’s time-independent numerical-continuation solutions.The baseline parameters are µ = 0.2, ρ = 5 10^-4, g = 9.81, Ac = 6.44 10^-3, Qc = 3.96 10^-3, and γ = 0.
  • 8.1. Validation of bending, stretching, gravity, inertia and collisions: The simulation sweeps fall height by prescribed container motion and records the coiling radius R for comparison with the reference curve.The radius is measured at the floor contact point relative to the nozzle axis, after transients disappear, and averaged over several periods.
  • 8.1. Validation of bending, stretching, gravity, inertia and collisions: The simulated radius follows the reference solution closely until folds, then transitions through transient or bifurcated branches with smaller radii.Upward and downward sweeps show different transition heights, while the simulation reproduces the reference curve’s meandering shape.
  • 8. Validation setup: The validation case has Π1 = 7000, Π2 = 7, and Π3 = 0, corresponding to significant but not extreme stretching as the thread descends.The radius decreases by a factor of order 2 over the considered fall-height range.
  • 8.2. Analysis of convergence: Convergence is assessed by varying time step ϵ and segment length ℓc at fixed ratio, using refined container and floor models to handle collisions.The coiling radius is evaluated after the initial transient and averaged over multiple periods.
  • 8.3. Validation of surface tension: With γ = 10^-3 and Π3 = 10.3 10^-3, simulations agree well with the reference curve after surface tension substantially changes the coiling radius.The comparison uses the same parameters as the zero-surface-tension validation except for γ.

9. Discussion

The method is validated against steady coiling and applied to transient thread dynamics, including jumps between coiling branches and moving-belt patterns. Its linear implicit scheme is stable but requires careful discretization and does not conserve angular momentum.

  • Transient regimes: Transient simulations capture jumps between steady-coiling branches, including smaller-radius solutions after increasing fall height and larger-radius solutions after decreasing it.These regimes are illustrated in figure 11.
  • Validation: Surface tension validation shows good agreement with a reference curve that includes γ = 10−3, while surface tension markedly changes the coiling radius.The comparison uses Π1 = 7000, Π2 = 7, and Π3 = 10.3 10−3.
  • The viscous sewing machine: The moving-belt simulation reproduces successive translated coiling, alternated loops, and meanders as belt velocity increases, matching the reported experimental sequence.At still larger velocities, oscillations disappear and the pattern becomes straight; more complex patterns arise at larger fall heights.
  • Limitations and perspective: The linear implicit scheme updates velocity from a locally linear viscous-force expression before updating position, improving stability relative to an explicit scheme.It is harder to implement than an explicit scheme and dissipates angular momentum, unlike the explored nonlinear implicit approach.
  • Discussion: The discrete model combines stretching, bending, twisting, inertia, large rotations, and surface tension, and is derived from a Lagrangian formulation with discrete twist.Its specialized implementation can avoid implementing the Hessian matrix needed for naturally curved elastic rods, while future work includes more general constitutive laws.

Appendix A. Equivalence with the constitutive equations of Ribe

The appendix reformulates Ribe’s steady-coiling analysis in the paper’s Lagrangian, moving-frame formalism and shows that the resulting bending-moment expressions agree with Ribe’s constitutive laws.

  • Related analysis: Ribe’s steady helical-coiling solutions describe viscous jets falling onto a plane through nonlinear boundary-value equations solved by numerical continuation.The analysis uses the frame rotating with the jet, where the centerline shape is stationary, and AUTO software for continuation.
  • Equivalence: The paper establishes equivalence between Ribe’s three-dimensional Stokes-derived constitutive laws and its own equations (65).Because the two formalisms differ, Ribe’s analysis must be reworded in the paper’s framework.
  • Twist-curvature representation: Ribe’s Eulerian twist-curvature vector is related to the paper’s Lagrangian variant by using s rather than S as the differentiation parameter.Its components are decomposed in the moving material frame, whose orientation varies with arc length.
  • Twist-curvature representation: The Eulerian kinematical twist and binormal curvature are obtained by decomposing the twist-curvature vector into transverse and tangential components.This parallels the decomposition used earlier in the paper for the Lagrangian representation.
  • Rotating-frame kinematics: In the coiling frame, steady shape and axial velocity determine the material angular velocity, which transforms to the laboratory frame by adding the coiling rotation.The laboratory-frame relation is ω = ωR + Ωez, where Ω is the coiling frequency.
  • Constitutive-law comparison: Substituting the resulting rotation gradient into constitutive law (65b) yields bending moments identical to those derived by Ribe.The rotation-gradient expression includes a term arising because the material frame is moving, and is assembled using equations (A.2) and (A.4).
Loading 1202.4971v2…