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

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 [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, 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 [22], 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 [29] have studied tethered sheets in fluids, confirming predicted scaling laws of Frey and Nelson [30] regarding thermal fluctuations in the presence of hydrodynamic interactions. Additionally, Dutta and Graham [31] have studied Miura-patterned sheets and observed various interesting dynamical regimes involving periodic tumbling, unfolding, and quasiperiodic limit cycles. Most recently, Yu and Graham [32] 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, 34, 35, 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 [38, 39, 40, 23, 41, 42, 43, 44, 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, $\kappa$, 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^{2}/\kappa$. From the classical theory of thin plates [53], both $\kappa$ and $Y$ are also related to the 3D Young’s modulus of the material, $E$, as $Y=Eh$ and $\kappa=Yh^{2}/[12(1-\nu^{2})]$, where $h$ is the thickness of the sheet and $\nu$ 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 [56]. While these complications and others, such as nonzero slip [43], 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., $\kappa/k_{B}T\gg 1$).

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.

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 $\bigtriangleup_{i}$ and $\bigtriangleup_{j}$ as:

where $\kappa$ is the bending rigidity and $\hat{\mathbf{n}}_{i}$ and $\hat{\mathbf{n}}_{j}$ are (consistently oriented) triangle normals [55, 57, 58, 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 $\tilde{\kappa}=\kappa/\sqrt{3}$ [59].

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 $\Delta 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 [60].

Harmonic bonds of the form

were applied along each of the edges of the triangulation with stiffness $k=1000\times 6\pi\eta\dot{\gamma}a$ and $r_{0}=2a$, where $\eta$ and $\dot{\gamma}$ are the viscosity and shear rate of the surrounding fluid. The continuum 2D Young’s modulus is related to $k$ via the expression $Y=2k/\sqrt{3}$ [55]. 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^{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 [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\times 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 $\hat{\mathbf{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 $\mathbf{x}\in\mathbb{R}^{3N}$ is a stacked vector of $N$ particle coordinates, $U$ is the total potential energy and $(\mathbf{L})_{mn}=\dot{\gamma}\delta_{m1}\delta_{n2}$ is the $3\times 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 [65] 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\to 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 $T^{-1}\approx\frac{\dot{\gamma}}{2\pi}\frac{h}{L}$ 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, 66, 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 $\phi$ about the flow axis and then by an angle $\theta=5^{\circ}$ about the vorticity axis (see Figure 2). Rotating by $\theta$ perturbs the sheet away from the stable state $(\theta,\phi)=(0^{\circ},0^{\circ})$. By symmetry, unique initial conditions are only generated by angles $\phi\in[0^{\circ},90^{\circ}]$. Simulations were conducted by integrating equation 5 via a forward Euler scheme with a timestep of at most $\dot{\gamma}\Delta t=5\times 10^{-4}$ using a custom plugin adapted from ref. [63] for the HOOMD-blue molecular simulation package [68] 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.

The dynamics of rigid ellipsoidal particles in simple shear flows was first worked out by Jeffery [69, 67] with details for non-axisymmetric particles later studied by Hinch and Leal [70]. The differential equation (in terms of matrix-vector products) governing the orientation vector, $\mathbf{p}$, of axisymmetric spheroidal particles is

where $\bm{\Omega}=(\mathbf{L}-\mathbf{L}^{\mathsf{T}})/2$ is the antisymmetric vorticity tensor, $\mathbf{E}=(\mathbf{L}+\mathbf{L}^{\mathsf{T}})/2$ is the symmetric strain-rate tensor, and $a_{1}$ and $a_{2}$ 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 $\phi=0^{\circ}$, 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. Figure 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.

Figure 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 [16] — is defined as

Here, $L$ is the characteristic radius of the sheet (equal to 58), and the numerical factor of $1/\pi$ 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 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.

In experimental systems, sheets may not be perfectly oriented with $\phi=0^{\circ}$, so we sought to explore the dynamical behavior of sheets with different initial orientations $\phi$ ranging from $0^{\circ}$ to $90^{\circ}$. 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 [75],

where the sum is over all triangles (i.e., all groups of three close-packed beads) of the mesh, $\hat{\mathbf{n}}_{i}$ is the normal of triangle $i$, the weight $w_{i}$ is the area of triangle $i$, and $\operatorname{dist}:(\mathbf{x},\mathbf{y})\mapsto\arccos(\mathbf{x}\dotproduct\mathbf{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_{\mathbf{x}}:S^{2}\to T_{\mathbf{x}}S^{2}$ is the inverse of the exponential map and maps points on the sphere onto the tangent plane at $\mathbf{x}$. Specifically, $\log_{\mathbf{x}}$ is calculated as

where $(\mathbf{I}-\mathbf{x}\mathbf{x}^{\mathsf{T}})(\hat{\mathbf{n}}-\mathbf{x})$ represents an orthogonal projection of the vector $(\hat{\mathbf{n}}-\mathbf{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 $\bar{\mathbf{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 ²²2There is some ambiguity if the sheet curls around itself more than $180^{\circ}$ since the Riemannian distance would not reflect this (i.e., distances between two points on the sphere are bounded between $0$ and $\pi$). 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..

Many numerical simulations of long duration (up to $\dot{\gamma}t=1000$) were conducted with sheets of varying bending rigidities and initial angles $\phi$. Several different dynamical behaviors were observed, as shown in Figure 6. 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 $\absolutevalue{\bar{\mathbf{n}}(t)\dotproduct\mathbf{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 $\phi\neq 0^{\circ}$ 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 $\dot{\gamma}t=1000$ were labeled in Figure 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\times 10^{-4}$ at $\phi=53^{\circ}$ in Figure 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 Figure 6 for large values of $S$ and values of $\phi$ near $90^{\circ}$. 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 $\dot{\gamma}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 $\phi=90^{\circ}$ around a value of $S$ that corresponds to the first predicted buckling transition at $\phi=0^{\circ}$ 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, 19]. Figure 7 shows summary statistics for representative sheets of various values of $S$ with similar initial angles $\phi$. Namely, the eigenvalues, $\lambda_{1}$ and $\lambda_{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(\bar{\mathbf{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. Figure 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 $\phi=53^{\circ}$ and $\phi=54^{\circ}$ in Figure 7. Additionally, the isolated spike during the initially tumbling trajectory of the stiffest sheet (red curve) at $\phi=53^{\circ}$ in Figure 7 represents a highly anisotropic transient crease that occurs as the sheet flips over itself before settling into the flat state.

Several conclusions can be gleaned from the data presented in Figure 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 ($\phi=4^{\circ}$) 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, $\lambda_{1}-\lambda_{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 $\phi=53^{\circ}$ and $\phi=54^{\circ}$ — 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 7). In fact, these continuously tumbling trajectories diverged away from each other exponentially quickly. This, of course, is the classical signature of chaos (see Figure 11 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 Figure 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.

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=\kappa/(\pi\eta\dot{\gamma}L^{3})$. Our findings on the behavior of athermal sheets can be summarized succinctly as follows:

1. 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. 2.

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

3. 3.

   Chaotic trajectories are associated with crumpled conformations.

4. 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.