Buckling, crumpling, and tumbling of semiflexible sheets in simple shear flow

Abstract

As 2D materials such as graphene, transition metal dichalcogenides, and 2D polymers become more prevalent, solution processing and colloidal-state properties are being exploited to create advanced and functional materials. However, our understanding of the fundamental behavior of 2D sheets and membranes in fluid flow is still lacking. In this work, we perform numerical simulations of athermal semiflexible sheets with hydrodynamic interactions in shear flow. For sheets initially oriented near the flow-vorticity plane, we find buckling instabilities of different mode numbers that vary with bending stiffness and can be understood with a quasi-static model of elasticity. For different initial orientations, chaotic tumbling trajectories are observed. Notably, we find that sheets fold or crumple before tumbling but do not stretch again upon applying greater shear.

---

1 Introduction

What do small 2D sheets do when subjected to fluid flow? Researchers have long studied problems of fluid–structure interactions ranging from flapping flags at macroscopic scales to tumbling polymer chains, such as DNA, at microscopic scales. This work follows a similar vein of research and, to our knowledge, represents one of the first studies of the dynamical behavior of athermal 2D sheets in shear flow at low Reynolds number. The quantification of different dynamical regimes explored here should help inform solution processing methods for nanomaterials and techniques to manipulate flexible materials via fluid flow. Such 2D materials of interest include graphene and graphene derivatives, clays, inorganic nanosheets, nacre-like materials, elastic membranes, colloidal membranes, 2D polymers, and 2D biological objects like kinetoplasts among others.

The dynamical behavior of semiflexible filaments and polymer chains dispersed in fluids is associated with a rich history of academic and industrial exploration. de Gennes famously studied the coil-stretch transition of polymers in extensional flow,¹ and subsequent theoretical² and experimental work^3–5 has looked at the behavior of polymers in shear flow. In particular, it has been shown that polymers stretch under shear but tumble somewhat erratically due to thermal fluctuations. Many additional simulation works^6–11 have also enhanced our understanding of the rheological properties of polymer chains. At perhaps a slightly larger length scale lies a significant body of work on fluid–structure interactions of flexible fibers and filaments.¹² Namely, fibers have been found to exhibit buckling^13–15 and morphological transitions¹⁶ as well as periodic and chaotic dynamical trajectories.^17–19

At high Reynolds number, one example of fluid-structure interactions with 2D objects that has attracted researchers' attention is the flapping of flags.^20,21 Relatively less work, though, has focused on the behavior of 2D semiflexible sheets in low-Reynolds-number flows, despite the increasing relevance of transition metal dichalogenides, graphene,²² graphene oxide,^23,24 and 2D polymers^25,26 in advanced technologies. Xu and Green^27,28 looked at sheets under shear and biaxial extensional flow, and Babu and Stark²⁹ have studied tethered sheets in fluids, confirming predicted scaling laws of Frey and Nelson³⁰ regarding thermal fluctuations in the presence of hydrodynamic interactions. Additionally, Dutta and Graham³¹ have studied Miura-patterned sheets and observed various interesting dynamical regimes involving periodic tumbling, unfolding, and quasiperiodic limit cycles. Most recently, Yu and Graham³² studied “compact-stretch” transitions of elastic sheets under extensional flow. However, there is still much to be learned about the dynamics of sheets subjected to fluid flow. Such fundamental understanding is particularly relevant now, as solution processing and liquid exfoliation^33–36 are commonly employed techniques in handling 2D materials, and nanosheet mechanical properties can be tuned by altering exfoliation protocols (e.g., solvent properties, the use of surfactant, etc.).^36,37 Furthermore, studies involving 2D materials whose functions depend significantly on coupled fluid motion are becoming more prevalent.^23,38–45

We consider asymptotically thin, isotropic, semiflexible sheets of size L immersed in a fluid at low Reynolds number. Sheets can be characterized by a bending rigidity, κ, and a 2D Young's modulus, Y, that have units of energy and force per length, respectively. The dimensionless ratio relating these two quantities is known as the Föppl–von Kármán (FvK) number and is equal to YL²/κ. From the classical theory of thin plates,⁵³ both κ and Y are also related to the 3D Young's modulus of the material, E, as Y = Eh and κ = Yh²/[12(1 − ν²)], where h is the thickness of the sheet and ν is the Poisson ratio. As reference, some examples of bending rigidities are provided in Table 1. For graphene, values can vary significantly from the theoretical microscopic bending rigidity (the first value) due to thermal fluctuation-induced stiffening, which is, in turn, experimentally challenging to measure and length-scale dependent.^54,55 It is also worth noting that many nanomaterials are produced containing multiple layers, and it has been shown that multilayer van der Waals materials exhibit intermediate behavior between classical thin plates and ideally lubricated stacks of sheets.⁵⁶ While these complications and others, such as nonzero slip,⁴³ are present for atomically thin nanomaterials, an athermal semiflexible membrane model may at least qualitatively describe the behavior of sheets with lateral dimensions much larger than the atomic scale and with (renormalized) bending rigidities much greater than the thermal energy (i.e., κ/k_BT ≫ 1).

Table 1 Approximate values of bending rigidities for some common materials and nanomaterials

Material	Bending rigidity (kBT|298K)	Literature references
Graphene	40, 280, 10^5	46, 47, 48
Graphene oxide	1	49
MoS2	370	50
Phospholipid bilayer	20	51
PMMA (100 nm film)	8 × 10^7	52

---

In this work, we quantify the behavior of high-FvK sheets subjected to shear flow using a numerical immersed boundary method. For different ratios of bending rigidity to shear energy (a dimensionless ratio we denote as S) and initial orientations, we find behavior ranging from buckling to transient tumbling and chaotic crumpling. Specifically, for initially flat sheets oriented near the flow-vorticity plane, quasi-1D buckling reminiscent of Euler buckling is observed, and a simple continuum elasticity model is presented to explain transitions in the buckling modes. The orientation of sheets that are stiff (relative to the shear strength) or initially oriented near the flow-vorticity plane is found to be well predicted by Jeffery's equations for thin oblate spheroids. However, deviations from the Jeffery orbits are observed for different initial orientations, with sheets of intermediate bending rigidity transiently tumbling before following a Jeffery-like trajectory and sheets of low bending rigidity crumpling into a compact structure and continuously tumbling in a chaotic manner. Summary statistics of sheet orientation (viz., the mean orientation and the orientational covariance matrix) are constructed and used to analyze the crumpling and chaotic motions quantitatively.

2 Model and methods

Hexagonal sheets of circumradius 58a were constructed by creating a surface triangulation with edges of length l = 2a and placing beads of radius a at each of the vertices (for a total of N = 2611 beads). These beads can be considered a geometric manifestation of the shortest resolvable length scale, a, in this coarse-grained sheet model. We caution readers to avoid thinking of a as a thickness since the sheets in this study behave hydrodynamically as though they are asymptotically thin, as discussed more below. We chose to use hexagonal sheets due to their symmetry and the density of the underlying triangular lattice of beads. Such a packing of beads is the closest possible in 2D, which is important for capturing self-avoidance and uniform hydrodynamic interactions, and represents a high degree of rotational symmetry, which makes the barrier to bending as a function of angle relatively uniform. These choices reflect a desire to create a discrete sheet model that tries to be as “isotropic as possible”.

Bending forces were modeled via dihedral forces over each pair of neighboring triangles Δ_i and Δ_j as:

U_bend(Δ_i,Δ_j) = κ(1 − n̂_i·n̂_j), (1)

where κ is the bending rigidity and n̂_i and n̂_j are (consistently oriented) triangle normals.^55,57–59 Note that such a model of bending may not be generally applicable for capturing bending forces of real materials (e.g., those that are nonlinear or plastically deform) or other models (e.g., the Helfrich-Canham model) at large curvatures. However, it is particularly desirable for efficiency and for capturing deviations from flat conformations. The equivalent continuum bending rigidity is given by .⁵⁹

Hard-sphere interactions between all of the beads (excluding those connected by a bond) were approximated via the pair potential,

where r is the distance between two interacting particles, and Δt is the integration timestep used. This pairwise potential displaces two overlapping particles to contact under Rotne–Prager–Yamakawa dynamics (discussed below) and accurately captures the thermodynamic properties of the hard sphere fluid.⁶⁰

Harmonic bonds of the form

were applied along each of the edges of the triangulation with stiffness k = 1000 × 6πη[small gamma, Greek, dot above]a and r₀ = 2a, where η and [small gamma, Greek, dot above] are the viscosity and shear rate of the surrounding fluid. The continuum 2D Young's modulus is related to k via the expression .⁵⁵ With these values of bond stiffness and bending rigidity, the Föppl–von Kármán (FvK) number of the system ranged in magnitude between 10⁴ and 10⁶, the latter of which is similar to that of a sticky note or a 100 nm-wide sheet of graphene. Thus, with these high FvK numbers and by scaling the in-plane stiffness with the flow strength, the sheets studied here can largely be considered inextensible, and we expect in-plane modes of motion to be generally irrelevant to the behavior observed.

For the surrounding fluid, we assume a low Reynolds number such that Stokes equations are valid. Hydrodynamic interactions between all of the beads in the sheet were included at the Rotne–Prager–Yamakawa (RPY) level.^61,62 That is, each bead produces a Stokes monopole and degenerate quadrupole but is “regularized” in such a way that ensures the mobility tensor remains positive definite for all (potentially overlapping) particle configurations. Although ostensibly a low-level of approximation for hydrodynamic interactions, it is important to note that the beads of the sheet act as regularized Stokeslets, so, in some sense, the use of large sheets with many beads employed here can be considered an approximation of the appropriate boundary integral for a smooth, continuous sheet.^63,64 Each 3 × 3 block of the RPY tensor mapping forces on particle j to the velocity of particle i is a function of the distance, r, between the particles and is given by:

where r̂ is a unit vector pointing from particle i to particle j. The particles were then advanced in time by integrating the following equation of motion:where is a stacked vector of N particle coordinates, U is the total potential energy and (L)_mn = [small gamma, Greek, dot above]δ_m1δ_n2 is the 3 × 3 velocity gradient tensor. Importantly, advecting the particles as “no-slip” point particles and not as rigid spheres (which would require more sophisticated integrators that operate at least at the stresslet level⁶⁵ and include translational-dipolar couplings) causes the sheet to behave as if it were asymptotically thin with respect to the shear flow in the sense that h/L → 0, where h is the thickness. This lack of thickness breaks the periodicity of Jeffery orbits and introduces stable flat states in the flow-vorticity plane, but otherwise does not significantly affect the trajectories away from these stable states compared to a very thin sheet. In fact, it is easy to show that the inverse of the Jeffery orbit period, T, of an oblate spheroid scales with the aspect ratio as for small values of h/L, so the periods of thin sheets can become arbitrarily long. Thus, we believe the phenomena reported here should be generalizable to more realistic sheets of small but nonzero aspect ratios.

Incorporating higher-order hydrodynamic interactions^65–67 poses several theoretical and computational challenges, the most important of which is the necessity to develop an elasticity model that properly constrains the relative rotation of beads within a sheet. Development of such a model in future work could be valuable.

Flat sheets initially constructed in the flow-vorticity plane were rotated by a varying angle ϕ about the flow axis and then by an angle θ = 5° about the vorticity axis (see Fig. 2). Rotating by θ perturbs the sheet away from the stable state (θ,ϕ) = (0°,0°). By symmetry, unique initial conditions are only generated by angles ϕ ∈ [0°,90°]. Simulations were conducted by integrating eqn (5)via a forward Euler scheme with a timestep of at most [small gamma, Greek, dot above]Δt = 5 × 10⁻⁴ using a custom plugin adapted from ref. 63 for the HOOMD-blue molecular simulation package⁶⁸ on graphics processing units (GPUs). This timestep was chosen to ensure numerical stability and was sufficiently small compared to the highest-frequency mode of motion in the system (i.e., the harmonic bonds between the beads). All simulations were conducted on NVIDIA GTX 1080 Tis and required thousands of GPU-hours.

3 Buckling

The dynamics of rigid ellipsoidal particles in simple shear flows was first worked out by Jeffery^67,69 with details for non-axisymmetric particles later studied by Hinch and Leal.⁷⁰ The differential equation (in terms of matrix–vector products) governing the orientation vector, p, of axisymmetric spheroidal particles iswhere Ω = (L − L^T)/2 is the antisymmetric vorticity tensor, E = (L + L^T)/2 is the symmetric strain-rate tensor, and a₁ and a₂. are the radii of the semi-axes that are parallel and orthogonal, respectively, to the axis of axisymmetry. In the limit of infinitely thin sheets, the rate-of-strain prefactor approaches −1, which is the value we used in comparing to our numerical simulations.

First, we considered sheets initially oriented with ϕ = 0°, which corresponds to a flat sheet whose normal is oriented in the flow-gradient plane. These sheets all began initially flat, buckled as they flipped in a way reminiscent of Euler buckling, and eventually flattened out again as they approached the stable flat state along the flow-vorticity plane. Note, again, that such a stable state only exists for infinitely flat sheets oriented in the flow-vorticity plane; sheets with finite thickness would eventually tumble and buckle again periodically, albeit with potentially long periods. Fig. 1 shows snapshots of a particular sheet buckling during its trajectory. The degree, or mode, of buckling depended on the bending stiffness and will be discussed more below.

Fig. 1 Snapshots of a buckling sheet with initial angle ϕ = 0° and dimensionless bending rigidity S = 4.61 × 10⁻⁴ at times [small gamma, Greek, dot above]t = 10, 10.75, 11.5, 12.25, and 13 from left to right. The axis labels 1, 2, and 3 represent the flow, gradient, and vorticity directions, respectively.

Fig. 2 Illustration of the initial angles, θ and ϕ, used for the initial conditions of the sheet as well as the length of the sheet, L, and the flow (1), gradient (2), and vorticity (3) directions of the imposed shear flow.

Fig. 3 shows the angle between the flow-vorticity plane and the line passing through the two ends of the sheets in the flow-gradient plane as a function of time and S, the dimensionless bending rigidity, compared to the angle predicted by the Jeffery orbit of an infinitely thin, rigid, oblate spheroid. This dimensionless bending rigidity—essentially the sheet's equivalent of the “elasto-viscous” number as it is sometimes called in the flexible filament literature¹⁶—is defined as

Here, L is the characteristic radius of the sheet (equal to 58), and the numerical factor of 1/π is related to the particular numerical discretization of the sheet (see Appendix). It is clear that the curves are all quite close to the Jeffery trajectory and become closer as bending stiffness increases (i.e., the sheet acts more like a rigid object). That there were not any significant deviations for the smallest value of S is noteworthy considering the fact that that particular sheet buckled to such a large extent that the hard–sphere interactions were engaged. Of course, this is also a symptom of the discretization employed as a finer discretization at the same value of S (i.e., larger number of beads) would allow for higher modes of buckling without hard-sphere contact. Nonetheless, Fig. 3 indicates that even sheets that deviate strongly from planarity exhibit a trajectory that can be approximated well by Jeffery's equations, at least for these initial conditions.

Fig. 3 Angle formed between the flow-vorticity plane and the line passing through the two ends of a sheet in the flow-gradient plane with initial orientation (ϕ,θ) = (0°,5°) and for various dimensionless bending rigidities, S. The prediction of Jeffery's equations for an infinitely thin oblate spheroid is shown with a black dotted line. See Fig. 1 for snapshots of a (buckling) sheet's orientation evolving in time.

Quasi-static elasticity model

A relatively simple 1D continuum model of elasticity can be used to predict the kind of buckling observed for sheets with initial orientations of ϕ = 0°. It is quasi-static in the sense that we assume the forces acting upon the sheet are solely a function of the angle θ(t) and that this angle is determined independently by the Jeffery orbit. This assumption seems to be particularly reasonable in light of the results presented in Fig. 3. Equating moments (see Appendix) and approximating the hydrodynamic stress on the sheet as that of an unperturbed linear flow (i.e., a local analysis), one arrives at the following quasi-static governing equation for small deviations, w, of the sheet from planarity along its length, x, in the flow-gradient plane:where S is the dimensionless bending rigidity from eqn (7), ŵ = w/L, x̂ = (L − x)/L, and L is the length of the sheet from the center (x̂ = 0) to the tip (x̂ = 1). q_⊥, q_‖, and are derived from stresses that are perpendicular or parallel to the sheet and, for a hexagonal sheet, are given by:

q_⊥(x̂) = f(x̂)x̂ sin² θ(‖ṗ(θ)‖/[small gamma, Greek, dot above] − 1) (10)

q_‖(x̂) = f(x̂)x̂ sin θ cos θ (11)

The piecewise nature of these expressions follows from the change in width of the hexagonal sheet along its length in the flow-gradient plane (see Appendix). We imposed the boundary conditions

ŵ(0) = ŵ′′(0) = ŵ′′(1) = ŵ′′′(1) = 0, (13)

where the primes indicate derivatives. These conditions represent symmetry of the sheet across x̂ = 0 and a free end with no applied moment or force, respectively. Similar quasi-static force balances have been applied to model the behavior of filaments subject to large deformations in flow.^12,13,71 Our model differs from these in that it accounts for heterogeneity of the hydrodynamic and elastic forces along the sheet in the flow-gradient plane and neglects the Lagrange multipliers typically employed to enforce inextensibility. Because the model is meant to identify the first deviations from flatness of the sheet only, as done in the classical analysis of Euler buckling, we believe these approximations are appropriate. It is also worth mentioning that, while this work considers hexagonal sheets, the above equations could be altered in a straightforward manner to model sheets of other shapes via modification of the function f.

Unlike Euler buckling, eqn (8) features an inhomogeneous term due to stresses perpendicular to the sheet. However, the solution to eqn (8) can easily be split into a homogeneous and particular solution as ŵ = ŵ^h + ŵ^p, indicating that it is the homogeneous part that is responsible for buckling instabilities. The model can be solved via a basis function expansion in ultraspherical polynomials⁷² as

Aw^h = S⁻¹Bw^h, (14)

where A represents the fourth-order derivative operator, B represents the right-hand side‡ of eqn (8), and w^h is a vector containing the coefficients of the polynomials. The boundary conditions are contained in the top four rows of A, and the top four rows of B are accordingly set to 0. Eqn (14) represents a generalized eigenvalue problem, where S (or S⁻¹) is a generalized eigenvalue. Solving eqn (14) with a basis of 100 polynomials yielded a series of buckling modes and buckling eigenvalues. The first five of these eigenvalue-mode pairs can be seen in Table 2.

Table 2 Predicted values of S and mode shapes corresponding to buckling modes of a hexagonal sheet with orientation ϕ = 0° at the angle of maximum in-plane stress, θ = π/4

S × 10^5	534.92	161.20	78.93	45.09	28.88
Mode

---

The success of this simple 1D model for sheets with initial orientations of ϕ = 0° and θ = 5° can be seen Fig. 4. The maximum in-plane stress is exerted on a flat sheet at an angle of θ = π/4. Neglecting transient effects, we used this value of θ to predict the growth of buckling modes before they ultimately decayed as the sheet reached the stable flat state at long times. Indeed, these values of S seem to faithfully predict the different buckling regimes observed. Sheets with varying dimensionless bending rigidity exhibited different modes of buckling at the apex of their trajectories (i.e., when they were most vertical and the point in time at which bending energy was approximately maximized). Some snapshots of different sheets within these different buckling regimes are shown in Fig. 4 along with images of their mean curvature, H, which was calculated with the “cotangent” formula that is commonly used in computational geometry.⁷⁴ The different modes attained are clearly apparent in these mean curvature diagrams.

Fig. 4 Bending energy over time as a function of dimensionless bending rigidity, S, for an initial orientation of (ϕ,θ) = (0°,5°). Vertical lines correspond to the first ten predicted buckling transitions from the quasi-static model at θ = π/4, the angle of maximum in-plane stress. Snapshots of sheets with different values of S and at dimensionless time ([small gamma, Greek, dot above]t = 11.5) corresponding to those of points labeled A–D in the bending energy diagram are shown with (signed) mean curvature, H, over the intrinsic coordinates of the sheet.

The bending energy, or the maximal bending energy, as a function of S shown in Fig. 5 is quite nonlinear and is inherently not described by this 1D quasi-static model. In fact, for the same reason, it is somewhat remarkable that the model does capture the different buckling modes of such a complex dynamical system. Some of this complicated dependence on S for the smallest values explored (S < 5 × 10⁻⁵), though, can be attributed to the interplay of buckling, strong hydrodynamic interactions, and self-avoidance via hard-sphere interactions, which, along with the smallest length scale of the discretization, places a bound on attainable bending energies. However, for values of S above the first predicted buckling transition (S = 5.35 × 10⁻³), the bending energy scales linearly with κ, as sheets do not buckle and are only perturbed slightly away from the flat state.

Fig. 5 Bending energy as a function of dimensionless bending rigidity, S, for an initial orientation of (ϕ,θ) = (0°,5°) at [small gamma, Greek, dot above]t = 11.4, the point at which sheets are oriented vertically to the flow. Vertical lines correspond to the first ten predicted buckling transitions from the quasi-static model at θ = π/4, the angle of maximum in-plane stress.

It is also interesting to note that the process of flattening out after [small gamma, Greek, dot above]t ≈ 11.4 proceeded via 2D modes of motion which likely cannot be well described by any 1D model (see Movie S2, for example, ESI†). This lack of time reversal symmetry is distinct from Jeffery orbits of rigid spheroidal particles, which are invariant under the transformation (t,L) → (−t, −L), and is also exhibited by the bending energy diagram in Fig. 4, for bending energy decreases more quickly during flattening than it increases during buckling.

4 Chaos and tumbling

In experimental systems, sheets may not be perfectly oriented with ϕ = 0°, so we sought to explore the dynamical behavior of sheets with different initial orientations ϕ ranging from 0° to 90°. To quantify such dynamical behavior, we focused on two summary statistics: the mean orientation and the orientational covariance matrix. The normal vector at each point of a regular surface lives on the unit sphere, S². Thus, the mean orientation of a sheet can be calculated via the weighted Fréchet mean,⁷⁵where the sum is over all triangles (i.e., all groups of three close-packed beads) of the mesh, n̂_i is the normal of triangle i, the weight w_i is the area of triangle i, and dist : (x,y) → arccos(x·y) is the Riemannian distance between two points on the sphere. Note that it is important to ensure that the normals are consistently oriented across the whole mesh. Additionally, with these choice of weights, the sum serves as a kind of finite-volume approximation of the continuous integral over a smooth, continuous surface. Like the typical arithmetic mean in Euclidean space, the Fréchet mean is a point on a manifold (in this case the unit sphere) that minimizes the squared distance to a set of points. The (weighted) orientational covariance matrix is given by:where log_x : S² → T_xS² is the inverse of the exponential map and maps points on the sphere onto the tangent plane at x. Specifically, log_x is calculated aswhere (I − xx^T)(n̂ − x) represents an orthogonal projection of the vector (n̂ − x) in Euclidean space (not intrinsically on the sphere). One can think of this process visually as unwrapping the path on the sphere between two points like an inextensible string onto the tangent plane of the first point, all while maintaining the string's original direction. The covariance matrix C̄ represents the spread of the orientations of the sheet about its mean orientation. A sheet that is flat would exhibit zero spread since all the normals would be identical, whereas a crumpled sheet would exhibit quite a bit of variance about the mean orientation since the normals across the sheet would point in different directions§.

Many numerical simulations of long duration (up to [small gamma, Greek, dot above]t = 1000) were conducted with sheets of varying bending rigidities and initial angles ϕ. Several different dynamical behaviors were observed, as shown in Fig. 6. Movies of representative trajectories are also included in the ESI.† Some sheets (shown with black dots) very closely tracked the trajectory predicted by the Jeffery orbit of an infinitely thin spheroid with an equivalent initial orientation. Specifically, we denote sheet trajectories as “quasi-Jeffery” if |n̄(t)·p(t)| > 0.99 for all times t (i.e., the mean normal is sufficiently close to the Jeffery orbit). Unsurprisingly, stiffer sheets (higher S values) exhibit a wider range of initial orientations that result in quasi-Jeffery trajectories. It is also important to note that the cutoff value of 0.99 did not significantly affect labeling. In other words, sheets that deviated from Jeffery orbits did so considerably. Although not explored in this work, we believe that the 2D buckling patterns exhibited by quasi-Jeffery sheets with ϕ ≠ 0° would likely be described well by a similar 2D quasi-static elasticity model.

Fig. 6 (left) Dynamical states of the sheet as a function of dimensionless bending rigidity, S, and initial angle ϕ. “QJ” denotes “quasi-Jeffery” (close to a Jeffery orbit), “IT” denotes “initial tumbling” (eventually reaching the stable flat state by [small gamma, Greek, dot above]t = 1000), and “CT” denotes “continuous tumbling”. (right) Example snapshots (every [small gamma, Greek, dot above]t = 12.5 units of time) of continuously tumbling, initially tumbling, and quasi-Jeffery trajectories with initial orientation ϕ = 54° and S values of 3.08 × 10⁻⁵, 6.15 × 10⁻⁴, and 3.08 × 10⁻³, respectively. The flow (1), gradient (2), and vorticity (3) directions are depicted at the bottom. Movies of representative trajectories are included in the ESI.†

Sheets that were not quasi-Jeffery but eventually reached the stable flat state by [small gamma, Greek, dot above]t = 1000 were labeled in Fig. 6 as “initial tumbling”. The rest of the sheets, which never reached the stable flat state, were denoted as “continuous tumbling”. Among the sheets that were labeled as initial tumbling, some exhibited erratic tumbling at the beginning of the trajectory while others flipped after long periods of time of ostensibly approaching the flat state (see S = 6.15 × 10⁻⁴ at ϕ = 53° in Fig. 7 as an example). Therefore, this “initial tumbling” regime can be viewed as a transitional regime between quasi-Jeffery and continuously tumbling trajectories. Some apparent outliers can be found in top-right corner of Fig. 6 for large values of S and values of ϕ near 90°. Some of these sheets oscillated back and forth in the flow seemingly indefinitely and in a stable manner. However, we believe that a duration of [small gamma, Greek, dot above]t = 1000 was not long enough to observe what would likely be a return to the stable flat state. That quasi-Jeffery orbits occurred for the inherently unstable initial orientation of ϕ = 90° around a value of S that corresponds to the first predicted buckling transition at ϕ = 0° is likely not coincidental and warrants further investigation.

Fig. 7 Crumpling and chaos in some representative examples of sheets. (A) The sum and difference of the eigenvalues, λ₁ and λ₂, of the orientational covariance matrix, C̄, for three different values of dimensionless bending rigidity, S, and three different initial orientation angles, ϕ. (B) Mean orientation as a function of time (represented by color) for S = 3.08 × 10⁻⁵ at slightly different initial orientation angles. The flow (1), gradient (2), and vorticity (3) directions are depicted at the bottom. A small change in initial orientation results in drastically different aperiodic trajectories, which is the hallmark of chaos.

The sheets that did tumble did so in a chaotic way without any regular periodicity. This behavior is similar to that observed with flexible filaments subjected to shear flow.^17–19Fig. 7 shows summary statistics for representative sheets of various values of S with similar initial angles ϕ. Namely, the eigenvalues, λ₁ and λ₂, of the orientational covariance matrix are presented since they are invariant with respect to rotations or a change of basis. Their sum (equal to tr(C̄)) is representative of the total variance of normals across the surface of the sheet, and their difference is representative of the degree of anisotropy of the distribution. Fig. 8 illustrates the physical meaning of these eigenvalues as they relate to the geometry of the sheet. For example, moderately uniform crumpling can be seen for the softest sheet (black curves) at ϕ = 53° and ϕ = 54° in Fig. 7. Additionally, the isolated spike during the initially tumbling trajectory of the stiffest sheet (red curve) at ϕ = 53° in Fig. 7 represents a highly anisotropic transient crease that occurs as the sheet flips over itself before settling into the flat state.

Fig. 8 Illustration of the physical meaning of the eigenvalues of the orientational covariance matrix C̄. An isotropically crumpled sheet is associated with two positive eigenvalues, a perfectly flat sheet is associated with two zero eigenvalues, and an anisotropically creased sheet is associated with one positive eigenvalue and another close to zero.

Several conclusions can be gleaned from the data presented in Fig. 7. First, it is clear that softer sheets (smaller S values) generally exhibit greater orientational variance than stiffer sheets, both for tumbling trajectories and for the quasi-Jeffery trajectory (ϕ = 4°) presented. This trend is entirely expected: softer sheets manifest greater degrees of crumpling, which directly corresponds to a greater variance of normals about the mean orientation. Second, it is apparent that the difference of the eigenvalues, λ₁ − λ₂, as a fraction of their sum was in general smallest for the softest sheet. This behavior indicates that the soft sheet tends to crumple more isotropically (similar crumpling in all directions), whereas stiffer sheets crumple more anisotropically (like folding a crease along a single direction). Finally, it is clear for ϕ = 53° and ϕ = 54°—the angles that produced tumbling—that a small change in initial conditions yielded drastically different trajectories, as evidenced by the pattern of the orientational covariance and the mean orientations themselves (panel B of Fig. 7). In fact, these continuously tumbling trajectories diverged away from each other exponentially quickly. This, of course, is the classical signature of chaos (see Fig. 10 and the Appendix for Lyapunov exponent estimates). It is important to note that although the sheets tumble in a circular way and the covariance seems to “ebb and flow” over time, we did not detect any significant signature of regular periodicity in power spectral densities. One can also see that for the sheets that initially tumble, the time at which the stable flat state is attained is quite erratic. Overall, the trends discussed for these few sheets hold for the many others that were studied and depicted as single dots in the state diagram of Fig. 6.

Interestingly, chaotic trajectories are often associated with crumpled configurations (i.e., higher orientational variances). This connection is not obvious a priori. Perhaps it is the more crumpled states of softer sheets that allow stronger and more complex fluid–bead and bead–bead interactions to occur and give rise to sensitive, chaotic dynamics.

5 Conclusions

The dynamical behavior of athermal 2D sheets immersed in a shear flow at low Reynolds number was investigated via immersed boundary simulations with a model semiflexible sheet at high Föppl–von Kármán numbers (i.e., softer bending relative to stretching modes of motion). The main governing dimensionless parameter of the system is the dimensionless bending rigidity, S = κ/(πη[small gamma, Greek, dot above]L³). Our findings on the behavior of athermal sheets can be summarized succinctly as follows:

1. For flat sheets initially oriented with a normal in the flow-gradient plane, transient buckling occurs and can be predicted as a function of S quite accurately using a simple 1D elasticity model.

2. Smaller values of S can result in chaotic, continuously tumbling trajectories, but not for all initial orientations.

3. Chaotic trajectories are associated with crumpled conformations.

4. Sheets that do not tumble chaotically (generally those lying near the flow-vorticity plane or large values of S) are described well by the equivalent Jeffery orbit of a rigid, spheroidal particle.

This delineation of the dynamical regimes of athermal sheets should inform the design of solution processing protocols of flexible 2D materials, where crumpling or buckling may or may not be desired for different applications. Specifically, future directions and applications include better understanding the influence of shear-induced morphological changes of dispersed nanosheets on bulk rheological properties (e.g., shear-dependent viscosity), the potential migration of the sheets in unbounded shear or in other flow geometries, and the design of precision flow systems to tune nanomaterial conformations. Additionally, the use of the “intrinsic” summary statistics employed in this work (mean orientation and the orientational covariance matrix) may prove to be useful in future theoretical and experimental studies of related systems.

Conflicts of interest

There are no conflicts to declare.

Appendix

A Derivation of quasi-static elasticity model

First, let us consider a 1D rod (see Fig. 9) under the assumptions of classical Euler–Bernoulli beam theory with a flexural rigidity of EI, where E is the 3D Young's modulus and I is the second moment of area of the cross-section.⁵³ The equation for a small deflection, w, of the rod in the y direction is given bywhere the right-hand side is given by the moments on the rod in the perpendicular and parallel directions (similar to a self-buckling analysis) and the functions q are given by

q_x = αη[small gamma, Greek, dot above]cos(θ)(L − x)sin(θ) (19)

q_y = αη([small gamma, Greek, dot above] − ‖ṗ‖)sin(θ)(L − x)sin(θ). (20)

Here, it is assumed that the force per unit length is proportional to the viscosity, η, and the velocity of an unperturbed linear shear flow with a dimensionless constant α that depends on the exact shape of the rod. The angle θ is approximated independently by the equivalent Jeffery orbit, p, of a rigid object, which is also used to adjust the net moment in the perpendicular direction. In this analysis, it is also assumed that there is no transient response, hence the “quasi-static” nature of the model.

Fig. 9 Schematic of a rod (or sheet) in shear flow.

Fig. 10 The norm of the difference between the vertex coordinates of two continuously tumbling sheets with initial conditions of (ϕ,θ = 5°) and (ϕ + 1°,θ = 5°) over time. For all trajectories shown, S = 3.08 × 10⁻⁵. An exponential divergence in trajectories can be seen for short times.

Differentiation once with respect to x yields

and differentiating once more yields The boundary conditions are The first two represent the free end of the rod at x = 0, which is subjected to no moment or force. The latter two reflect the symmetry of the behavior of the rod about the center streamline.

Nondimensionalizing the equation using ŵ = w/L and x̂ = (L − x)/L yields:

whereis the elasto-viscous number.¹⁶

For a 2D sheet extending in and out of the plane of the paper (the z direction), deforming purely in a 1D way, and experiencing stresses independent of the z direction, EI can be considered an effective modulus proportional to , where κ is the bending rigidity of the 2D sheet. Likewise, q_x and q_y can be considered as effective 1D stresses where for some constant β (with units of length⁻¹) that depends on the hydrodynamics of the sheet. Thus, the relevant dimensionless quantity becomes S, defined as

instead of B. It is interesting to note that for varying widths, W, in the z direction, S scales like W/L⁴. For widths approaching the bead radius, the elasto-viscous L⁴ scaling for filaments is recovered.

A.1 Determination of constants for bead model. First, as mentioned in the main text, the bending rigidity for a triangulated sheet with dihedral forces is times greater than the bending rigidity of the equivalent continuum sheet.⁵⁹ Thus, in mapping S to simulation data, the numerator of S should be multiplied by . Additionally, one can calculate β based on the discretization of beads employed. Since the beads (of radius a = 1) are 2D close-packed, the stress on the sheet due to a velocity, u, in a given area, A, can be approximated aswhere is the density of close-packed spheres in 2D. Therefore, β for our model must be which does indeed have units of length⁻¹. The integral of β contributes a factor of a/L to the definition of S. However, where the proportionality factor of L/a reflects the density of dihedrals on the discrete sheet and derives from the fact that there are (2L/2a) bonds in the z direction. This factor of L/a cancels out the contribution of a/L from β. With these considerations, S is independent of the bead length scale a and is given by

A.2 Accounting for changing width. As one travels along the x direction, the width of a hexagonal sheet in the z direction changes. That is, the width of a hexagonal sheet as a fraction of its total width is given byOne can account for this change of width in the effective 1D stress constant α and flexural rigidity EI by simply multiplying β and κ respectively by this (dimensionless) piecewise function, rendering both quantities x-dependent. The stresses q_y and q_x, in turn, gain additional x dependence due to this changing width. Now,

q_x = αf(x)η[small gamma, Greek, dot above]cos(θ)(L − x)sin(θ) (31)

q_y = αf(x)η([small gamma, Greek, dot above]− ‖ṗ‖)sin(θ)(L − x)sin(θ). (32)

It is these expressions that lead to eqn (9)–(12) in the main text.

A.3 Verification of scaling. Fig. 11 is essentially the same as Fig. 5 but shows the bending energy attained when vertically oriented for differently sized sheets of lengths L = 48, 52, and 58 (with numbers of beads equal to 1801, 2107, and 2611, respectively). In these additional simulations, the bond stiffness was varied in order to maintain a constant FvK number. As one can see, the proposed L³ scaling of the dimensionless bending rigidity, S, indeed holds as all of the curves collapse onto each other.

Fig. 11 Bending energy as a function of dimensionless bending rigidity, S, for an initial orientation of ϕ = 0° and θ = 5° at [small gamma, Greek, dot above]t = 11.4, the point at which sheets are oriented vertically to the flow. The different colors represent differently sized sheets in the simulation. Vertical lines correspond to the first ten predicted buckling transitions from the quasi-static model at θ = π/4, the angle of maximum in-plane stress.

B Maximal Lyapunov exponents

Fig. 10 shows the norm of the difference between the vertex coordinates of sheets of dimensionless rigidity S = 3.08 × 10⁻⁵ with slightly different initial conditions (i.e., where ϕ differs by 1°). This value of S corresponds to trajectories that are continuously tumbling for all of the orientations featured in Fig. 10. At each timestep, the coordinates were translated in space so that the center of mass of the sheet was located at the origin. Given that sheets do not stretch to a significant extent, this translation ensures that the coordinates of the sheet live in a compact state space. One can see that, at least for short times, the small perturbation in initial orientation ϕ causes trajectories to diverge from each other exponentially. Maximal Lyapunov exponents⁷⁶ can be estimated by fitting the data to the functional form

‖Δx_t‖₂ = ‖Δx₀‖₂e^{λ_max[small gamma, Greek, dot above]t}. (33)

Here, Δx_t is the difference between the vertex coordinates (after subtracting the center of mass) of two sheets with different initial conditions at time t. Performing a least-squares linear regression on the logarithm of the norm differences at times 2.25 ≤ [small gamma, Greek, dot above]t ≤ 4.75 for the three groups of initial angles shown in Fig. 10 yields maximal Lyapunov exponents of 0.31, 0.41, and 0.56 for ϕ near 36°, 54°, and 72°, respectively. Thus, the maximal Lyapunov exponents are dependent on initial orientation and increase as ϕ increases. It should be noted that these Lyapunov exponents are dimensionless; in real units of time, one can see that they scale with the shear rate.

C Movies

Movies of sheet dynamics can be found in the ESI,† and include the following:

S1. Buckling of a sheet with S = 3.08 × 10⁻⁵ and initial orientation ϕ = 0°.

S2. Quasi-Jeffery trajectory of a sheet with S = 3.08 × 10⁻⁵ and initial orientation ϕ = 4°.

S3. Continuously tumbling trajectory of a sheet with S = 3.08 × 10⁻⁵ and initial orientation ϕ = 54°.

S4. Buckling of a sheet with S = 6.15 × 10⁻⁴ and initial orientation ϕ = 0°.

S5. Quasi-Jeffery trajectory of a sheet with S = 6.15 × 10⁻⁴ and initial orientation ϕ = 4°.

S6. Initially tumbling trajectory of a sheet with S = 6.15 × 10⁻⁴ and initial orientation ϕ = 54°.

S7. Quasi-Jeffery trajectory (and imperceptible buckling) of a sheet with S = 3.08 × 10⁻³ and initial orientation ϕ = 0°.

S8. Quasi-Jeffery trajectory of a sheet with S = 3.08 × 10⁻³ and initial orientation ϕ = 54°.

Acknowledgements

K. S. S. was supported by the Department of Energy Computational Science Graduate Fellowship program under grant No. DE-FG02-97ER25308. M. S. S. was supported by the Department of Energy, Office of Science, Basic Energy Sciences under grant No. DE-FG02-08ER46488 Mod 0008. J. W. S. was supported by NSF Career Award No. CBET-1554398. The authors would like to thank M. Ekiel-Jeżewska, G. McKinley, P. Doyle, and N. Fakhri for helpful discussions.

Notes and references

1. P. G. de Gennes, J. Chem. Phys., 1974, 60, 5030–5042.
2. R. G. Winkler, Phys. Rev. Lett., 2006, 97, 128301.
3. C. M. Schroeder, R. E. Teixeira, E. S. G. Shaqfeh and S. Chu, Phys. Rev. Lett., 2005, 95, 018301.
4. D. E. Smith, H. P. Babcock and S. Chu, Science, 1999, 283, 1724–1727.
5. M. Harasim, B. Wunderlich, O. Peleg, M. Kröger and A. R. Bausch, Phys. Rev. Lett., 2013, 110, 108302.
6. A. Nikoubashman, A. Milchev and K. Binder, J. Chem. Phys., 2016, 145, 234903.
7. A. Nikoubashman and M. P. Howard, Macromolecules, 2017, 50, 8279–8289.
8. C.-C. Huang, G. Sutmann, G. Gompper and R. G. Winkler, Europhys. Lett., 2011, 93, 54004.
9. I. S. Dalal, N. Hoda and R. G. Larson, J. Rheol., 2012, 56, 305–332.
10. T. Munk, O. Hallatschek, C. H. Wiggins and E. Frey, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2006, 74, 041911.
11. P. S. Doyle, E. S. G. Shaqfeh and A. P. Gast, J. Fluid Mech., 1997, 334, 251–291.
12. A. Lindner and M. Shelley, in Fluid–Structure Interactions in Low-Reynolds-Number Flows, ed. C. Duprat and H. A. Stone, Royal Society of Chemistry, 2015, ch. 5, pp. 168–192.
13. Y.-N. Young and M. J. Shelley, Phys. Rev. Lett., 2007, 99, 058303.
14. V. Kantsler and R. E. Goldstein, Phys. Rev. Lett., 2012, 108, 038103.
15. N. Quennouz, M. Shelley, O. du Roure and A. Lindner, J. Fluid Mech., 2015, 769, 387–402.
16. Y. Liu, B. Chakrabarti, D. Saintillan, A. Lindner and O. du Roure, Proc. Natl. Acad. Sci. U. S. A., 2018, 115, 9438–9443.
17. A. M. Słowicka, E. Wajnryb and M. L. Ekiel-Jeżewska, J. Chem. Phys., 2015, 143, 124904.
18. A. M. Słowicka, H. A. Stone and M. L. Ekiel-Jeżewska, Phys. Rev. E, 2020, 101, 023104.
19. J. LaGrone, R. Cortez, W. Yan and L. Fauci, J. Non-Newtonian Fluid Mech., 2019, 269, 73–81.
20. M. J. Shelley and J. Zhang, Annu. Rev. Fluid Mech., 2011, 43, 449–465.
21. M. Argentina and L. Mahadevan, Proc. Natl. Acad. Sci. U. S. A., 2005, 102, 1829–1834.
22. P. Liu, A. T. Liu, D. Kozawa, J. Dong, J. F. Yang, V. B. Koman, M. Saccone, S. Wang, Y. Son, M. H. Wong and M. S. Strano, Nat. Mater., 2018, 17, 1005–1012.
23. F. Tardani, W. Neri, C. Zakri, H. Kellay, A. Colin and P. Poulin, Langmuir, 2018, 34, 2996–3002.
24. D. Parviz, S. D. Metzler, S. Das, F. Irin and M. J. Green, Small, 2015, 11, 2661–2668.
25. J. W. Colson and W. R. Dichtel, Nat. Chem., 2013, 5, 453–465.
26. P. Payamyar, B. T. King, H. C. Öttinger and A. D. Schlüter, Chem. Commun., 2016, 52, 18–34.
27. Y. Xu and M. J. Green, J. Chem. Phys., 2014, 141, 024905.
28. Y. Xu and M. J. Green, J. Polym. Sci., Part B: Polym. Phys., 2015, 53, 1247–1253.
29. S. B. Babu and H. Stark, Eur. Phys. J. E: Soft Matter Biol. Phys., 2011, 34, 136.
30. E. Frey and D. L. Nelson, J. Phys., 1991, 1, 1715–1757.
31. S. Dutta and M. D. Graham, Soft Matter, 2017, 13, 2620–2633.
32. Y. Yu and M. D. Graham, Soft Matter, 2021, 17, 543–553.
33. J. Stafford, A. Patapas, N. Uzo, O. K. Matar and C. Petit, AIChE J., 2018, 64, 3246–3276.
34. A. Ambrosi and M. Pumera, Chem. Soc. Rev., 2018, 47, 7213–7224.
35. J. N. Coleman, M. Lotya, A. O'Neill, S. D. Bergin, P. J. King, U. Khan, K. Young, A. Gaucher, S. De, R. J. Smith, I. V. Shvets, S. K. Arora, G. Stanton, H.-Y. Kim, K. Lee, G. T. Kim, G. S. Duesberg, T. Hallam, J. J. Boland, J. J. Wang, J. F. Donegan, J. C. Grunlan, G. Moriarty, A. Shmeliov, R. J. Nicholls, J. M. Perkins, E. M. Grieveson, K. Theuwissen, D. W. McComb, P. D. Nellist and V. Nicolosi, Science, 2011, 331, 568–571.
36. V. Nicolosi, M. Chhowalla, M. G. Kanatzidis, M. S. Strano and J. N. Coleman, Science, 2013, 340, 1226419.
37. E. Varrla, C. Backes, K. R. Paton, A. Harvey, Z. Gholamvand, J. McCauley and J. N. Coleman, Chem. Mater., 2015, 27, 1129–1139.
38. S. A. Mallory, C. Valeriani and A. Cacciuto, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2015, 92, 012314.
39. T. Gibaud, C. N. Kaplan, P. Sharma, M. J. Zakhary, A. Ward, R. Oldenbourg, R. B. Meyer, R. D. Kamien, T. R. Powers and Z. Dogic, Proc. Natl. Acad. Sci. U. S. A., 2017, 114, E3376–E3384.
40. F. Ding, J. Liu, S. Zeng, Y. Xia, K. M. Wells, M.-P. Nieh and L. Sun, Sci. Adv., 2017, 3, e1701212.
41. P. Davidson, C. Penisson, D. Constantin and J.-C. P. Gabriel, Proc. Natl. Acad. Sci. U. S. A., 2018, 201802692.
42. A. Laskar, O. E. Shklyaev and A. C. Balazs, Sci. Adv., 2018, 4, eaav1745.
43. C. Kamal, S. Gravelle and L. Botto, Nat. Commun., 2020, 11, 2425.
44. M. Pezzulla, E. F. Strong, F. Gallaire and P. M. Reis, Phys. Rev. Fluids, 2020, 5, 084103.
45. A. R. Klotz, B. W. Soh and P. S. Doyle, Proc. Natl. Acad. Sci. U. S. A., 2020, 117, 121–127.
46. R. Nicklow, N. Wakabayashi and H. G. Smith, Phys. Rev. B: Condens. Matter Mater. Phys., 1972, 5, 4951–4962.
47. N. Lindahl, D. Midtvedt, J. Svensson, O. A. Nerushev, N. Lindvall, A. Isacsson and E. E. B. Campbell, Nano Lett., 2012, 12, 3526–3531.
48. M. K. Blees, A. W. Barnard, P. A. Rose, S. P. Roberts, K. L. McGill, P. Y. Huang, A. R. Ruyack, J. W. Kevek, B. Kobrin, D. A. Muller and P. L. McEuen, Nature, 2015, 524, 204–207.
49. P. Poulin, R. Jalili, W. Neri, F. Nallet, T. Divoux, A. Colin, S. H. Aboutalebi, G. Wallace and C. Zakri, Proc. Natl. Acad. Sci. U. S. A., 2016, 113, 11088–11093.
50. J.-W. Jiang, Z. Qi, H. S. Park and T. Rabczuk, Nanotechnology, 2013, 24, 435705.
51. U. Seifert, Adv. Phys., 1997, 46, 13–137.
52. Polymer Handbook, ed. J. Brandrup, E. H. Immergut and E. A. Grulke, Wiley, New York, 4th edn, 1999.
53. L. D. Landau and E. M. Lifshitz, Theory of Elasticity, Elsevier, Amsterdam, 3rd edn, 2008, vol. 7.
54. D. Nelson and L. Peliti, J. Phys., 1987, 48, 1085–1092.
55. M. J. Bowick, A. Košmrlj, D. R. Nelson and R. Sknepnek, Phys. Rev. B, 2017, 95, 104109.
56. G. Wang, Z. Dai, J. Xiao, S. Feng, C. Weng, L. Liu, Z. Xu, R. Huang and Z. Zhang, Phys. Rev. Lett., 2019, 123, 116101.
57. X. Bian, S. Litvinov and P. Koumoutsakos, Comput. Methods Appl. Mech. Eng., 2020, 359, 112758.
58. A. Guckenberger and S. Gekle, J. Phys.: Condens. Matter, 2017, 29, 203001.
59. G. Gompper and D. M. Kroll, J. Phys., 1996, 6, 1305–1320.
60. D. M. Heyes and J. R. Melrose, J. Non-Newtonian Fluid Mech., 1993, 46, 1–28.
61. J. Rotne and S. Prager, J. Chem. Phys., 1969, 50, 4831–4837.
62. H. Yamakawa, J. Chem. Phys., 1970, 53, 436–443.
63. A. M. Fiore, F. Balboa Usabiaga, A. Donev and J. W. Swan, J. Chem. Phys., 2017, 146, 124116.
64. J. W. Swan and G. Wang, Phys. Fluids, 2016, 28, 011902.
65. A. M. Fiore and J. W. Swan, J. Chem. Phys., 2018, 148, 044114.
66. E. Wajnryb, K. A. Mizerski, P. J. Zuk and P. Szymczak, J. Fluid Mech., 2013, 731, R3.
67. S. Kim and S. J. Karrila, Microhydrodynamics: Principles and Selected Applications, Dover Publications, Mineola, N.Y, 2005.
68. J. A. Anderson, J. Glaser and S. C. Glotzer, Comput. Mater. Sci., 2020, 173, 109363.
69. G. B. Jeffery, Proc. R. Soc. London, Ser. A, 1922, 102, 161–179.
70. E. J. Hinch and L. G. Leal, J. Fluid Mech., 1979, 92, 591–607.
71. L. E. Becker and M. J. Shelley, Phys. Rev. Lett., 2001, 87, 198301.
72. S. Olver and A. Townsend, SIAM Rev., 2013, 55, 462–489.
73. L. Trefethen, Approximation Theory and Approximation Practice, SIAM, 2013.
74. M. Meyer, M. Desbrun, P. Schröder and A. H. Barr, in Visualization and Mathematics III, ed. G. Farin, H.-C. Hege, D. Hoffman, C. R. Johnson, K. Polthier, H.-C. Hege and K. Polthier, Springer Berlin Heidelberg, Berlin, Heidelberg, 2003, pp. 35–57.
75. X. Pennec, J. Math. Imaging Vis, 2006, 25, 127.
76. S. H. Strogatz, Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering, Perseus Books, Cambridge, Mass, 2nd edn, 2001.

---

Footnotes

† Electronic supplementary information (ESI) available: Movies of sheet dynamics. See DOI: 10.1039/d0sm02184a

‡ The piecewise nature of eqn (8) was not treated specially. Smooth polynomials of high degree can approximate piecewise continuous functions with minimal error.⁷³

§ There is some ambiguity if the sheet curls around itself more than 180° since the Riemannian distance would not reflect this (i.e., distances between two points on the sphere are bounded between 0 and π). Such a scenario, though, is only possible with highly crumpled sheets, which would still be reflected in these statistics as a large variance about the mean orientation.

---