Buckling, Crumpling, and Tumbling of Semiflexible Sheets in Simple Shear Flow

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 in the flow-gradient 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, which is in distinct contrast to polymer chains.


INTRODUCTION
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 [1], and subsequent theoretical [2] and experimental work [3][4][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][7][8][9][10][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 [12]. Namely, fibers have been found to exhibit buckling [13][14][15] and morphological transitions [16] as well as periodic and chaotic dynamical trajectories [17,18].
At high Reynolds number, one example of fluidstructure interactions with 2D objects that has attracted researchers' attention is the flapping of flags [19,20]. 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 [21], graphene oxide [22,23], and 2D polymers [24,25] in advanced technologies. Xu and Green [26,27] looked at sheets under shear and biaxial extensional flow, and Babu and Stark [28] have studied tethered sheets in fluids, confirming predicted scaling laws of Frey and Nelson [29] regarding thermal fluctuations in the presence of hydrodynamic interactions. Additionally, Dutta and Graham [30] have studied Miurapatterned sheets and observed various interesting dynamical regimes involving periodic tumbling, unfolding, and quasiperiodic limit cycles. Most recently, Yu and Graham [31] studied "compact-stretch" transitions of elastic sheets under extensional flow. However, there is still * jswan@mit.edu 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 [32][33][34][35] 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.) [35,36]. Furthermore, experimental studies involving 2D materials whose functions depend significantly on coupled fluid motion are becoming more prevalent [22,[37][38][39][40][41][42][43].

Material
Bending Energy ( kBT | 298 K ) Graphene 40 [44], 280 [45] We consider thin, isotropic, semiflexible sheets of size L immersed in a fluid at low Reynolds number. Sheets can be characterized by a bending stiffness, κ, 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 Y L 2 /κ. From the classical theory of thin plates [51], both κ and Y are also related to the 3D Young's modulus of the material, E, as Y = Eh and κ = Y h 3 /[12(1 − ν 2 )], where h is the thickness of the sheet and ν is the Poisson ratio. As reference, some examples of bending energies are provided in Table I. For graphene, values can vary significantly from the theoretical microscopic bending stiffness (the first value) due to thermal fluctuation-induced stiffening, which is, in turn, experimentally challenging to measure and length-scale dependent [52,53]. It is also worth noting that many nanomaterials are produced containing multiple layers, and it has been shown that mutlilayer van der Waals materials exhibit intermediate behavior between classic thin plates and ideally lubricated stacks of sheets [54].
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 energy to shear energy (a dimensionless ratio we denote as S), we find behavior ranging from deterministic quasi-1D buckling to chaotic tumbling. A simple continuum elasticity model is presented to explain this buckling, and summary statistics of sheet orientation (viz., the mean orientation and the orientational covariance matrix) are constructed and used to analyze chaotic motions.

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). 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 (C 6 ), 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: where κ is the bending energy andn i andn j are (consistently oriented) triangle normals [53,[55][56][57]. Such a model of bending may not be generally applicable for capturing bending forces due to large curvature (under the Helfrich-Canham model, for example) but is particularly well suited for the naturally flat sheets considered here. The equivalent continuum bending energy is given byκ = κ/ √ 3 [57]. 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 time step 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 [58]. Harmonic bonds of the form were applied along each of the edges of the triangulation with stiffness k = 1000 × 6πηγa and r 0 = 2a, where η andγ are the viscosity and shear rate of the surrounding fluid. With these values of bond stiffness and bending energies, the Föppl-von Kármán (FvK) number of the system ranged in magnitude between 10 4 and 10 6 , 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 [59,60]. 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 [61,62]. 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: wherer 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 x ∈ R 3N is a stacked vector of N particle coordinates, U is the total potential energy and (L12)mn =γδm1δn2 is the 3 × 3 velocity gradient tensor. Importantly, advecting the particles as point particles and not as finitely sized, rigid spheres (which would require more sophisticated integrators that operate at the stresslet level [63]) causes the sheet to behave as if it had zero thickness with respect to the shear flow. This lack of thickness breaks the periodicity of Jeffery orbits and introduces stable flat states, 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 T −1 ≈γ 2π h L for small values of the thickness, h, 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.
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. 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 equation 5 via a forward Euler scheme with a timestep of at mostγ∆t = 5 × 10 −4 using a custom plugin adapted from ref. [61] for the HOOMD-blue molecular simulation package [64] on graphics processing units (GPUs). Specifically, all simulations were conducted on NVIDIA GTX 1080 Tis and required thousands of GPU-hours.

BUCKLING
The dynamics of rigid ellipsoidal particles in simple shear flows was first worked out by Jeffery [65,66] with details for non-axisymmetric particles later studied by Hinch and Leal [67]. The differential equation governing the orientation, p, of axisymmetric spheroidal particles iṡ where Ω = (L12 − L T 12 )/2 is the antisymmetric vorticity tensor, E = (L12 + L T 12 )/2 is the symmetric strain-rate tensor, and a1 and a2 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 center streamline of the flow. Note, again, that such a stable state only exists for infinitely flat sheets; sheets with finite thickness would eventually tumble and buckle again periodically. Figure 1 shows snapshots of a particular sheet (κ = 15) buckling during its trajectory. The degree, or mode, of buckling depended on the bending stiffness and will be discussed more below. Figure 2 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 energy, compared to the angle predicted by the Jeffery orbit of an infinitely thin, rigid, oblate spheroid. This dimensionless bending energy -essentially the sheet's equivalent of the "elasto-viscous" number as it is sometimes called in the flexible filament literature [16] -is defined as Here, L is the characteristic radius of the sheet (equal to 58), and the numerical factor of 6/6π 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, Figure  2 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.

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 Figure 2. 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 energy from equation 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 r = − 1 x dy q (y) 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 θ cos θ r (x) = − sin θ cos θ 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 conditionŝ where the primes indicate derivatives. These conditions represent symmetry of the sheet acrossx = 0 and a free end with no applied moment or force, respectively. In-plane tension forces, which are often applied in 1D filament modeling as Lagrange multipliers to maintain inextensibility [12], were neglected but could easily be incorporated. Unlike Euler buckling, equation 8 features an inhomogeneous term due to stresses perpendicular to the sheet. However, the solution to equation 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 [68] as where A represents the fourth-order derivative operator, B represents the right-hand side [69] Table II. 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.
The success of this simple 1D model for sheets with initial orientations of φ = 0 • can be seen Figure 3. 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 stiffness 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 Figure 3 along with images of their mean curvature, H, which was calculated with the "cotangent" formula that is commonly used in computational geometry [70]. The different modes attained are clearly apparent in these mean curvature diagrams.  Vertical lines correspond to the first ten predicted buckling transitions from the quasi-static model at θ = π/4, the angle of maximum in-plane stress.
The bending energy, or the maximal bending energy, as a function of S shown in Figure 4 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 −5 ), 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.32 × 10 −3 ), the bending energy scales linearly with κ, as sheets do not buckle and are only perturbed slightly away from the flat state. It is also interesting to note that the process of flattening out afterγt ≈ 11.4 proceeded via 2D modes of motion which likely cannot be well described by any 1D model (see Movie S2, for example). 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 Figure 3, for bending energy decreases more quickly during flattening than it increases during buckling.

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 2 . Thus, the mean orientation of a sheet can be calculated via the weighted Fréchet mean [71], n = arg min where the sum is over all triangles (i.e., all groups of three close-packed beads) of the mesh,ni is the normal of triangle i, the weight wi 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 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 2 → TxS 2 is the inverse of the exponential map and maps points on the sphere onto the tangent plane at x. 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 its direction. The covariance matrixC represents the spread Figure 5. Left) Dynamical states of the sheet as a function of dimensionless bending stiffness, S, and initial angle φ. "QJ" denotes "quasi-Jeffery" (close to a Jeffery orbit), "IT" denotes initial tumbling (eventually reaching the stable flat state bẏ γt = 1000), and "CT" denotes "continuous tumbling". Right) Example snapshots (everyγ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 −5 , 6.15 × 10 −4 , and 3.08 × 10 −3 , respectively. Movies of representative trajectories are included in the SI.
of the orientations of the sheet about its mean orientation. A sheet that is flat would exhibit zero spread, whereas a crumpled sheet would exhibit quite a bit of variance about the mean orientation [72]. Many numerical simulations of long duration (up toγt = 1000) were conducted with sheets of varying bending energies and initial angles φ. Several different dynamical behaviors were observed, as shown in Figure 5. Movies of representative trajectories are also included in the SI. 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.
Sheets that were not quasi-Jeffery but eventually reached the stable flat state byγt = 1000 were labeled in Figure 5 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 −4 at φ = 53 • in Figure 6 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 Figure 5 for large values of S and values of φ near 90 • . Some of these sheets seemed to oscillate back and forth in the flow indefinitely. However, we believe that a duration ofγt 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.
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,18]. Figure 6 shows summary statistics for representative sheets of various values of S with similar initial angles φ. Namely, the eigenvalues, λ1 and λ2, 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. Several conclusions can be gleaned from the data presented in Figure 6. 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, λ1 − λ2, 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 Figure 6). This, of course, is the classical signature of chaos. 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 Figure 5.
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.

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 energy, S = 6κ/(6πηγL 3 ). 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., sheardependent viscosity) 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.

Determination of Constants for Bead Model
First, as mentioned in the main text, the bending stiffness for a triangulated sheet with dihedral forces is √ 3 times greater than the bending stiffness of the equivalent continuum sheet [57]. Thus, in mapping S to simulation data, the numerator of S should be multiplied by √ 3. 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 as where √ 3π/6 is the density of close-packed spheres in 2D. Therefore, β for our model must be 6π √ 3/(6a), which does indeed have units of length −1 . With these two considerations, S is given by S = 6κ 6πηγL 3 . (A12)

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 by x ≥ L 2 . (A13) One 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 qy and qx, in turn, gain additional x dependence due to this changing width. Now, It is these expressions that lead to equations 9-12 in the main text.