# 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

- Rayleigh (1914) L. Rayleigh, “On the theory of long waves and bores,” Proceedings of the Royal Society of London. Series A (1914).
- Watson (1964) E. Watson, “The radial spread of a liquid jet over a horizontal plane,” Journal of Fluid Mechanics (1964).
- Ellegaard et al. (1998) 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).
- Barnes and Woodcock (1958) G. Barnes and R. Woodcock, “Liquid rope-coil effect,” Am. J. Physics, 26, 205–209 (1958).
- Habibi et al. (2010) M. Habibi, Y. Rahmani, D. Bonn, and N. M. Ribe, “Buckling of liquid columns,” Phys. Rev. Lett., 104, 074301 (2010).
- Habibi et al. (2008) 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).
- Kaye (1963) A. Kaye, “A bouncing liquid stream,” Nature, 197, 1001–1002 (1963).
- Herczynski et al. (2011) A. Herczynski, C. Cernuschi, and L. Mahadevan, “Painting with drops, jets, and sheets,” Physics Today, 31, 32–36 (2011).
- Chiu-Webster and Lister (2006) 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).
- Taylor (1968) G. I. Taylor, “Instability of jets, threads and sheets of viscous fluid,” in Proceedings of the 12th International Congress of Applied Mechanics (1968).
- Ribe et al. (2006) 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).
- Morris et al. (2008) S. W. Morris, J. Dawes, N. Ribe, and J. Lister, “Meandering instability of a viscous thread,” Physical Review E (2008).
- Blount and Lister (2011) M. J. Blount and J. R. Lister, “The asymptotic structure of a slender dragged viscous thread,” J. Fluid Mech., 674, 489–521 (2011).
- Bergou et al. (2010) M. Bergou, B. Audoly, E. Vouga, M. Wardetzky, and E. Grinspun, “Discrete viscous threads,” Transactions on Graphics, 29, 116 (2010).
- Audoly et al. (2011) 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.
- Ribe (2004) N. M. Ribe, “Coiling of viscous jets,” Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 460, 3223–3239 (2004).
- Mahadevan et al. (1998) L. Mahadevan, W. Ryu, and A. Samuel, “Fluid rope trick investigated,” Nature, 392, 140 (1998).
- Ribe et al. (2006) 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).
- Trouton (1906) 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).
- Buckmaster et al. (1975) J. D. Buckmaster, A. Nachman, and L. Ting, “The buckling and stretching of a viscida,” Journal of Fluid Mechanics, 69, 1–20 (1975).
- Bergou (2010) M. Bergou, Discrete Geometric Dynamics and Artistic Control of Curves and Surfaces, Ph.D. thesis, Columbia University (2010).
- Welch et al. (2012) R. L. Welch, B. Szeto, and S. W. Morris, “Frequency structure of the nonlinear instability of a dragged viscous thread,” arXiv:1201.0817v1 (2012).