# A numerical investigation of the fluid mechanical sewing machine

## Abstract

A thin thread of viscous fluid falling onto a moving belt generates a surprising variety of patterns depending on the belt speed, fall height, flow rate, and fluid properties. Here we simulate this experiment numerically using the Discrete Viscous Threads method that can predict the non-steady dynamics of thin viscous filaments, capturing the combined effects of inertia, stretching, bending and twisting. Our simulations successfully reproduce nine out of ten different patterns previously seen in the laboratory, and agree closely with the experimental phase diagram of Morris et al. (2008). We propose a new classification of the patterns based on the Fourier spectra of the longitudinal and transverse motion of the point of contact of the thread with the belt. These frequencies appear to be locked in most cases to simple ratios of the frequency of steady coiling obtained in the limit of zero belt speed. In particular the intriguing ‘alternating loops’ pattern is produced by combining the first five multiples of .

###### Contents:

## I Instabilities of viscous threads

### i.1 Introduction

A thin stream or jet of liquid falling onto a fixed surface is one of
the simplest situations in fluid mechanics, yet it can generate a
remarkable range of phenomena. Fast jets produce hydraulic jumps,
which are circular when the viscosity is very low Rayleigh (1914); Watson (1964) and
polygonal when it is somewhat higher Ellegaard *et al.* (1998). Thin streams
of very viscous fluid can exhibit steady coiling Barnes and Woodcock (1958) or
rotatory folding Habibi *et al.* (2010), and under some conditions coiling
produces spiral waves of air bubbles in the thin fluid layer spreading
over the surface Habibi *et al.* (2008). Finally, thin streams of
non-Newtonian fluid can exhibit the Kaye (”leaping shampoo”) effect
in which the stream rebounds episodically from the pile of previously
deposited fluid Kaye (1963).

A further degree of complexity is introduced if the source of the jet
and the surface onto which it falls are in relative motion. This is
the case when a home cook lays down ‘squiggles’ of icing or molten
chocolate on a cake, or when an artist lets paint dribble onto a
canvas from a moving paintbrush, a technique used to great effect by
Jackson Pollock Herczynski *et al.* (2011). An analogous situation involving
many interacting jets is the ‘spunbonding’ process of non-woven fabric
production, in which thousands of threads of molten polymer solidify
and become entangled as they fall onto a moving belt, producing a
fabric with a random texture.

Here we use a numerical approach to study an idealized model for these processes: the continuous fall of a single thread of viscous fluid onto a belt moving with a constant velocity in its own plane (Figure 1).

This system was first studied experimentally by Chiu-Webster and Lister Chiu-Webster and Lister (2006) (henceforth CWL), who called it the ‘fluid mechanical sewing machine’ on account of the stitch-like patterns traced on the belt by the thread. The complexity and diversity of these patterns testifies to a rich nonlinear dynamics and bifurcation structure. The appeal of the system is further increased by the theoretical and numerical challenges involved in modeling it.

In CWL’s experiments, viscous fluid (corn syrup) with density , surface tension coefficient and viscosity was ejected at a volumetric rate from a vertical nozzle of diameter , from which it fell a distance onto a belt moving at speed (fig. 1). The experiments were conducted by varying and for several different combinations of values of , and . When greatly exceeded a fall height-dependent critical value , the fluid thread had the form of a steady dragged catenary. As the belt speed was decreased towards , the lowermost part of the thread evolved into a backward-facing ‘heel’, which became unstable to periodic meandering when . Further decrease of the belt speed led to a series of bifurcations to more complex patterns (fig. 5), ending with the establishment of steady coiling for . CWL successfully predicted the shape of the steady dragged catenary using a ‘viscous string’ model that neglected bending stresses. However, this solution ceases to exist when the extensional axial stress at the bottom of the thread becomes zero, corresponding to the incipient formation of a heel in which the axial stress is compressional. Because a state of axial compressive stress is a necessary condition for the buckling of a slender body Taylor (1968), CWL interpreted the onset of meandering as a buckling instability of the heel.

Ribe et al. Ribe *et al.* (2006) carried out a numerical linear
stability analysis of the dragged catenary state to predict the
critical belt speed and the frequency for the onset
of meandering, using a more complete ‘viscous rod’ theory
incorporating bending and twisting of the filament. The numerical
predictions of and thereby obtained agree closely
with the experimental measurements of Chiu-Webster and Lister (2006).
Ribe et al. Ribe *et al.* (2006) also documented a close
correspondence between incipient meandering and finite-amplitude
steady coiling on a motionless () surface, such that
is nearly identical to the steady coiling frequency for any
given fall height . Moreover, the critical belt speed is
nearly identical to the vertical (free-fall) speed of the fluid
at the bottom of the thread, indicating that meandering sets in when
the belt is no longer moving fast enough to carry away in a straight
line the fluid falling onto it.

More extensive experiments were conducted by Morris et al.
Morris *et al.* (2008), using a carefully engineered apparatus with
silicone oil as the working fluid for better stability and
reproducibility. They determined a complete phase diagram for the
patterns as a function of and for a particular set of values
of the fluid viscosity and the flow rate . They showed that
the observed amplitude of meandering as a function of the belt speed
is consistent with a Hopf bifurcation, and proposed a simple model to
predict it based on the hypothesis that the contact point moves at
constant speed relative to the belt. Finally, they proposed a generic
set of amplitude equations which they used to characterize the
alternating loops (which they called ‘figure-of-eight’) and translated coiling patterns.

Most recently, Blount and Lister Blount and Lister (2011) performed a detailed asymptotic analysis of a slender dragged viscous thread, focussing on the structure of the heel. They showed that the lowermost part of the thread can exhibit three distinct dynamical regimes depending on whether the belt speed is greater than, nearly equal to, or less than the free-fall speed . Their asymptotic stability analysis of these steady states indicates that meandering sets in when the horizontal reaction force at the belt begins to be slightly against the direction of belt motion, corresponding to the heel ‘losing its balance’.

As the above summary indicates, our current theoretical understanding
of the fluid-mechanical sewing machine is essentially limited to the
initial bifurcation from the steady dragged configuration to
meandering. In this paper, we push forward into the fully nonlinear
regime with the help of a new computational
algorithm Bergou *et al.* (2010); Audoly *et al.* (2011)
that permits for the first time robust numerical modeling of arbitrary
non-stationary dynamics of viscous threads. After describing the
method briefly, we use it to generate a complete phase diagram of
sewing-machine patterns that reproduces all the major features of the
diagram determined experimentally by Morris et al.
Morris *et al.* (2008). We then perform a detailed Fourier analysis
of the motion of the contact point for each of the patterns simulated,
and propose a new classification of them based on the spectral content
of the motions of the contact point in two orthogonal directions.

For most of the patterns studied, we find that the frequencies present in the spectra of the contact point motion are multiples of the steady coiling frequency , indicating that the dynamics of the sewing machine are closely related to those of steady coiling. Accordingly, we set the stage for our study with a brief summary of the latter in the next section.

### i.2 Steady coiling

In steady coiling, the contact point of the thread with the surface moves with a constant angular velocity along a circle of radius (Figure 2a).

In most cases, the thread comprises two distinct parts: a long, roughly vertical ‘tail’ which deforms primarily by stretching under gravity, and a helical ‘coil’ in which the deformation is dominated by bending and (to a lesser extent) twisting. Thus the radius of the thread within the coil is nearly constant. By conservation of mass, the axial speed of the fluid in the coil is .

Steady coiling can occur in several distinct dynamical regimes
characterized by different balances of the viscous, gravitational, and
inertial forces acting on the thread Mahadevan *et al.* (1998); Ribe (2004); Ribe *et al.* (2006). These regimes appear clearly on plots of
the coiling frequency vs. the fall height. For convenience, we
define a dimensionless fall height and a dimensionless
coiling frequency :

(1) |

Figure 2b shows as a function of for the
parameters of the experiments of Morris et al.
Morris *et al.* (2008), viz., ,
,
, and
. For , coiling occurs in
a gravitational (G) regime. Inertia is negligible everywhere in the
thread, which is governed by a balance between gravity and the viscous
forces that resist stretching (in the tail) and bending (in the coil).
At intermediate heights , a complex
inertio-gravitational (IG) regime appears, in which the coiling
frequency is a multivalued function of the fall height. The
centrifugal force now becomes important in the tail, which behaves as
a distributed pendulum with an infinity of eigenmodes whose
corresponding eigenfrequencies are proportional to the simple pendulum
frequency . The first three of these frequencies are
shown by the dotted lines in Figure 2b. When one of these
eigenfrequencies is close to the frequency set by the coil, the tail
enters into resonance with the latter, giving rise to resonance peaks
that appear as rightward-facing ‘bumps’ in the curve . For
large heights , coiling occurs in an inertial (I) regime
in which the viscous bending force in the coil is balanced by
inertia Mahadevan *et al.* (1998).
The tail in this regime is almost perfectly vertical, and is
controlled by a balance between gravity, the viscous stretching force,
and the axial momentum flux. Finally, there is also a viscous (V)
regime in which both gravity and inertia are negligible everywhere in
the thread, but this is only observed when both and are much
smaller than in the experiments of CWL and Morris et al.

In a typical laboratory experiments on steady coiling, the parameters and and the fluid properties , and are held fixed while is varied. Each such experiment is completely defined by the values of the three dimensionless groups

(2) |

As an example, , , and for
all the experiments of Morris et al. Morris *et al.* (2008). The
(dimensionless) functional dependence of the coiling frequency on the
other parameters now takes the general form

(3) |

The effect of surface tension is measured by the parameter . Surface tension modifies the coiling frequency quantitatively but introduces no essentially new dynamics, as can be seen by comparing the solid () and dashed () curves in Fig. 2b.

The continuation method used to generate the curves in Fig. 2b can be used for steady coiling because the flow is stationary in a co-rotating reference frame that moves with the contact point. No such reference frame exists for the sewing machine configuration. We therefore require a different numerical method, which is described in the next section.

## Ii Numerical method

Our numerical experiments of the sewing machine were set up using the
computational method of Discrete Viscous Threads (henceforth DVT) originally
described in a conference
paper Bergou *et al.* (2010), and
presented in detail in a recent
paper Audoly *et al.* (2011).
To the best of our knowledge, DVT is the
only available numerical method that is fast and robust enough to be
applicable to the sewing machine geometry while retaining all the
relevant modes of deformation, namely stretching, twisting and
bending. For a thin thread the stretching modulus varies as the square of the thread’s radius, while the bending and twisting moduli are proportional to the fourth power of the radius. As a result the dynamics of
thin threads is a nonlinear and numerically stiff problem. The DVT
method addresses this difficulty by introducing an elaborate spatial
discretization of the equations based on ideas from differential
geometry. The method allows simulations to be run for arbitrary mesh
sizes, even very coarse ones, with optimal stability. This contrasts
with conventional discretization schemes which are typically stable
for sufficiently small mesh sizes only, the maximum mesh size being in
practice a small fraction of the smallest length scale in the flow,
here the small size of the coiled region at the bottom.

Below we briefly present the principles of the DVT method, introduce adaptive mesh refinement which provides a tremendous speed-up of the calculations when gravity stretches the tail significantly, validate the code against known solutions for steady coiling, explain the details of the numerical procedure, and finally present our numerical results.

### ii.1 Smooth case: the equations for thin viscous threads

The numerical code makes use of the Lagrangian approach and the viscous thread is discretized as a polygonal line. A mass is assigned to each vertex, forces are set up on these masses, and the motion of each mass is obtained by time integration of the fundamental law of dynamics. The discrete forces are designed in such a way that the motion of the polygonal line is equivalent to that of a thin viscous thread in the smooth limit. The key element of the numerical method is the expressions for these discrete viscous forces. To prepare the derivation, we start by reformulating the smooth case, usually expressed in Eulerian variables, in the Lagrangian framework. We introduce a Lagrangian coordinate that marks cross-sections and follows them during motion. This Lagrangian coordinate is defined as the arc-length in an imaginary reference configuration where the thread is a cylindrical tube of constant radius. It plays the same role as the vertex index in the discrete case.

At a particular time , the configuration of the thread is defined by its center-line and an orthonormal triad . This triad allows one to keep track of twisting, because the directions of and follow the orientation of a cross-section as it spins about the center-line. Thin threads deform is such a way that cross-sections remain approximately planar and perpendicular to the centerline. As a result, the center-line tangent and the unit vector are aligned. Denoting by the axial stretch factor based on the reference configuration, given by the norm of , we have:

(4) |

Since the triad is orthonormal, its time evolution defines a rigid-body rotation for any particular value of . The associated instantaneous angular velocity is such that

(5) |

One can take advantage of the fact that the vectors and
must remain aligned by equation (4)
to capture the twisting motion of the thread using a single parameter,
the angular spinning velocity defined by . In
this view, which we call the centerline/spin
representation Audoly *et al.* (2011),
the material velocity is a secondary variable which can
be reconstructed by the equation

(6) |

where the time derivative is denoted using a dot. In a viscous thread, the fundamental kinematical quantities are the strain rates, defined by

(7) |

where is the center-line velocity. Here captures the strain rate associated with axial stretching, while the vector captures in a combined manner the strain rates for the twisting and bending modes. The strain rate for stretching, , is related to the time evolution of center-line stretch by the formula .

For a thin thread, the internal stress is described by the resultant of the viscous forces over a particular cross-section, and their moment with respect to the origin of the cross-section. These internal force and moment vectors play the same role as the tensor in 3D continuum mechanics. For the special case of a thin thread geometry, Stokes’s constitutive law states that stress is proportional to the rate of deformation: for the stretching mode, we have

(8) |

and for the bending and twisting modes,

(9) |

Here is the dynamic viscosity, the cross-sectional area, is the moment of inertia, is the radius, is the tangent projection operator, and is the normal projection operator. The radius is a dependent variable which is reconstructed by the incompressibility condition

(10) |

where is the radius in the configuration of reference (note
that in the reference configuration by definition). The
stretching modulus was originally derived by
Trouton Trouton (1906),
and the bending modulus can be found for instance
in Buckmaster *et al.* (1975).

These equations are complemented by the balance of linear and angular momentum, known as the Kirchhoff equations for thin rods:

(11) | ||||

(12) |

Following an approximation introduced by Kirchhoff himself which is valid for thin threads, we have neglected the rotational inertia in the second equation. The vector is the resultant of the external forces per unit reference length . The weight of the thread and the surface tension are taken into account by setting

(13) |

where is the cross-sectional area in the reference configuration, the acceleration of gravity, and is the surface tension at the fluid-air interface. The last term in equation (13) is the net force on the center-line due to surface tension at the fluid-air interface, the argument in the derivative being the axial force due to the capillary overpressure . Note that there is no need to consider a linear density of applied torque in equation (12) for the problem at hand.

### ii.2 A variational view: Rayleigh potentials

The equations of
motion (11–12)
and the constitutive
law (8–9)
can be discretized in a natural
manner Audoly *et al.* (2011)
if they are first re-written in terms of a Rayleigh potential. The
Rayleigh potential yields the power dissipated by
viscous forces as a function of the vertex velocity and spinning velocity . For a thin thread, it has three contributions
corresponding to the stretching, bending and twisting modes,
. Note that the Rayleigh
potential also depends on the current configuration
but this dependence will be implicit in our
notations. As an illustration, consider the stretching contribution
. It only depends on the vertex velocities,
and reads

(14) |

where in the integrand is given by
equation (7), and is the
infinitesimal arc-length in current configuration. The expressions
for the twist and bending contributions can be found in
reference Audoly *et al.* (2011).

The Rayleigh dissipation potential is useful as it captures the effect of the internal viscous stress in a compact mathematical form. Indeed in the equations of motion (11–12) the net viscous force and the net viscous moment in the left-hand sides can be shown to be equivalent to a density of applied force

(15) |

and a density of applied twisting torque

(16) |

In these equations, the right-hand sides denote functional derivatives, as the dissipation potential takes the functions and as arguments. The hat notation expresses the fact that the derivatives have to be calculated formally first, and evaluated with using the real motion next.

### ii.3 Discretization

In the discrete case, the center-line of the thread is represented by a polygonal chain of particles . The length and unit tangent of an edge are defined by

(17) |

We consider viscous threads having a circular cross-section. As a result there is no need to keep track of the absolute orientation of the cross-sections during motion. Twist is taken into account through the instantaneous angular velocity of spin of an edge, noted . This is an unknown of the motion, for which we will derive an equation. Representing rotations by a single degree of freedom is beneficial for the simulation. The angular velocity vector is a secondary quantity in the simulation, which is reconstructed from the spinning velocity by an equation similar to the smooth equation (6).

The generalized velocity of a viscous thread is obtained by complementing the vertex velocities with the spinning velocities of the edges . The dynamics of the thread is specified by a differential equation involving the position , velocities , , as well as the acceleration . Rotational inertia is neglected and so does not enter into the equation. This differential equation is derived next.

As in the smooth case, internal viscous stress is captured by means of a discrete Rayleigh dissipation potential which is the sum of three contributions, . As an illustration, the stretching contribution is defined in close analogy with equation (14) by

(18) |

where is a discrete stretching modulus
defined by analogy with equation (8),
and is a discrete axial strain rate defined by analogy
with equation (7). For an expression of the
twist and bending contributions and
, we refer the reader
to Audoly *et al.* (2011).
They make use of discrete notions of curvature and twist, based on
ideas from discrete differential geometry.

In analogy with the smooth case, we write the equations of motion as

(19) | ||||

(20) |

The first equation is a balance of linear momentum for the vertices
and is associated with the positional degrees of freedom
, while the second equation is a balance of twisting
torque at each edge, associated with the spinning degrees of freedom
. Here, is the diagonal matrix filled
with the mass of the vertices, and
are the viscous forces and twisting
torques representing the internal stress in the thread, and
combines the weight and the net force on vertices due to
surface tension. For a detailed derivation of the discrete surface
tension forces,
see Audoly *et al.* (2011).
As in equation (12) for the smooth
case, we have neglected rotational inertia in the right-hand side of
equation (20): we have checked that
this has negligible effect on the simulation when the thread is thin
enough.

Our discretization of the thread is based on the Rayleigh potentials, and the discrete viscous forces and moments are defined as in equations (15) and (16) by:

(21) | ||||

(22) |

The equations of the present section provide a complete system of equations for the dynamics of a discrete viscous thread. For the time discretization, we use a semi-implicit Euler scheme, which provides good stability even for quite large time-steps (by semi-implicit, we mean that we linearize the equations near the current configuration at every time-step, before applying an implicit scheme). The treatment of boundary conditions is explained in section II.5.

### ii.4 Adaptive mesh refinement

The DVT method uses a Lagrangian grid which is advected by the flow. In sewing machine experiments, gravity can typically stretch the centerline by a factor 5 to 10 over the course of the descent. In the absence of refinement, this makes the grid very inhomogeneous: to capture the small-scale features near the belt one has to use an extremely fine mesh size near the nozzle. Overall, a large number of degrees of freedom are required and the simulation is inefficient. In addition, important inhomogeneities in edge lengths makes the time-stepping problem ill-conditioned and robustness is affected. To overcome these difficulties, we have set up adaptive mesh refinement.

In our implementation of refinement, edges are subdivided whenever their length exceeds a nominal length, which is a prescribed function of the distance to the belt. In the upper part of the belt, this nominal length is equal to twice the initial segment length, defined by the periodic release of new (Lagrangian) vertices from the nozzle. To better resolve the coil region, this nominal segment length decreased near the belt according to a prescribed exponential profile. This profile was adjusted in such a way that the coil region always includes a sufficient number of vertices, typically 10 to 30, with a final edge length typically below , and that the interval between subsequent subdivisions of a given edge is always larger than two simulation steps.

Whenever an edge was marked as needing subdivision, a new vertex was inserted. We computed the properties of the new vertex and of the two new edges as follows: the Lagrangian coordinate of the new vertex is obtained by linear interpolation, the mass stored in the former edge is equally split among its children, the position and velocity of the new vertex are calculated by an interpolation of order 4, the cross-sectional area , the spinning velocities, and the viscosities of the new edges are obtained by second-order interpolation. These interpolation orders were chosen in such a way that the bending and twist forces remain continuous upon subdivision. We used the steady coiling geometry to adjust the subdivision parameters and to validate the subdivision scheme.

### ii.5 Emission from the nozzle, capture by the belt

We found that the implementation of boundary conditions was a key point to successfully reproduce the patterns and the phase diagram of the experimental sewing machine. We tried simple implementations first, and could reproduce the curves for the frequency of steady coiling as well as the simplest stitch patterns, but failed to reproduce entire regions of the phase diagram. Further examination revealed the presence of spurious oscillations in the calculated acceleration in the steady coiling geometry (), even though the coiling frequency was correctly predicted. We removed these spurious oscillations by a more careful account of the boundary conditions both at the nozzle and at the belt, as explained below. Suppressing these oscillations appeared to be sufficient to bring the numerical predictions in close agreement with the experimental ones.

New vertices need to be emitted periodically from the nozzle. In a first implementation of the clamped boundary conditions there, we considered an infinite string of fluid particles which were moved with the prescribed ejection velocity , until they passed below the nozzle and were released. The position of the first vertex clamped inside the nozzle varies abruptly as a new vertex is released, and this was the cause of unwanted oscillations. They were suppressed by considering a string of two particles in the nozzle, with a fixed position relative to the nozzle; the injection velocity is then modelled by steadily increasing the length of the second edge, and periodically inserting a new vertex in third position.

Impact with the belt was first handled by detecting penetration of
vertices into the belt at the end of every time-step, and constraining
their velocity to match the belt’s velocity at any subsequent time.
This also induces large unwanted fluctuations in the acceleration,
which can be interpreted by the fact that the vertical momentum
resulting from the collision is not transferred to the thread until
the following time step. The oscillations were removed by using a
technique known as roll-back. Then, every time-step involves an
iteration whose aim is to determine the set of vertices undergoing a
collision during the time-step: whenever an unexpected collision takes
place, the step is discarded, time is rolled back, and a new time-step
is computed, forcing additional vertices to land on the belt at the
end of the time-step. An additional difficulty in the implementation
of roll-back in the context of a linearized implicit scheme, is that a
only small number of particles can be captured at every step. We
circumvented this difficulty by using adaptive time refinement. Such
an adaptive time refinement is presumably not needed if a fully
(non-linear) implicit scheme is used, such as the one presented in
Ref. Bergou *et al.* (2010).

### ii.6 Validation

The numerical code was validated by comparing its predictions of the steady coiling frequency to the predictions of the continuation method of Ribe Ribe (2004). The agreement is excellent for all fall heights (Figure 3). The hysteresis of the dynamic simulation in the range is physical, but inaccessible to the continuation method because the latter records all steady-state solutions regardless of their stability.

## Iii Simulation results

Our dynamic simulations of the sewing machine patterns, based on the
DVT method were carried out with non-dimensional quantities. This is
achieved by setting the density , the viscosity , and the
gravity to the value 1. This choice implies that both the length
scale and the time scale of the
problem introduced in equation (1) are equal to 1. The three
other physical parameters namely the nozzle diameter, the flow rate
and the surface tension were chosen to match the values of
, and prescribed by Morris et
al. experiments Morris *et al.* (2008): , and . The simulations were initiated
from a vertical thread of uniform radius comprising 172 segments of
equal length, falling from a height (gravitational
regime) onto a belt at rest. To avoid dealing with a shock when the
thread hits the belt, the simulation was started with the bottom end
of the thread clamped into the ground. The simulation was run until
the radius settled to a steady profile as a function of the elevation,
and a steady state of coiling was established. Next the space was sampled by slowly varying or in
turn.

As an illustration of the capabilities of our numerical technique, we show in Figure 4 a simulation of the ’translated coiling’ pattern that occurs for relatively low belt speeds.

Fig. 4a shows a three-dimensional view of the falling thread and the trace it lays down on the belt. The simulation time is for one period of this pattern using a Intel Core i7 processor and of DDR3 memory.

Fig. 4b shows the trajectory of the contact point in the
frame of the nozzle (solid line). Interestingly, the belt motion
breaks the symmetry of steady coiling not only in the longitudinal
-direction, but also in the transverse -direction.
Fig. 4b also shows the velocity of the contact point
relative to the moving belt as a function of position along the
trajectory (thin lines and green arrows). The magnitude of the
relative velocity varies significantly over a period, in contrast to
the meandering pattern for which it is nearly uniform
Morris *et al.* (2008). The relative speed is maximum at point A
where the contact point is moving upstream along the belt, and very
small at C where the motion is downstream. Finally,
fig. 4c shows the amount of viscous power dissipated by
the various modes, as a function of arc-length along the thread. The
curves cross each other at a height that corresponds to
the transition from the bending-dominated coil, to the
stretching-dominated tail. Thanks to adaptivity, the coil is well
resolved and the curves for the viscous power dissipation remain
smooth there, even though they vary rapidly.

In addition to the dimensionless parameters , and in equation (2) that describe the fluid properties and the ejection conditions, the patterns depend on the dimensionless fall height in equation (1), and the dimensionless belt speed

(23) |

Our simulations were carried out by varying and for
fixed values of and corresponding to the
experiments of Morris et al.Morris *et al.* (2008). Some of our
simulations were done with a surface tension parameter matching the
experimental value ; for reasons of numerical stability, however, most of the simulations used .

Figure 5 summarizes all the types of patterns that were encountered when scanning the () plane in the simulations, together with their experimental equivalents Chiu-Webster and Lister (2006).

The simulation reproduces nine out of the ten patterns reported by
Morris et al.Morris *et al.* (2008) and observed by CWL
Chiu-Webster and Lister (2006). The only missing pattern, the slanted
loops, will be discussed later on. To ease the understanding, the
name ’figure-of-8’ is replaced by ’alternating loops’ which more
accurately describes the pattern in Figure 5C.

In Figure 6 a phase diagram lists all the numerical patterns encountered in the space for .

The effect of surface tension is included, in
the simulation. For comparison, the patterns observed experimentally
by Morris et al.Morris *et al.* (2008) are shown by dots. The
agreement between the simulations and the experiments is remarkable:
all four patterns that were observed experimentally in this region of
the parameter space, namely translated coiling, alternating loops,
meanders and catenary, are recovered. In addition the location of the
boundaries separating the different patterns are reproduced
accurately.

The inclusion of surface tension gives rise to numerical instabilities
for heights approximately above , which we have not
been able to overcome by decreasing the mesh size or the time-step.
This is why there is no simulation data shown in the lower right
portion of Figure 6. Since surface tension is not
expected to modify qualitatively the dynamics of threads (see
Fig. 2) we investigated the case of larger fall heights
without any surface tension (). Five new patterns were
observed for larger fall heights, as shown in Figure
7, namely double coiling, double meanders, stretched
coiling, W-pattern, disordered pattern, despite the absence of surface
tension in the simulations. The new portion of the phase diagram,
, is very similar to that reported by Morris et
al.Morris *et al.* (2008), shown in the inset in
Fig. 7. In both diagrams, the alternating loops
pattern disappears at a critical height, beyond which there is a
substantial height ‘window’ containing only simple patterns (catenary,
translated coiling, and meanders). When the height is increased,
three patterns having a complex aspect (disordered pattern, stretched
coiling and the double meanders) all appear together at nearly the
same height. Finally, for some values of the height, disordered
patterns, shown in grey in the diagram, occur in two separate ranges of
the belt speed, with stretched coiling in between.

It is instructive to compare the numerical and experimental phase diagrams with the curves of frequency vs. height for steady coiling, calculated with same value of surface tension ( for the simulations, for the experiments). These curves are shown below each phase diagram in Figure 7. The comparison reveals that some of the more complicated patterns (stretched coiling, W-pattern, disordered pattern) appear at heights close to that for the onset of the multivalued IG regime in the steady coiling. In the steady coiling geometry, it is known that the multivalued regime is caused by the competition of several ‘viscous pendulum’ modes. This suggests that the complex patterns of the sewing machine are produced by the non-linear interaction of these pendulum modes.

Despite their similarities, the numerical (N) and experimental (E)
phase diagrams in Figure 7 exhibit some systematic
differences. In E, double coiling (pink) first appears at the same
height as the disordered (grey) and stretched coiling (yellow),
whereas in N it appears at significantly greater heights. Double
meanders (purple) have a common boundary with the catenary pattern
(black) in E, but occur only for significantly lower belt speeds in N.
In N, the catenary state can transition to disorder (grey) over a
range of heights, unlike in E. The range of belt speeds for double
coiling is significantly wider in N than in E. Finally, in N the
W-pattern is observed sporadically and for greater heights than in the
diagram in Fig. 7. Some of these differences are due
to the absence of surface tension in the simulations, and to the fact
that collisions of the viscous thread with its trace are not
prevented. Another explanation for the discrepancies may be the fact
that Morris et al. Morris *et al.* (2008) performed their pattern
recognition visually, whereas we used a more quantitative automatic
procedure based on Fourier decomposition. This is described in the
following section, where we review each of the individual patterns in
detail and propose a systematic classification of them.

## Iv Analysis of the patterns

To illustrate our pattern analysis, we fix the belt speed and vary the fall height, thereby following the horizontal line AB in the phase diagram of Fig. 7. In order of increasing height, the patterns seen along this line are meanders, alternating loops, translated coiling, and double coiling. We track the spectral content of these patterns continuously as they change smoothly or bifurcate. To do so, we focus on the trajectory of the contact point of the thread with the belt. Let and be its longitudinal and transverse coordinates in the laboratory (nozzle) reference frame, and define . Let and be the coordinates (also relative to the nozzle) at time of a material point that was laid down on the belt at time , and let be a generic point in the trace. The advection by the belt is expressed by:

(24) |

This equation means that the pattern seen on the belt is obtained by unfolding the motion of the contact point motion , as illustrated in Fig. 9a.

The numerical traversal of the line AB in Figure 7 required about 78500 time steps of size . We performed a Fourier analysis of the motion over a sliding window of 2000 time steps (the FFT spectrum was computed every 500 steps). The spectra obtained in this way typically comprise several well-defined peaks whose frequencies can be determined precisely (Fig. 9b).

Let , , etc. be the peak frequencies of the motion in the -direction, and , , etc. be those for the motion in the -direction. Because the fall height is slowly changing with time during the simulation, each observed frequency or can be plotted on a diagram to provide a ‘portrait’ of the evolving frequency content of the contact point motion, as a function of the fall height . The result is shown in Figure 8. The principal observed frequencies and are indicated in red and blue, respectively. Also shown for reference is the steady coiling frequency for the same fluid properties and ejection parameters (solid line), together with several multiples (1/3, 2/3, 1/2, 4/3, 5/3, 2) of that frequency (dashed lines).

The first pattern () is meandering, which is
characterized by two frequencies with a ratio
. At the lowest height where meandering first appears, the meandering frequency
is very close to the steady coiling frequency
for the same height, as expected on theoretical groundsRibe *et al.* (2006); Blount and Lister (2011). The next pattern () is the alternating loops, which has a rich spectrum
involving five multiples of . Translated coiling appears
next (), and is characterized by a single
frequency very close to the steady
coiling frequency. Finally, double coiling () has two frequencies and .

Figure 8 shows that the stitch patterns are combinations of motions in two orthogonal directions with frequencies closely related to the steady coiling frequency . Accordingly, we now change our perspective and classify the patterns based on their frequency content rather than on the shape they lay down on the belt. This frequency analysis is used to set up an efficient tool for identifying the patterns and assembling the numerical phase diagram automatically. In addition it leads to a simple kinematic model that provides a unified description of the different patterns.

### iv.1 Spectral signature of a pattern

We illustrate our method using the example of meandering which in most cases is the first pattern to bifurcate from the catenary state as the belt speed decreases. The red and blue lines in figure 9b show typical spectra of the motion of the contact point in the longitudinal and transverse -directions, respectively. The amplitude of the transverse motion is much greater than that of the longitudinal motion, and the frequency of the latter is exactly twice that of the former, . This suggests that a meander pattern can be synthesized by retaining only the two main frequencies, viz.

(25) |

where and are the amplitudes of the longitudinal
and transverse motions, respectively. Here the phase difference is
required to reproduce the symmetry of the pattern. A similar
two-frequency model was used by Morris et al. Morris *et al.* (2008)
to analyze weakly non-linear meanders.
Figure 9c1 shows the contact line trajectory in the
frame of the nozzle implied by (25) with
, and Figure 9c2 shows the
corresponding meander pattern obtained by advecting the motion
(25) in the -direction with and
.
Based on Figure 9, we define as a ‘meander’ any pattern
whose longitudinal motion, compared to the transverse motion, has
twice the frequency, a much smaller amplitude, and a phase shift of
.

Generalizing the above example, we show that all the sewing machine patterns can be represented by a superposition of a few harmonic motions of the form:

(26) |

Two terms in each direction are sufficient to reproduce the main features of the patterns, and . In Eqn. (26), and are the amplitudes of the components of the motion with frequencies and , and and are the phases relative to the highest-frequency mode. We now show how each of the sewing machine patterns can be characterized in terms of the parameters that appear in Eqn. (26).

### iv.2 Identification of the patterns

The identity of each pattern is determined not by the absolute values
of the parameters in Eqn. (26), but rather by the
dimensionless groups that can be formed from them, namely frequency
ratios, amplitude ratios, and the relative phases . In the
following we identify the characteristic values of these groups for
each of the patterns in turn. For ease of reference, the results are
summarized in Table 1.

1. Translated coiling. This pattern occurs for and low belt speeds
(Figure 7). Figure 10A shows a
simulation of this pattern (upper left) and the corresponding Fourier
spectra of the longitudinal and transverse motions of the contact
point (upper right), and the steady coiling frequency .
The longitudinal and transverse motions have similar amplitudes and
are characterized by a single dominant frequency that is very close to the steady coiling frequency
for the same fall height . The peak
frequency deviates from its original value as the belt speed
increases. The amplitudes in both directions remain equal and an
almost circular shape is created. The panels at lower left show the
reconstructed motion in the frame of the nozzle (right) and on the
belt (left), calculated using (26) with ,
and .
Note that the experimental pattern shifts upwards or downwards as the
belt speed is increased; this shift does not affect the spectral
context, but could be taken into account by including a purely
imaginary constant in Eqn. (26).

Patterns | V | ||||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|

Translated Coiling | - | - | 0 | - | 0 | - | 1. | - | 1. | - | .5 | ||

Meanders | 2 | - | - | 0 | - | - | .2 | - | 1. | - | 1.4 | ||

Alternating loops | 1 | - | 1/2 | - | 0 | 0 | 1 | - | .5 | .5 | 0.33 | ||

Double Coiling | 1/2 | 1/2 | 0 | 0 | .5 | 1.5 | .1 | 1.5 | .5 | ||||

Double Meanders | 1/2 | - | - | - | 0 | - | 1. | 0 | 0 | 1.5 | .75 | ||

Stretched Coiling | 1 | 1 | 0 | 0 | 1. | .1 | .5 | .1 | .6 | ||||

W-pattern | 1 | 1 | 0 | 0 | 1. | .2 | .2 | .5 | .7 | ||||

Catenary | - | - | 0 | - | - | - | 0 | - | - | - | 0 | - | 1 |

Disorder | - | - | - | - | - | - | - | - | - | - | - | - | - |

2. Meanders. On Figure 7 this pattern is
seen for and a range of
intermediate belt speeds. The typical Fourier spectra of the
meandering pattern were previously shown in Figure 9.
The pattern is a superposition of a single longitudinal, and a single
transverse harmonic motions, with frequency ratio
, amplitude ratio
, and a relative phase (with
the convention ). Near the catenary/meander boundary
; farther from the boundary, the
meandering frequency departs significantly from . The
regular symmetrical meanders correspond to a phase difference of
between the two directions (fig. 9) However,
the pattern may deform into a bean-like shape in certain cases. This
situation was reproduced kinematically by reducing the phase difference
to a value close to .

3. Alternating loops. This pattern was called
‘figure-of-eight’ by CWLChiu-Webster and Lister (2006) and Morris et
al.Morris *et al.* (2008), but we prefer to call it ‘alternating
loops’. The domain of this pattern is an elongated ‘bubble’
sandwiched between translated coiling and meandering at relatively low
fall heights (Figure 7). This pattern
displays a remarkably rich frequency spectrum with five principal
peaks (Figure 10B-c). In contrast to meandering, the
motion with the largest amplitude is longitudinal, with a frequency
that locks onto the frequency (see
also Figure 8). The next largest peaks correspond to a
transverse motion with frequencies and
, both amplitudes being very close.
The harmonics in the longitudinal direction, and
in the transverse direction are also visible.

The frequencies of all five peaks can be written compactly as

Even though the spectra in Fig. 10B-b shows five peaks,
the kinematic model in Eqn. (26) generates an almost
identical pattern if one retains only the three main contributions
, , and with an
amplitude ratio , , and phases and
(Figure 10B-c). We used these characteristics of the
three main frequencies as a criterion for automatic detection of
alternating loops.

4., 5. Double coiling and double meanders. The Fourier spectra for these patterns are shown in Figures 10C and 11A, respectively. Both the longitudinal and transverse components have peaks at two frequencies and . The origin of these frequencies is clear from the uppermost part of Figure 8, which shows them as functions of fall height for double coiling. The range of fall heights in question is within the inertio-gravitational regime of steady coiling, for which the curve is multivalued in specific height ranges (Figure 2). The portion of that curve has two stable branches corresponding to different ‘pendulum’ modes of the tail: a lower branch (labelled LB in Figure 8) with , and an upper branch (UB) with 2.1-2.2. Figure 8 shows that the higher double coiling frequency stays locked to the upper branch of , which is the only stable one when . The lower frequency, by contrast, follows a ‘phantom’ branch with frequency that is very nearly a direct continuation of the lower branch of to greater fall heights. This behavior is possible because the ratio of the frequencies of the upper and lower branches happens to be quite close to 2.0.

Although double coiling and double meanders have the same frequency
content, they are distinguished by the relative amplitudes and phases of
the transverse and longitudinal motions. For double coiling, the amplitudes
of the two motions are the same at both frequencies
(, ), and
the relative phases are and .
For double meandering, by contrast, the transverse motion is
dominated by the frequency
while the longitudinal motion is dominated by
. The relative phases are
while .

6., 7. Stretched coiling and the W-pattern. These patterns
occur predominantly in the range of heights corresponding to
inertio-gravitational coiling (right-hand side of Figure
7). Their typical Fourier spectra are shown in Figure
11B and C, respectively. Like double coiling and double
meanders, their characteristic signature is the () frequency couple. But whereas double coiling and
double meanders are dominated by transverse motion at the frequency
, stretched coiling and the W-pattern are
dominated by longitudinal motion at the frequency .
The difference between stretched coiling and the W-pattern is only due to different relative amplitudes of the transverse and longitudinal motions denoted in Figure 11B-b and C-b. The direction is dominant in both cases but the peak at in the -direction is much closer for the W-pattern than for the stretched coiling. Theses differences cause the change of invagination between the two patterns (Figure 11B-a and C-a).
The phase difference between the longitudinal and transverse motions is in both cases.

8. Disorder. Disordered patterns appear in several parts of
the phase diagram (grey in Figure 7), primarily at
heights within the inertio-gravitational coiling regime. The typical
Fourier spectra of these patterns is very rich, with more than four
peaks in both the longitudinal and transverse directions with
comparable and strongly time-dependent amplitudes
(Fig. 12-i and corresponding FFT). Such
patterns are not transient between two steady patterns, as the
aperiodic behavior persists indefinitely in time.

9. Catenary. The catenary is obtained when the point of
contact is at rest in the nozzle frame, which happens in the upper
region of the phase diagram. The FFT spectrum is then empty.

In addition to the patterns discussed above, slanted
loops Chiu-Webster and Lister (2006) were also reported by Morris et
al. Morris *et al.* (2008) in a very narrow region of their phase
diagram. Slanted loops is a pattern wherein a buckle is periodically
laid down on the belt, and subsequently closes up into a loop when the
thread touches itself and coalesces. Our numerical simulations do not
account for self-contact of the thread, nor for surface
tension-mediated coalescence. This probably explains why we observed
slanted loops transiently only, as shown in Fig. 12-ii. A
similar argument may also explain why the W-pattern is observed for
significantly larger fall heights in the simulations than in the
experiments: considering self-contact of the thread would certainly
favor its existence over the stretched coiling pattern in the
simulation.

We observe that CWL Chiu-Webster and Lister (2006) reported yet another pattern, side kicks. Side kicks consist of small heaps of fluid regularly spaced along an otherwise perfectly straight trace. We suggest that this pattern is a limiting case of stretched coiling in which the amplitude of the transverse motion becomes very small relative to that of the longitudinal motion. Side kicks have not been reported in the experimental phase diagram of Morris and al., the one we attempt to reproduce in the simulations; consistently, we have observed this pattern in our simulation only transiently.

## V Discussion

The simulations of the fluid-mechanical sewing machine presented here were performed using a new numerical algorithm, Discrete Viscous Threads (DVT). The essential idea of DVT is to start from a complete geometrical and kinematical description of the thread in the discrete setting, and push the discrete approach as far as possible; in particular a discrete representation of viscous stress is built based on a variational view (Rayleigh potentials). This approach leads to a code that is robust even for quite large mesh size. This numerical method is used to simulate complex unsteady behavior of a viscous thread and offers good efficiency. As an example, the traversal (78500 time steps) of the line AB in Figure 7 required on a Intel Core processor. The accuracy of the method is demonstrated by its ability to reproduce the curve of steady coiling frequency vs. height, as predicted by an independent continuation method (Figure 2), and by the close agreement between the calculated and experimentally determined phase diagrams for the sewing machine patterns (Figure 6).

At each time step in an unsteady DVT simulation, any desired
kinematical or dynamical variable can be calculated as a function of
arc-length, i. e. at every vertex along the thread’s centerline.
Examination of these functions provides insights into the thread’s
dynamics. We saw an example in Figure 4c, which showed
the local rates of viscous dissipation of energy due to deformation by
stretching (or compression), bending, and twisting. The figure
reveals that the thread is divided into two distinct parts: a ‘tail’
in which the dissipation is dominated by stretching, and a ‘coil’ in
which it is dominated by bending and twisting, but with a significant
contribution from compression. The differential equations describing
bending are of higher order than those describing stretching, so
Figure 4c implies that the coil is an ‘inner’ solution or
boundary layer whose presence is required by the need to satisfy all
the relevant boundary conditions at the thread’s contact point with
the belt. While the boundary-layer character of the coil region has
long been recognized for steady coiling Mahadevan *et al.* (1998); Ribe (2004), our simulations open up the possibility of studying the
associated non-steady dynamics.

Another benefit of the simulations is to allow exploration of regions of parameters space that are inaccessible in the real world but provide new insights into the dynamics of the thread. In Figure 4c, the rates of viscous dissipation for the bending and twisting modes were added together. A real thread is made up of an incompressible fluid, and the ratio of the bending to the twisting modulus is always (both moduli are in particular proportional to the fourth power of the thread’s radius). In the simulation, it is possible to investigate the relative importance of bending and twisting in the thread’s dynamics by setting the twisting modulus to zero, taking in equation (9), while keeping the bending modulus unchanged. We performed additional DVT simulations for a perfectly twist-compliant sewing machine. The resulting phase diagram shows only minor differences with Figure 7, showing that twist plays a negligible role, relative to bending, in the selection of the stitch patterns.

The experiments and simulations of the viscous sewing machine reveal in the first instance a great diversity of patterns, whose relations to one another are not evident. Our goal was to characterize the patterns in a more unified way. This is possible by going beyond a visual identification and computing for each pattern the Fourier spectra of the longitudinal and transverse components of the motion of the thread’s contact point with the belt. We showed that each pattern has a distinct spectral signature consisting of isolated peaks at a small number of well-defined frequencies. The patterns differ from each other in the values of those frequencies, in their relative amplitudes, and in their distributions among the longitudinal and transverse modes.

A closer look shows that the frequencies in the spectra are closely related to the frequency of steady coiling of a thread falling on a motionless () surface from the same height. The precise nature of the relationship depends on the pattern considered. For meanders, the frequency of the transverse motion at onset, but then deviates significantly from as the belt speed is decreased beyond the critical value (Figure 8). In all the other patterns, however, the frequencies are locked to in some way. In stretched coiling, the dominant frequency of both the longitudinal and transverse motions is . Still more complicated are double meanders, double coiling, stretched coiling, and the W-pattern, for which the dominant frequencies are and . The presence of these frequencies reflects the nonlinear interaction of the two lowest modes of inertio-gravitational coiling, whose unforced frequencies differ by approximately a factor of 2 (Figure 8). Among periodic patterns, the richest spectral content is achieved by the alternating loops pattern, for which the five dominant frequencies are multiples of (Figure 8).

We proposed a simple kinematic model whereby each pattern is reconstructed by a superposition of a few frequencies, with appropriate amplitudes and relative phases in the longitudinal and transverse directions. This makes it possible to set up automated recognition of the patterns, and leads a classification of the patterns within a unified descriptive framework. The next step is to elucidate the physical mechanisms responsible for these simple spectral signatures. In future work we plan to investigate the similarities between the sewing machine and low-dimensional oscillator models with non-linear forcing, using the direct DVT simulations as a starting point to look into these analogies.

We are very grateful to Miklòs Bergou and Eitan Grinspun as the
collaborative development of the DVT method has made the present work
possible. Preliminary numerical results concerning the viscous sewing
machine
have been obtained in collaboration with them Bergou (2010); Bergou *et al.* (2010). We would like to
thank E. Grinspun for suggestions on the manuscript. During the
preparation of this manuscript, we learnt of an independent but
related experimental work by R. L. Welch, B. Szeto, and S. W. MorrisWelch *et al.* (2012).

### References

- L. Rayleigh, “On the theory of long waves and bores,” Proceedings of the Royal Society of London. Series A (1914).
- E. Watson, “The radial spread of a liquid jet over a horizontal plane,” Journal of Fluid Mechanics (1964).
- C. Ellegaard, A. E. Hansen, A. Haaning, K. Hansen, A. Marcussen, T. Bor, J. L. Hansen, and S. Watanabe, “Creating corners in kitchen sink flows,” Nature, 392, 767–768 (1998).
- G. Barnes and R. Woodcock, “Liquid rope-coil effect,” Am. J. Physics, 26, 205–209 (1958).
- M. Habibi, Y. Rahmani, D. Bonn, and N. M. Ribe, “Buckling of liquid columns,” Phys. Rev. Lett., 104, 074301 (2010).
- M. Habibi, P. C. F. Mller, N. M. Ribe, and D. Bonn, “Spontaneous generation of spiral waves by a hydrodynamic instability,” Europhys. Lett., 81, 38004 (2008).
- A. Kaye, “A bouncing liquid stream,” Nature, 197, 1001–1002 (1963).
- A. Herczynski, C. Cernuschi, and L. Mahadevan, “Painting with drops, jets, and sheets,” Physics Today, 31, 32–36 (2011).
- S. Chiu-Webster and J. Lister, “The fall of a viscous thread onto a moving surface: a ’fluid-mechanical sewing machine’,” Journal of Fluid Mechanics (2006).
- G. I. Taylor, “Instability of jets, threads and sheets of viscous fluid,” in Proceedings of the 12th International Congress of Applied Mechanics (1968).
- N. M. Ribe, J. R. Lister, and S. Chiu-Webster, “Stability of a dragged viscous thread: Onset of ”stitching” in a fluid-mechanical ”sewing machine”,” Physics of fluids, 18, 124105 (2006a).
- S. W. Morris, J. Dawes, N. Ribe, and J. Lister, “Meandering instability of a viscous thread,” Physical Review E (2008).
- M. J. Blount and J. R. Lister, “The asymptotic structure of a slender dragged viscous thread,” J. Fluid Mech., 674, 489–521 (2011).
- M. Bergou, B. Audoly, E. Vouga, M. Wardetzky, and E. Grinspun, “Discrete viscous threads,” Transactions on Graphics, 29, 116 (2010).
- B. Audoly, N. Clauvelin, M. Bergou, P.-T. Brun, E. Grinspun, and M. Wardetzky, “Simulating the dynamics of thin viscous threads,” (2011), Submitted to the Journal of Computational Physics.
- N. M. Ribe, “Coiling of viscous jets,” Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 460, 3223–3239 (2004).
- L. Mahadevan, W. Ryu, and A. Samuel, “Fluid rope trick investigated,” Nature, 392, 140 (1998).
- N. Ribe, H. Huppert, M. Hallworth, M. Habibi, and D. Bonn, “Multiple coexisting states of liquid rope coiling,” Journal of Fluid Mechanics, 555, 275–297 (2006b).
- F. R. S. Trouton, “On the coefficient of viscous traction and its relation to that of viscosity,” Proceedings of the Royal Society of London, A, 77, 426–440 (1906).
- J. D. Buckmaster, A. Nachman, and L. Ting, “The buckling and stretching of a viscida,” Journal of Fluid Mechanics, 69, 1–20 (1975).
- M. Bergou, Discrete Geometric Dynamics and Artistic Control of Curves and Surfaces, Ph.D. thesis, Columbia University (2010).
- R. L. Welch, B. Szeto, and S. W. Morris, “Frequency structure of the nonlinear instability of a dragged viscous thread,” arXiv:1201.0817v1 (2012).