Toward Nanomechanical Models of Liquid-Phase Exfoliation of Layered 2D Nanomaterials: Analysis of a π − peel Model

Introduction

Graphene and other 2D nanomaterials in the form of atomically thin nanosheets promise unexpected performance in many applications. The nanosheets can be embedded in nanocomposites to make them conductive (Santagiuliana et al., 2018), improve barrier properties or increasing strength and toughness (Rafiee et al., 2010). Or they can be suspended in liquid solvents to produce conductive inks for printed electronics and high-performance coatings (Torrisi and Carey, 2018). While these applications are currently tested in small scale applications, to reach true market impact it is paramount to produce large quantities of 2D nanosheets cheaply, and with control over thickness, lateral area and amount of defects (Ferrari et al., 2015).

A very promising technique for the large-scale production of 2D nanosheets is liquid-phase exfoliation (Hernandez et al., 2008; Coleman, 2009; Yi and Shen, 2015). This technique is relatively simple. It consists in subjecting microparticles of layered 2D nanomaterial (each microparticle being composed of hundreds or thousands of layers) to large mechanical forces that detach the layers. Several different variants of the technique exist, the main ones being turbulent exfoliation (e.g., in rotor-stator mixers, Paton et al., 2014 or microfluidisation devices, Paton et al., 2017), wet ball milling (Knieke et al., 2010; Zhao et al., 2010), and sonication (Alaferdov et al., 2014). These techniques have in common the fact that the microparticles are initially suspended in a liquid (in wet ball milling the particles are wet by a liquid, but a thin layer of liquid is still present). In addition to transmit mechanical stresses, the liquid enables to reduce the adhesion between the layers, prevent reaggregation, and make the mechanical action less aggressive (Shen et al., 2015). Bottom-up synthesis methods, such as Chemical Vapor Deposition in its different variants (Aïssa et al., 2015), are promising for producing high-quality 2D nanomaterials particularly suitable for devices. However, to produce 2D nanomaterials very cheaply for large-scale applications such as nanocomposites, inks or coatings, liquid-phase exfoliation is difficult to beat.

Optimizing liquid-exfoliation processes requires addressing a delicate balance: the mechanical stresses applied to the particles by the fluid must be sufficiently high to delaminate the particles, but not much higher. Excessive stresses can fragment the nanosheets, producing small area sheets of low intrinsic value, or damage the sheets (Johnson et al., 2015). Reaching the right stress level is thus paramount. However, particle-level stresses cannot be controlled directly. What can be controlled are large scale flow variables, such as the mixing power or, equivalently, the average shear rate the suspension is subject to. These are the macroscopic variables that can be controlled by the user in the production process. It would be highly valuable if analytical formulas relating these macroscopic variables to microscopic exfoliation thresholds and time scales were available. Developing these formulas requires an understanding of the flow physics and deformation mechanics at the particle level (Figure 1).

FIGURE 1

Figure 1. Optimizing liquid-phase exfoliation processes (left, reproduced with permission from Paton et al., 2014) requires models to link large-scale flow variables to the micromechanics of exfoliation.

Quite surprisingly, despite the growing importance of liquid-phase processing in the production of graphene and other 2D nanomaterials, the development of theoretical models for liquid-phase exfoliation is at its infancy. Two theoretical models have appeared recently which seem to be relevant. An exfoliation model based on a sliding deformation was proposed by Chen et al. (2012), and later extended by Paton et al. (2014). In this model, the shear forces exerted by the fluid are assumed to balance the rate at which the total surface energy (including solid-solid, solid-liquid, and liquid-liquid interaction energy) changes with respect to the sliding distance. The model was used to describe the dependence of the critical shear rate for exfoliation on the size of the suspended plate-like particles, the viscosity of the fluid, and the energy of adhesion. Chen-Paton's model is insensitive to the mechanical property of the particle: the bending rigidity or young Modulus of the particles do not appear in the model. The assumptions of this model were not stated with sufficient clarity to rigorously assess its validity from a nanomechanical perspective. Prior to the work of Chen, in 2009, a model of exfoliation was developed by Borse and Kamal (2009) to model the exfoliation of multilayer clay particles in polymers. Unlike Chen/Paton's model, Borse's model is sensitive to the mechanical properties of the particle. The formulation of Borse's model is based on Kendall's theory for peeling of elastomers (Kendall, 1975). Borse correctly identified that because the area of contact between the sheets is large, a simple balance between adhesive forces and shear forces gives values of the shear stress too high to be realistic. As a consequence, one must hypothesize that the debonding of the layers is due to stress concentration at the microscopic crack tip of a pre-existing flaw in the inter-layer interface.

Borse's model includes in the fracture mechanics formulation the work done by the fluid forces, the stretching energy associated to the extension of each layer, and the adhesion energy associated to the van der Waals forces between the layers. For the case of constant edge load (identical to the one analyzed by Kendall), Borse and collaborators analyzed the values of the critical shear stress as a function of the peeling angle and the Young modulus of the exfoliated sheet.

We adopt a view similar to that of Borse, and consider exfoliation due to extension of an initial flaw in the inter-layer interface. The extension is caused by a peeling process, whereby the forces driving peeling are hydrodynamic in nature. We envision that the fluid opens a pre-existing crack of length a. The opening angle will depend on the ratio of hydrodynamic and bending forces. If this ratio is small, the opening angle will be small. If this ratio is much larger than one, the flap will turn by 180, and the direction of pulling will be parallel to the direction of propagation of the crack, as illustrated in Figure 2b. This configuration is analogous to the one considered in the “π − peel” mechanical test to measure adhesion (Lin et al., 2002). For brevity, in the current paper, we refer to this configuration as “π-peel” configuration.

FIGURE 2

Figure 2. (a) Peeling deformation of Boron Nitride following wet ball milling; Scanning Electron Microscopy image reproduced with permission from Li et al. (2011). (b) Schematic of the “π-peel” configuration. We assume that the inter-layer interface is debonded over a length a − 2R ≈ a, where R is the radius of curvature of the fold. The local shear flow, whose linear profile is illustrated in the sketch, produces a tangential stress on the flap of order $\mu \dot{\gamma }$.

In this paper we analyze in detail a mathematical model of exfoliation for this “π-peel” configuration, by rigorously justifying the fluid and solid mechanics aspects of the problem. Particularly, we discuss the parameter values for which this configuration may be observed. In contrast to the model by Borse et al., we assume that each layer is perfectly inextensible. Our work takes inspiration from the work on the “inextensible fabric” approximation discussed in Roman (2013), but we recast our results in the context of 2D nanomaterials processing, so that a direct comparison with experimental data from the literature can be made.

The problem we are tackling is one of the first explorations in terms of mechanics of a very complex fluid-structure interaction problem, which needs to be analyzed from different perspectives to be fully understood. Rather than analyzing in depth a sophisticated model, our aim with the paper is to test the prediction of simple models which will enable us to identify the theoretical directions that could bring us close, in terms of orders of magnitude, to the published experimental data for the critical shear rates.

Dimensional Analysis of the Exfoliation Problem

Before considering a particular exfoliation configuration, let us first analyze the general problem of exfoliation from the point of view of dimensional analysis. The critical shear rate for exfoliation, $\dot{\gamma }$, depends on the properties of the fluid (density ρ, viscosity μ, and wetting properties), the mechanical and geometrical properties of the layered micro-particle (particle length L, particle width W, Young modulus of each layer E, bending rigidity of each layer D₀, total number of layers N, flaw size a), and the mechanical and geometrical properties of each inter-layer interface (adhesion energy Γ, area of contact, presence of flaws, etc.).

In exfoliation problems, the particle Reynolds number is typically small, $Re=\rho \dot{\gamma }{L}^2/\mu \ll 1$, so ρ is not an important parameter. Among the mechanical properties, the adhesive properties of each inter-layer, parameterized by Γ and the bending rigidity of each layer D₀ are likely to be dominant controlling parameters. The Young modulus of graphene is huge (E ~ 1TPa, Lee et al., 2008) while its bending rigidity is low (D₀~7eV, Lindahl et al., 2012), so we may regard graphene as an inextensible membrane of finite bending rigidity.

The above parameters suggest a functional relationship of the form

$\begin{array}{ll}\dot{\gamma }=f(\Gamma ,D_0,N,n,\mu ,a,L,W)& (1)\end{array}$

Some simplifications are possible, under the following hypotheses and observations:

• In the zero Reynolds number limit appropriate for colloidal particles, the pressure and viscous stresses scale as $\mu \dot{\gamma }$, so the viscosity enters into the problem only multiplied by the shear rate;

• It can be assumed that the forces induced by the fluid scale proportionally to W. This is not strictly true unless W ≫ L, while typically W ~ L. However, if edge effects are neglected this is a reasonable approximation if the goal is to obtain order of magnitude estimates;

• The direct dependence on n can be neglected if one assumes $\frac{n}{N}\ll 1$. The indirect dependence of the problem on n is still present, because the bending rigidity of the flap depends on n. The bending rigidity of multilayer graphene scales roughly as n³ for n ≥ 3, so for multilayers we can write $D(n)\cong {D_{0}n}^3$, where D₀ ≅ 20eV is the extrapolated value of the single graphene sheet (Sen et al., 2010);

• The dependence on L can be ignored, at least as a first approximation, because the applied hydrodynamic forces on the flap and the “resistive” forces due to elasticity and adhesion depend primarily on a.

Under the hypotheses above, using a and D to make the other parameters dimensionless, Equation (1) can be written as

$\begin{array}{ll}\frac{\dot{\gamma }\mu a^3}{D}=f_1\left(\frac{\Gamma a^2}{D}\right)& (2)\end{array}$

where f₁ is a non-dimensional function of its argument. Equation (2) shows that the non-dimensional critical shear rate, $\frac{\dot{\gamma }\mu a^3}{D}$, is a unique function of the non-dimensional adhesion energy, $\frac{\Gamma a^2}{D}$.

To discuss the comparison of analytical and experimental results, it is useful to assume for f₁ a power-law relationship, for which Equation (2) becomes

$\begin{array}{ll}\frac{\dot{\gamma }\mu a^3}{D}~{\left(\frac{\Gamma a^2}{D}\right)}^{\xi }& (3)\end{array}$

The values the exponent ξ can attain are constrained by physical consideration. It is expected that an increase in flaw size will make the material weaker, so $\dot{\gamma }$ must decrease with increasing a. This requires ξ < 3/2. In addition, an increase in adhesion energy must translate into a larger fluid shear stress. This is only possible if ξ > 0.

An interesting possibility is the case ξ = 1. For this choice of the exponent, the bending rigidity becomes an irrelevant parameter. A complete independence on the bending rigidity can only be plausible if the bending rigidity is so low, that its precise value does not matter (or nearly so).

Exfoliation in the “π-peel” Configuration

Figure 2a, from Li et al. (2011), shows the surface of a nanomaterials after wet ball milling. Flexible layers of nanomaterials have been partially peeled off due to strong shear forces, leaving a fold of very small radius of curvature. A model for this situation can be developed by considering a continuum sheet partially detached from a rigid “mother particle” (Figure 2b). We assume that all the sheets have the same length L. The total thickness of the microparticle is h ≪ L. The geometry of the peeled flap is composed of a curved fold, of radius of curvature R, and a flat portion (from point A to point B in the schematic). The flat portion is subject to a tangential shear stress $\mu \dot{\gamma }$.

In the case of ball milling, the shear stress is created by a relative velocity U between the milling balls acting on a small gap d between the balls, leading to a shear rate $\dot{\gamma }\cong U/d$. In the case of a multilayer particle suspended in a turbulent flow, the ambient shear rate instantaneously “seen” by the particle is the result of the work done by the largest turbulence structures on the smallest, dissipative Kolmogorov eddies (Tennekes and Lumley, 1972). For microparticles smaller than the Kolmogorov scale, equilibrium between the energy rate input P (e.g., the mixing power) and the viscous dissipation occurring at the scale of the particles gives $\dot{\gamma }\cong \sqrt{P/(V\mu )}$, where V is the liquid volume (see, e.g., Varrla et al., 2014 for an application to graphene).

The length of the detached layer is a. Because R ≪ a, the flap length is also approximately equal to a. The tangential force per unit width acting on the flap is $\mu \dot{\gamma }a$.

We need to better specify certain fluid mechanics assumptions. First of all, the model assumes that the particle is aligned with the flow. In fact, plate-like particles rotate when suspended in a shear flow. However, the rotation is very slow when the plate-like particle is nearly aligned with the flow and the aspect ratio is large, with a rate of rotation of the order of $h/L\dot{\gamma }$ (Jeffery, 1922). For most of its rotation period, a plate-like particle of large aspect ratio such as a multilayer microparticle can be considered to be aligned with the flow. A second aspect not included in the model above is the effect of normal hydrodynamic stresses. For a plate-like particle aligned with the flow, normal stresses of the order of $\left(\frac{h}{L}\right)\mu \dot{\gamma }$ act on the surfaces of the particle parallel to the flow direction (e.g., on the flap surface between A and B) (Singh et al., 2014). Because $\frac{h}{L}\ll 1$, along the surface of the flap these stresses are subdominant with respect to the tangential stresses and are thus neglected in the current model. The normal stresses acting on the curved fold, in the direction of the flow, are instead of $O(\mu \dot{\gamma })$. The effect of normal stresses on the model predictions will be considered in section Normal Load on the Curved Fold. We assume that L is smaller than any scale of the flow, so that the flow field around the particle is smooth. For particles in turbulence, this requirement translates to L ≪ η_K, where η_K is the Kolmogorov scale (Landau and Lifshitz, 1986).

The critical value of $\dot{\gamma }$ for exfoliation can be calculated by considering the instantaneous equilibrium between the external work, the change in adhesion energy and the change in bending energy for an inextensible sheet. The analysis is similar to that in Roman (2013), but now the external force is not constant. The velocity of point A (or B) is twice the velocity of the advancing peeling front. As a consequence, as the peeling front moves by an amount da the work done by the external force is $(W\mu \dot{\gamma }a)2da$, where W is the width. The change in adhesion energy is WΓda, where Γ is the adhesion energy per unit area (also called work of separation). Calling $dE(\dot{\gamma })$ the infinitesimal change in adhesion energy corresponding to da, the critical value of $\dot{\gamma }$ for crack initiation satisfies:

$\begin{array}{ll}dE(\dot{\gamma })+W\Gamma da=2W\mu \dot{\gamma }ada& (4)\end{array}$

The bending energy is denoted $E(\dot{\gamma })$ to highlight that this quantity is a function of the shear rate. Considering a curvilinear coordinate s with origin at the crack tip, the bending energy is $E=\frac{1}{2}DW{{\int }_0^a\left(\dot{\theta }\right)}^{2}ds$, where θ(s) is the rotation angle at s and $\dot{\theta }=\frac{d\theta }{ds}$ is the local curvature. The bending energy integral is dominated by the curvature 1/R at the crack tip, which has an extent of the order of 2R. Thus, the order of magnitude of the bending energy is $W\frac{D}{R}$. The function θ(s) can be calculated by considering the equation for the Elastica (Audoly and Pomeau, 2010) for a fixed value of a, which in our case reads

$\begin{array}{ll}\ddot{\theta }+\frac{\mu \dot{\gamma }a}{D}sin(\theta )=0& (5)\end{array}$

Here D is the bending rigidity of the elastic element. Equation (5) is obtained from the equation of equilibrium to rotation for an infinitesimal element of flap, $\frac{dM}{ds}=(F(s)\times t)·e_z$ (Landau and Lifshitz, 1986). In this expression, $M=D\dot{\theta }$ is the moment of the internal stresses, F is the internal force per unit width, t is the local tangent vector, and e_z is the unit vector oriented along the width direction. We have assumed that the only external (hydrodynamic) forces act from point A to point B. Hence, in the curved portion of the flap, the equation of equilibrium to translation $\frac{d\text{F}}{ds}=0$ requires that F is a constant vector (Landau and Lifshitz, 1986). Such constant vector can be easily calculated by noting that point A the force per unit width F is equal to the tension $\cong \mu \dot{\gamma }a{\text{e}}_x$ (the unit vector e_x being parallel to the flow and pointing in the flow direction). As a consequence, $F(s)\cong \mu \dot{\gamma }a{e}_x$ and $(F(s)\times t)·e_z\cong \mu \dot{\gamma }asin\theta$. Inserting this last expression into the equation of equilibrium to rotation gives (5).

The parameter $\frac{\mu \dot{\gamma }a}{D}$ appearing in Equation (5) has the dimensions of the square of a reciprocal length. Because there are no other characteristic lengths in the problem, we anticipate that the curvature of the fold R scales as ${\left(\frac{D}{\mu \dot{\gamma }a}\right)}^{1/2}$.

A solvable first-order non-linear equation can be obtained by multiplying Equation (5) by $\dot{\theta }$, and integrating with respect to s. The result is $\frac{1}{2}{(\dot{\theta })}^2-\frac{\mu \dot{\gamma }a}{D}cos(\theta )=c_1$, where c₁ is a constant. For values of s corresponding to the region from point A to point B, the rotation angle is constant and equal to θ = π. Hence, $c_1=\frac{\mu \dot{\gamma }a}{D}.$ Using this value and the trigonometric identity (1 + cosθ) = 2cos²θ/2, we obtain $\dot{\theta }=2\sqrt{\frac{\mu \dot{\gamma }a}{D}}cos\left(\frac{\theta }{2}\right)$. This is a separable equation whose solution is θ(s′) = 2arcsin(tanh s′), where $s^{\prime }=s/\sqrt{\frac{D}{\mu \dot{\gamma }a}}$. The corresponding bending energy is

$\begin{array}{ll}E=2W\sqrt{\mu \dot{\gamma }aD}& (6)\end{array}$

For $s\ll \sqrt{\frac{D}{\mu \dot{\gamma }a}}$ we have θ(s) ≅ 2s/$\sqrt{\frac{D}{\mu \dot{\gamma }a}}$. Hence, the radius of curvature near the crack tip is of the order of ${\left(\frac{D}{\mu \dot{\gamma }a}\right)}^{1/2}$, as anticipated.

Inserting expression (6) into the energy balance (4) gives

$\begin{array}{ll}{a}^{-1/2}\sqrt{\mu \dot{\gamma }D}+\Gamma =2\mu \dot{\gamma }a& (7)\end{array}$

This expression yields the critical shear rate as a function of the parameters a, D and μ.

It is convenient to recast Equation (7) in terms of the non-dimensional shear rate $\frac{\dot{\gamma }\mu a^3}{D}$and non-dimensional adhesion energy $\frac{\Gamma a^2}{2D}$ (introduced in section Dimensional Analysis of the Exfoliation Problem):

$\begin{array}{ll}\frac{\Gamma a^2}{2D}=\frac{\dot{\gamma }\mu a^3}{D}-\frac{1}{2}\sqrt{\frac{\dot{\gamma }\mu a^3}{D}}& (8)\end{array}$

In contrast to the constant edge force case (Roman, 2013), in our case the bending energy depends on a and the term dE is not zero, giving rise to the square root term on the right-hand side of Equation (8).

From Equation (8), and comparing with Equation (3), the limit $\frac{\mu \dot{\gamma }{a}^3}{D}\to \infty$ gives a power-law exponent ξ = 1, exactly. The effect of the bending energy term is to increase the critical exfoliation value obtained in this asymptotic limit by an amount that depends on the square-root of the shear rate.

Conditions for “π-peel”

For the analytical solution (8) to be valid, the flap length must be long in comparison to the radius of curvature of the fold (Roman, 2013). The condition R ≪ a gives $\frac{\mu \dot{\gamma }{a}^3}{D}\gg 1$. Equation (9) also requires the flap geometry at equilibrium to assume a shape similar to that in Figure 2. For the flap to bend to such an extent the ratio of viscous forces ($~\mu \dot{\gamma }Wa)$to bending forces (~DW/a²) must be large. This, again, gives $\frac{\mu \dot{\gamma }{a}^3}{D}\gg 1$. Hence, the condition $\frac{\mu \dot{\gamma }{a}^3}{D}\gg 1$ simultaneously ensures a good separation of scale between the fold and the flap, and that the assumption of nearly tangential viscous force on the flap is valid.

Figure 3 illustrates how the radius of curvature of the fold, estimated as $R\cong \frac{1}{2}{\left(\frac{D}{\mu \dot{\gamma }a}\right)}^{1/2},$ varies with the length of the flap, for three typical shear rates. Rather than plotting R we plot 2R, which gives a measure of the maximum height of the folded region. The bending rigidity is set to D = 10⁻¹⁸J, close to that of single-layer graphene (Lindahl et al., 2012). In Figure 3A, a dynamic viscosity μ = 0.001 Pa · s is assumed, typical of aqueous solvents and NMP (Paton et al., 2014). In Figure 3B, the viscosity is increased to μ = 1 Pa · s. For μ = 0.001 Pa · s, a good separation of scales between R and a occurs if the shear rate is at least $\dot{\gamma }=10^6{s}^{-1}$. As the viscosity increases to μ = 1 Pa · s, the “π-peel” configuration regime occurs for smaller values of shear rate (Figure 3B).

FIGURE 3

Figure 3. Two times the radius of curvature of the fold vs. length of the flap for (A) μ = 0.001 Pa · s and (B) μ = 1 Pa · s. The bending ridigity is assumed to be D = 10⁻¹⁸J, corresponding to ~6eV. The quantity 2R corresponds approximately to the maximum height of the fold.

The radius of curvature of the flap increases strongly with the number of layers n. Assuming a cubic relation $D=D_0{n}^3$ yields the results shown in Figure 4. For this figure, a = 200 nm and $\dot{\gamma }=10^6{s}^{-1}$; results are shown for three different viscosities. Even considering relatively few layers, for R to be sufficiently small in comparison to a, either the viscosity must be relatively large or, if the viscosity is comparable to that of water, the shear rate must be larger than 10⁶s⁻¹.

FIGURE 4

Figure 4. Two times the radius of curvature of the fold vs. number of layers for a = 200 nm, $\dot{\gamma }=10^6{s}^{-1}$and D = 10⁻¹⁸J.

We can call strongly-stressed nanosheets, 2D material nanosheets for which the scale separation between R and a is complete, and mildly-stressed nanosheets, nanosheets for which R is larger than a, but not by a very large factor. Correspondingly we have two approximations for the critical shear rate. For strongly stressed nanosheets the bending rigidity contribution is negligible and

$\begin{array}{ll}\dot{\gamma }\approx \frac{\Gamma }{2\text{μa}}& (9)\end{array}$

This equation is, strictly speaking, a good approximation when $\frac{\mu \dot{\gamma }{\text{a}}^3}{\text{D}}\to \infty$ (which is equivalent to Γa²/D → ∞). For mildly stressed sheets, solving Equation (8) for $\dot{\gamma }$ yields

$\begin{array}{ll}\frac{\dot{\gamma }\mu {\text{a}}^3}{\text{D}}=\frac{\Gamma a^2}{2D}\left(1+\frac{1}{4(\Gamma a^2/D)}+\frac{\sqrt{1+\frac{8\Gamma a^2}{D}}}{\frac{4\Gamma a^2}{D}}\right)& (10)\end{array}$

This approximation is more accurate than (9) when $\frac{\Gamma a^2}{D}$ (or equivantly $\frac{\mu \dot{\gamma }{\text{a}}^3}{\text{D}})$ is large but finite. For order of magnitude estimates, both (9) and (10) are acceptable. Note that Γ is the work required to separate two surfaces, so Γ/2 is equal to the surface tension of the solid.

Equation (9) is consistent with Kendall's theory for peeling of an extensible thin sheet [Borse's model for exfoliation of clays (Borse and Kamal, 2009) is based on Kendall's theory, so it reaches the same conclusions]. In the case of a constant edge load F applied at an angle ϕ, Kendall calculated that the peeling force for an elastic film of thickness d and Young's modulus E satisfies ${\left(\frac{F}{W}\right)}^2\frac{1}{2dE}+\left(\frac{F}{W}\right)\left(1-cos\varphi \right)=\Gamma$. For E → ∞ (inextensible sheet) and ϕ = π, this equation recovers Equation (9) when F = μaW. Incidentally, the order of magnitude scaling suggested by Equation 9 is identical to that obtained by Paton et al. (2014) for sliding of parallel, rigid platelets, although their configuration is different. In Paton's formulation, the energy term is given by the change in surface energy corresponding to a change in overlap length x between the platelets. Paton's argument for obtaining the critical sliding force is that the overlap area is Wx, and the corresponding surface energy is ΓWx (apart from constant energy terms that do not affect the force). Thus, the sliding force d/dx(Γx) per unit width is constant, and equal to Γ. This results seems to contradict Kendall's theory, which shows that for ϕ = 0 (as in a lap shear joint configuration) F is proportional to $\sqrt{\Gamma }$, not Γ. A possible explanation for the discrepancy is the neglect in Paton's model of the straining of molecular bonds as two rigid sheets slide by a distance comparable to the crystal lattice (Xu et al., 2011) and the disregard of non-uniformities in the interfacial stress (Pugno, 2010).

Normal Load on the Curved Fold

Models (9) and (10) neglect the normal force on the curved fold acting in the direction of the flow. This distributed force pushes the fold in the same direction as the tangential force pulling the flap, so we expect a reduction in the critical shear rate. We can better quantify this statement. The normal force on the flap f can be estimated to be of the order of $\mu \dot{\gamma }RW$. Without attempting to solve Equation (5) exactly by including a distributed force, we can approximate the problem by assuming that f is concentrated in the mid-point of the fold, s ≅ R/2. At this location the internal force F(s) will have a discontinuity: for s > R/2, $F\approx \mu \dot{\gamma }a{e}_x$ as before; for s < R/2, F will increase to $\approx \mu \dot{\gamma }(a+R)e_x$.

The equation for the flap shape for s ≤ R/2 is thus identical to Equation (5), except that in this case (a + R) replaces a. In terms of the non-dimensional curvilinear coordinate $S=s\sqrt{\mu \dot{\gamma }a/D}$, the system to be solved is thus:

$\begin{array}{ll}\frac{d^2\theta }{d{S}^2}+\left(1+\frac{R}{a}\right)sin(\theta )=0forS\le \frac{1}{2}& (\mathrm{11a})\end{array}$

$\begin{array}{ll}\frac{d^2\theta }{d{S}^2}+sin(\theta )=0forS>\frac{1}{2}& (\mathrm{11b})\end{array}$

Figure 5 illustrates a numerical solution of the system above, comparing the rotation angle θ(S) obtained for f = 0 with the corresponding value obtained for $f=\mu \dot{\gamma }RW.$A relatively large ratio $\frac{R}{a}=0.7$is chosen to make the effect of f clearer on the graph. The deviation in θ(S) due to f is small, and decays as S increases. Because of this small deviation, it is useful to decompose θ as $\theta ={\theta }_0+{\theta }^{\prime }$ where θ₀ is the unperturbed solution. Figure 6A illustrates how θ′ varies with S for $\frac{R}{a}=0.5$. The perturbation has a maximum near S = 1, and decays to negligible values for S ≅ 6. Figure 6B shows how the maximum value of θ′ changes with R/a. This figure shows that θ′ ∝ R/a when R/a is sufficiently small. The numerical data suggests

FIGURE 5

Figure 5. Effect of a point load of magnitude $f=\mu R\dot{\gamma }W$ acting on a fold of radius R for R/a = 0.7.

FIGURE 6

Figure 6. (A) Perturbation in the local angle due to the point load for R/a = 0.5. (B) Maximum value of |θ′(S)| vs. R/a'. The dashed line is |θ′(S)| = 0.07R/a.

$\begin{array}{ll}{\theta }^{\prime }(S)=\frac{kR}{a}g(S)& (12)\end{array}$

where k ≈ 0.07 and g is a non-dimensional function whose maximum is 1. Up to first order in θ′ the bending energy can be written as

$\begin{array}{ll}E\cong \frac{1}{2}D\frac{W}{R}{\displaystyle {\int }_0^{\infty }}{\left(\frac{d{\theta }_0}{dS}\right)}^{2}dS+D\frac{W}{R}{\displaystyle {\int }_0^{\infty }}\left(\frac{d{\theta }_0}{dS}\right)\left(\frac{d{\theta }^{\prime }}{dS}\right)dS& (13)\end{array}$

The first integral is the bending energy contribution appearing in Equation 8. The second integral is, to leading order, the change in E due to f. By observing that that the integral from 0 to ∞ converges for values of S of order 1, the order of magnitude of this second term can be estimated as

$\begin{array}{ll}D\frac{W}{R}{\displaystyle {\int }_0^{\infty }}\left(\frac{d{\theta }_0}{dS}\right)\left(\frac{d{\theta }^{\prime }}{dS}\right)dS~D\frac{W}{R}k\left(\frac{R}{a}\right)& (14)\end{array}$

This expression shows that the bending energy contribution originating from the normal load on the fold is, to leading order, independent of R and $\dot{\gamma }$.

Using Equation (14), Equation (7) can be written as

$\begin{array}{ll}{a}^{-1/2}\sqrt{\mu \dot{\gamma }D}+\left(\Gamma -c_{L}k\frac{D}{a^2}\right)=2\mu \dot{\gamma }a& (15)\end{array}$

where the constant c_L is a numerical prefactor. By comparing Equations (15) and (7), we can see that the leading-order effect of f is to reduce the critical shear rate. The net effect is analogous to reducing Γ, by an amount that depends on D and a.

The second-order term neglected in Equation (13), together with the corresponding external work done by f, would give rise to a bending energy contribution scaling as DWk²R/a. This correction would translate to a term $O\left(\frac{k^2}{\sqrt{\frac{\dot{\gamma }\mu a^3}{D}}}\right)$ in Equation 10, much smaller than the retained bending energy terms.

We have initially assumed that the normal force on the fold is $\mu \dot{\gamma }RW$. The correct prefactor, and thus c_L, will depend on the detailed fluid dynamics of the problem. The drag force per unit width on a cylinder of radius R attached to a wall in shear flow is $4\pi \mu \dot{\gamma }R$ (Davis and O'Neill, 1977). If the drag on the fold is assumed to be half that on a cylinder, then $f=2\pi \mu \dot{\gamma }RW$. Based on this estimate, the numerical prefactors should be 2π larger than assumed so far.

The effect of the normal force on the fold is typically small, because R ≪ a. However, it could become important if the tangential force on the flap is reduced. Carbon 2D nanomaterials in contact with water exhibit relatively large slip velocities [with slip lengths in the range 10–80 nm (Thomas and McGaughey, 2008; Falk et al., 2010)], meaning that the no-slip condition at the solid-liquid interface is not satisfied identically. If a large velocity slip is present, the tangential stress on the flap will be reduced from the value $\dot{\gamma }\mu$ assumed in our model. On the other hand, the value of f–being related to a normal force—should not depend strongly on the slip length (a finite normal force on the fold is expected even for infinite slip lengths). In this situation, the effect of the normal force on the fold could become significant.

Modeling the force using a constant distributed force in the horizontal direction, rather than a point force, does not change the order of magnitude of the estimate presented in the current section. Indeed, in this case $\frac{dF}{ds}=-\mu \dot{\gamma }{e}_x$, which integrated gives $F=\mu \dot{\gamma }(s^{*}+a-s)e_x$, where s^* ≅ 2R is the location where the flap starts becoming horizontal. Inserting into the equation for the moment and normalizing gives $\frac{d^2\theta }{d{S}^2}+\left(1+\frac{s^{*}}{a}-\frac{s}{a}\right)sin(\theta )=0$ for s < s^*. The essential difference with Equation (11a) is the term containing sinθ multiplied by s; this difference makes the equation not easily solvable by direct integration (Rohde, 1953). However, because the term in parenthesis is smaller than $\left(1+\frac{s^{*}}{a}\right)$, and s^* ≅ 2R, the equation is bound from above by Equation (11a), if we choose $f=2\mu \dot{\gamma }RW$ instead of $f=\mu \dot{\gamma }RW$. Differences between the point force and distributed force prediction are thus within the uncertainties, which we have discussed above, in the value of the drag force f.

Effect of the Solvent on Adhesion Strength

The parameter Γ is the energy of adhesion per unit area of contact. In vacuum or in an inert gas, the work to separate the surfaces is twice the surface tension of the solid, Γ = 2γ_so (Lawn, 1993).

It is well-known that the presence of a liquid can reduce adhesion. Johnson et al. (1971) carried out lubricated adhesion experiments with rubber spheres, both in air and in liquids. They found that immersion of the surfaces in water reduced the adhesion between the spheres. When the contact was immersed in a 0.01 molar solution of sodium dodecyl sulfate (SDS) the results closely agreed with the Hertz theory down to the lowest loads measured, indicating that adhesion was practically suppressed. Haidara et al. (1995) studied the adhesion of semi spherical PDMS lenses on flat PDMS surfaces. They also found that the presence of surfactants reduced the size of the contact region. The deformations resulting on contacting small (1–2 mm) semispherical lenses of elastomeric poly-(dimethylsiloxane) (PDMS) with the flat sheets of this material were measured in air and in mixtures of water and methanol by Chaudhury and Whitesides (1991). They found that the adhesion between PDMS surfaces was strongest in water, and decreased as the hydrophobicity of the medium decreased. van Engers et al. (2017) using a method conceptually similar to that of Johnson, Kendall and Roberts measured experimentally the graphene-graphene interfacial energy γ_so, obtaining 115 mN/m in inert gas, 83 in water, 29 in sodium cholate, a surfactant. The value obtained for graphene in an inert gas is close to the value 110 mN/m previously reported for graphene-graphite interaction by Wang et al. (2016). Shih et al. (2010) carried out molecular dynamics of the interaction between rigid graphene sheets in a variety of solvents, and obtained a minimum in the interaction energy of 250kJmol⁻¹nm⁻² for water and around 100kJmol⁻¹nm⁻² for N-methylpyrrolidone (NMP). Because the reduction in cohesive stresses can be related to a reduction in the depth of the potential energy well for the molecular bonds near the crack tip (Stoloff and Johnston, 1963; Lawn, 1993), this last data suggests that the values of the surface tension for water and NMP are of comparable order of magnitude (despite the fact that NMP is considered a much better solvent for graphene than water!).

From the data above there is evidence to suggest that: (i) the intrinsic value of Γ corresponding to vacuum or an inert gas is in the range 0.20–0.25 N/m; (ii) “good” solvents can reduce Γ, but probably not by several orders of magnitude. To compare with experimental data we will assume that Γ in typical liquid solvents ranges from 0.01 N/m to 0.1 N/m, with the lower value being characteristic of a good solvent.

Critical Shear Rate: Comparison With Experimental Data

We have seen that if the microscopic peeling configuration is as in Figure 2 the order of magnitude of the critical shear rate is given by $\dot{\gamma }\approx \frac{\Gamma }{2\text{μa}}$, i.e., Equation 3 with ξ = 1. Figure 7A compares $\dot{\gamma }=\frac{\Gamma }{2\text{μa}}$ curves against experimental data for μ = 10⁻³Pa · s and two surface energy values: Γ = 0.1N/m and Γ = 0.01N/m. Each experimental case is denoted by a horizontal row of symbols corresponding to several values of a. This was done because the value of a in experiments is unknown. Experimental cases correspond to turbulent exfoliation in a rotor-stator mixer (Paton et al., 2014), turbulent exfoliation in microfluidization (Karagiannidis et al., 2017; Paton et al., 2017) and wet ball milling (Knieke et al., 2010). In the case of wet ball milling, a typical stress energy SE~0.134 μJ was reported for ball diameter d_gm = 100 μm and rotation rate 1,500 rpm. The critical stress energy can be converted to a stress τ of the order of 10⁵Pa assuming that SE is dissipated within a contact region of volume ${~d}_{gm}^3$. With $\tau =\mu \dot{\gamma }$ and viscosity μ = 10⁻³Pa · s, a stress of 10⁵Pa corresponds to an equivalent critical shear rate value of 10⁸s⁻¹. This is the estimated value reported in Figure 7A.

FIGURE 7

Figure 7. (A) Critical shear rate for exfoliation as a function of flaw size (in nanometers) for different scaling exponents, assuming D = 10⁻¹⁸J. Experimental references: Knieke et al. (2010), Paton et al. (2014, 2017), and Karagiannidis et al. (2017). (B) In the “cleavage” configuration the fluid stresses act normal to the flap, and the opening angle ϕ is small.

The graph shows that for realistic values of Γ and for a range of plausible values of flaw size (we can estimate that for particles of length L = 1μm the flaw size ranges roughly from 50 nm to 500 nm) the expression $\dot{\gamma }\approx \frac{\Gamma }{2\mu a}$ largely overestimates the experimental data for turbulent exfoliation (i.e., exfoliation in rotor-stator mixers or microfluidisation). Even for Γ = 0.01 N/m, which is smaller than any measured value for the adhesion energy of graphene (see discussion in section Effect of the Solvent on Adhesion Strength), the critical shear rate is much larger than what reported for rotating mixer and microfluidization. Instead, the values predicted by $\dot{\gamma }\approx \frac{\Gamma }{2\text{μa}}$are reasonably close to those estimated for ball milling, a technique that allows large stresses to be produced.

An Alternative Model: Flow-Induced Cleavage (for Small Opening Angles)

From Figure 7A it can be seen that if an exponent ξ = 1/2 is assumed, a much closer agreement with the experimental data for rotating mixer and microfluidisation approaches is achieved. An exponent of 1/2 can be obtained if we assume that the bending of the flap is caused by a normal stress on the flap of order of $\mu \dot{\gamma }$, as in the “cleavage” configuration illustrated in Figure 7B. In this case the flap would be subject to a force ~ μ$\dot{\gamma }aW$ acting on a lever arm ~ a. Equating the corresponding external moment ~ μ$\dot{\gamma }{a}^{2}W$ to the bending moment $DW\sqrt{\Gamma /D}$ required for fracture initiation (Obreimoff, 1930) leads to $\frac{\mu \dot{\gamma }{a}^3}{D}~{(\Gamma a^2/D)}^{1/2}$. The same result can be obtained by considering the analytical solution for the displacement of a cantilever beam of length a subject to a constant normal load $\mu \dot{\gamma }$ (Timošenko, 1940):

$\begin{array}{ll}w=\frac{\mu \dot{\gamma }}{24D}{(x-a)}^2(x^2+2ax+3a^2)& (16)\end{array}$

The bending energy per unit width $U=\frac{D}{2}\int {\left(\frac{d^{2}w}{d{x}^2}\right)}^{2}dx$is quadratic in $\mu \dot{\gamma }$. The strain energy release rate $G=\frac{dU}{da}$, with units of energy per unit area, is proportional to $\frac{{{(\mu \dot{\gamma })}^{2}a}^4}{D}$. Equating the strain energy release rate to Γ (Griffith's criterion, Lawn, 1993) gives again the scaling

$\begin{array}{ll}\frac{\mu \dot{\gamma }{a}^3}{D}~{(\Gamma a^2/D)}^{1/2}& (\mathrm{17a})\end{array}$

An intuitive explanation for why the critical shear rate in the “cleavage” configuration is smaller than in the “π -peel” configuration is that even if the forces applied to the flap are the same in both cases, the lever arm in the “π -peel” configuration is of $O\left(\frac{R}{a}\right)$ smaller than in the “cleavage” configuration. For a given value of a, reaching the critical bending moment in the “cleavage” configuration therefore requires as smaller value of $\dot{\gamma }$.

Equation (17a) is valid provided that the displacement of the flap is small. For larger displacements, one needs to account for two factors, which are extensively discussed in a recent paper by the author (Salussolia, 2019): (i) even in the case of uniform pressure applied to the flap, the direction of load depends on the normal to the flap, which itself depends on the shape of the flap (i.e., the load is follower); (ii) as the flap displacement increases, the pressure on the flap increases, essentially because more area of the flap is exposed to flow. Factor (ii) makes the dependence on Γa²/D weaker, reducing effectively the power-law exponent significantly below 1/2. The inclusion of large-deformation effects does not change this conclusion, although it changes somewhat the value of the critical shear rate (without a change in order of magnitude of this quantity).

Time Scales of Microstructural Rearrangement: Toward Exfoliation Kinetics

Once the critical shear rate is reached, the fracture will propagate and exfoliation will occur. But exfoliation is not an instantaneous event. If this was the case, at a critical shear rate one would obtain complete exfoliation of all the microplates. This is not observed in practice. For instance, Paton and co-workers found that the concentration of exfoliated material follows a rather slow kinetics, apparently governed by a power-law of time (Paton et al., 2014). While explaining the emergence of this power-law requires tools that go beyond simple micromechanics, it is instructive to consider examples of dissipation processes occurring in the vicinity of the crack tip that could make the exfoliation mechanics depend on time.

An important dissipative process is due to the viscosity of the solvent. As the crack of the interlayer interface propagates, fluid must be drawn in from the surroundings toward the crack tip. Because the gap distance between the layers is small in the crack region, the viscous dissipation can be substantial (Rieutord et al., 2005; Lister et al., 2013). The rate of dissipation per unit volume of fluid is

$\begin{array}{ll}\epsilon ~\mu {\left(\frac{U(x)}{h(x)}\right)}^2& (\mathrm{17b})\end{array}$

where h(x) is the gap height at a position x, and U(x) is the characteristic fluid velocity at x (Eggers and Fontelos, 2015). The total power per unit width dissipated in the fluid can be estimated to be of the order of (Rieutord et al., 2005).

$\begin{array}{ll}\frac{P_d}{W}={\displaystyle {\int }_{x_{min}}^{x_{max}}}{\displaystyle {\int }_0^h}\epsilon dxdy~\mu \frac{V^2}{H^2}(H{\ell }_D)& (18)\end{array}$

where V is the crack tip speed. Here H and ℓ_D are the characteristic height and length of the wedge-like region in the vicinity of the crack tip where the bulk of the viscous dissipation occurs. The high-dissipation region is limited by the coordinates x_min and x_max. If there is no external force applied to the flap, the power to drive the crack motion against the viscous dissipation is provided by the adhesion force, $\frac{P_d}{W}=\Gamma V$, leading to (Rieutord et al., 2005)

$\begin{array}{ll}V~\frac{\Gamma }{\mu }\frac{H}{{\ell }_D}& (19)\end{array}$

In the “π -peel” case examined in the current paper, an external force per unit width $\mu \dot{\gamma }a$ is applied to the flap and the corresponding driving power is $\frac{P_d}{W}=(2\mu \dot{\gamma }a-\Gamma )V$. Hence, the crack velocity, for the case in which viscous dissipation is the only effect resisting the motion of the crack, can be estimated to be of the order of

$\begin{array}{ll}{V}_{viscous}~\frac{\left(2\mu \dot{\gamma }a-\Gamma \right)}{\mu }\frac{H}{{\ell }_D}& (20)\end{array}$

The inverse of the geometrical ratio $\frac{H}{{\ell }_D}$ is equal to the integral of h(x) from x_min to x_max (Rieutord et al., 2005). Because the region near the crack tip is thin and slender, we have $\frac{H}{{\ell }_D}\ll 1$. Thus, $V_{viscous}\ll \frac{(2\mu \dot{\gamma }a-\Gamma )}{\mu }$.

A second important dissipative process is caused by the time-dependent rupture of the molecular bonds in the adhesion zone. A model for such process assumes that the rupture rate is governed by Maxwell-Boltzmann statistics. For small deviations from Griffith's condition, this model gives an exponential dependence of the crack velocity on the difference between the strain energy release rate G and the energy of adhesion per unit area of contact 2γ_so (Lawn, 1993):

$\begin{array}{ll}V={2v}_0{a}_{0}exp\left(-\frac{U^{\circ }}{kT}\right)sinh\left(\alpha \frac{G-2{\gamma }_{so}}{kT}\right)& (21)\end{array}$

Here, a₀ is the characteristic lattice dimension, v₀ is a molecular frequency (typically a few THz), α is an activation area, kT is the thermal energy, γ_so is the surface tension of the solid (section Effect of the Solvent on Adhesion Strength) and ΔU° is a quiescent activation energy. In the “π -peel” configuration, the external work $2\mu \dot{\gamma }a$ takes the place of G. Using the definition Γ = 2γ_so the following expression crack velocity expression is obtained:

$\begin{array}{ll}{V}_{br}~V_{br,0}sinh\left(\alpha \frac{2\mu \dot{\gamma }a-\Gamma }{kT}\right)& (22)\end{array}$

where

$\begin{array}{ll}{V}_{br,0}={2v}_0{a}_{0}exp\left(-\frac{\Delta U^{\circ }}{kT}\right)& (23)\end{array}$

The molecular velocity v₀a₀ is of the order of 100 m/s−1 km/s, i.e., sonic velocities. The energy barrier ΔU^∘ is typically taken to be of the order of 10kT (Brochard-Wyart and de Gennes, 2003). Setting ΔU^∘ = 25kT, gives values of V_br,0 in the range 1−10 nm/s. The effect of the solvent (through chemically-assisted bond rupture) can be accounted for through the dependence of Γ on the solvent, and by appropriately modifying the activation terms (Stoloff and Johnston, 1963; Lawn, 1993).

We can use Equations (20) and (22) to make some considerations regarding the relative importance of viscous dissipation and dissipation due to bond rupture. Let us assume that we are carrying out exfoliation using a shear rate of the order of the critical one, a situation that is expected to hold in practice. Because $\frac{2\mu \dot{\gamma }a}{\Gamma }=O(1)$ and $\frac{H}{{\ell }_D}<1$, an upper bound for V_viscous is the “capillary velocity” $\frac{\Gamma }{\mu }$. This quantity is ~100 m/s for solvents having the viscosity of water (using a reference value Γ = 0.1 N/m). The velocity due to (interlayer) molecular bond rupture can be estimated to be, typically, larger than this value. This is because for a₀ in the nanoscopic range the energy scale Γα is typically larger than the energy barrier ΔU° (taking α = 1nm² and Γ = 0.1 N/m, we get Γα ≈ 250kT). As a result, the growing hyperbolic sine term is much larger than the decaying term $exp\left(-\frac{\Delta U^{\dot{}}}{kT}\right)$. Noting that the “capillary velocity” $\frac{\Gamma }{\mu }$ is roughly comparable in magnitude to v₀a₀, we have V_viscous ≪ V_br.

Thus, viscous dissipation can be very important. The relative importance with respect to dissipation due to bond rupture (which we estimate to be subdominant) will depend on the very specific values of α and ΔU°. Recent measurements of self-tearing and peeling of graphene sheets from solid substrates suggest a method to calculate these parameters (Annett and Graham, 2016).

In the regime of viscous dominated dissipation, the characteristic time required to complete peeling of a flap from a particle of lateral size L is

$\begin{array}{ll}{T}_v=\frac{L}{V_{viscous}}~\frac{\mu }{(2\mu \dot{\gamma }L-\Gamma )}\frac{{\ell }_D}{H}& (24)\end{array}$

For values of the applied shear rate much larger than the critical one the viscous peeling time scales proportionally to the inverse shear rate, $T_v~{\dot{\gamma }}^{-1}\frac{{\ell }_D}{2H}$. This time scale depends on the viscosity only through the dependence of the ratio $\frac{{\ell }_D}{2H}$ on hydrodynamics. Because $\frac{{\ell }_D}{H}\gg 1$ we expect $T_v\gg {\dot{\gamma }}^{-1}.$ The viscous peeling time scale is proportional to the local shear rate, with a large numerical prefactor.

In this section, we have focused on the “π -peel” configuration. The arguments proposed can be extended to the “cleavage' configuration by replacing the work of the external force per unit width with the strain energy release rate.

Conclusions: Frontiers in the Hydrodynamics of 2D Nanomaterials

We have analyzed theoretically a simple model of hydrodynamic peeling, by focusing mostly on the “π - peel” configuration (Figure 2). Based on our analysis, the shear rates that one needs to apply to the fluid for this configuration to result in exfoliation are rather large, of the order of 10⁸ − 10⁹s⁻¹. These large shear rates are available in exfoliation methods that produce large stresses on the particle, such as ball milling, but not in most published exfoliation approaches in which the shear rate is produced by conventional turbulence (as in rotor-stator shear mixing or microfluidisation approaches). The fact that in rotor-stator or microfluidisation approaches exfoliation is observed for shear rates in the range 10⁴ − 10⁶s⁻¹ suggests that alternative microscale exfoliation configurations are likely to dominate the exfoliation micromechanics in these approaches, at least in the initial stages of exfoliation.

We suggest that, due to its stronger dependence on the flaw size a, a “cleavage” configuration (Figure 7B) can yield critical shear rate values much smaller than those predicted for “π -peel,” closer to the values observed experimentally. While for “π -peel” the critical shear rate decays as $\frac{1}{a}$(Equation 9), for cleavage under a constant pressure $\mu \dot{\gamma }$ the critical shear rate decays as $\frac{1}{a^2}$ (Equation 17a). We anticipate that if the pressure is not constant, but depends on the configuration of the wedge creating the fracture, the dependence on a can be even stronger. This effect results in even lower values of the critical shear rate. The analysis of this case is the subject of a separate paper (Salussolia, 2019).

This initial study on the hydrodynamics of 2D materials exfoliation only begins to uncover the complexity of the exfoliation process as seen at the scale of each particle. To better understand the micromechanics of liquid-phase exfoliation, future computational investigations must consider the fully coupled fluid-structure interaction problem and include atomistic details. In the wedge region, near the crack tip, the liquid is strongly confined, with gaps of nanometric dimensions. Hence the problem severely challenges the use of continuum approaches. For carbon nanomaterials in contact with water, hydrodynamic slip characterizes the flow behavior at the solid-liquid surface (Tocci et al., 2014; Striolo et al., 2016), with slip lengths of the order of a few tens of nanometers (Thomas and McGaughey, 2008; Falk et al., 2010). Our work provides estimates for geometric quantities—such as the shear rate dependent radius of curvature of the fold R–that are obtained neglecting molecular effects. These quantities can be used as a reference to evaluate the limitations of continuum treatments. Molecular dynamics studies for the configurations investigated in the current paper could provide insights into the general process of liquid intercalation under flow, a topic that could lead to the design of improved solvents for the exfoliation of 2D nanomaterials.

Another important topic is the prediction of the kinetics of exfoliation. We have estimated that the time scale for microscopic peeling can in many cases be controlled by viscous dissipation. For “π -peel” at shear rates largely exceeding the critical one, the time scale for viscous peeling has been found to scale proportionally to ${\dot{\gamma }}^{-1}$ (with a large prefactor, $\frac{{\ell }_D}{H}\gg 1$, see Equation 24). This time scale has to be compared to two dynamic time scale: the time scale of particle rotation (of the order of ${\dot{\gamma }}^{-1}\Lambda$, where Λ is the aspect ratio, Challabotla et al., 2015), and the time scale of permanence in a turbulent structure (which can be estimated to be of the order of ${\dot{\gamma }}^{-1},$ Babler et al., 2012). The inter-play of these different time scales should give rise to a rich phase diagram. Uncovering the features of this diagram, through a comparison of multi-particle simulations and experiments (Voth, 2015), should be of interest for researchers investigating the statistical physics of complex systems, and is a necessary step toward predicting the yield of industrial-scale exfoliation processes.

Theoretical modeling of exfoliation processes in sheared liquids is an important and currently largely unexplored frontier in carbon nanomaterials research. Understanding how graphene “breaks” under the action of fluid dynamic forces has implications not only for large-scale graphene production but also for quantifying the time-evolving size and thickness distribution of graphene (or any other 2D material) during its transport. In industrial liquids such as lubricants and paints, the performance of a graphene additive whose size depends on time—because of fragmentation or exfoliation—will also depend on time. Size and shape are dominant variables affecting how a nanoparticle interacts with a biological cell, and the prediction of these variables during flow is thus essential to evaluate toxicological effects on human health or the environment (Wick et al., 2014). Our work on the modeling of graphene hydrodynamics, funded by the European Research Council, lays basic theoretical building blocks that are necessary to develop the next generation of predictive multiscale software for fragmentation, exfoliation and liquid processing of 2D nanomaterials.

Author Contributions

LB designed the research, carried out the analysis, and wrote the paper.

Conflict of Interest

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

Funding from ERC Starting Grant FLEXNANOFLOW (no. 715475) is gratefully acknowledged. The author would like to thank Dr. E. Barbieri, and G. Salussolia for useful scientific discussions.

References

Aïssa, B., Memon, N. K., Ali, A., and Khraisheh, M. K. (2015). Khraisheh recent progress in the growth and applications of graphene as a smart material: a review. Front. Mater. 2:58. doi: 10.3389/fmats.2015.00058

CrossRef Full Text | Google Scholar

Alaferdov, A. V., Gholamipour-Shirazi, A., Canesqui, M. A., Yu Danilov, A., and Moshkalev, S. A. (2014). Size-controlled synthesis of graphite nanoflakes and multi-layer graphene by liquid phase exfoliation of natural graphite. Carbon 69, 525–535. doi: 10.1016/j.carbon.2013.12.062

CrossRef Full Text | Google Scholar

Annett, J., and Graham, L. W. (2016). Cross Self-assembly of graphene ribbons by spontaneous self-tearing and peeling from a substrate. Nature 535, 271–275. doi: 10.1038/nature18304

CrossRef Full Text | Google Scholar

Audoly, B., and Pomeau, Y. (2010). Elasticity and Geometry: From Hair Curls to the Non-linear Response of Shells. Oxford, UK: Oxford University Press

Google Scholar

Babler, M. U., Biferale, L., and Lanotte, A. S. (2012). Lanotte. Breakup of small aggregates driven by turbulent hydrodynamical stress. Phys. Rev. 85:025301. doi: 10.1103/PhysRevE.85.025301

CrossRef Full Text | Google Scholar

Borse, N. K., and Kamal, M. R. (2009). Estimation of stresses required for exfoliation of clay particles in polymer nanocomposites. Polymer Eng. Sci. 49, 641–650. doi: 10.1002/pen.21211

CrossRef Full Text | Google Scholar

Brochard-Wyart, F., and de Gennes, P.-G. (2003). Unbinding of adhesive vesicles. Comp. Rendus Phys. 4, 281–287. doi: 10.1016/S1631-0705(03)00048-3

CrossRef Full Text | Google Scholar

Challabotla, N. R., Zhao, L., and Andersson, H. I. (2015). Orientation and rotation of inertial disk particles in wall turbulence. J. Fluid Mech. 766:R2. doi: 10.1017/jfm.2015.38

CrossRef Full Text | Google Scholar

Chaudhury, M. K., and Whitesides, G. M. (1991). Direct measurement of interfacial interactions between semispherical lenses and flat sheets of poly (dimethylsiloxane) and their chemical derivatives. Langmuir 7, 1013–1025. doi: 10.1021/la00053a033

CrossRef Full Text | Google Scholar

Chen, X., Dobson, J. F., and Raston, C. L. (2012). Vortex fluidic exfoliation of graphite and boron nitride. Chem. Commun. 48, 3703–3705. doi: 10.1039/c2cc17611d

PubMed Abstract | CrossRef Full Text | Google Scholar

Coleman, J. N. (2009). Liquid-phase exfoliation of nanotubes and graphene. Adv. Funct. Mater. 19, 3680–3695. doi: 10.1002/adfm.200901640

CrossRef Full Text | Google Scholar

Davis, A. M. J., and O'Neill, M. E. (1977). Separation in a slow linear shear flow past a cylinder and a plane. J. Fluid Mech. 81, 551–564. doi: 10.1017/S0022112077002225

CrossRef Full Text | Google Scholar

Eggers, J., and Fontelos, M. A. (2015). Singularities: Formation, Structure, and Propagation. Cambridge: Cambridge University Press.

Google Scholar

Falk, K., Sedlmeier, F., Joly, L., Netz, R. R., and Bocquet, L. (2010). Molecular origin of fast water transport in carbon nanotube membranes: superlubricity versus curvature dependent friction. Nano Lett. 10, 4067–4073. doi: 10.1021/nl1021046

PubMed Abstract | CrossRef Full Text | Google Scholar

Ferrari, A. C., Bonaccorso, F., Fal'ko, V., Novoselov, K. S., Roche, S., Bøggild, P., et al. (2015). Technology roadmap for graphene, related two-dimensional crystals, and hybrid systems. Nanoscale 7, 4598–4810. doi: 10.1039/C4NR01600A

PubMed Abstract | CrossRef Full Text | Google Scholar

Haidara, H., Chaudhury, M. K., and Owen, M. J. (1995). A direct method of studying adsorption of a surfactant at solid-liquid interfaces. J. Phys. Chem. 21, 8681–8683. doi: 10.1021/j100021a037

CrossRef Full Text | Google Scholar

Hernandez, Y., Nicolosi, V., Lotya, M., Blighe, F. M., Sun, Z., De, S., et al. (2008). High-yield production of graphene by liquid-phase exfoliation of graphite. Nat. Nanotechnol. 3:563. doi: 10.1038/nnano.2008.215

PubMed Abstract | CrossRef Full Text | Google Scholar

Jeffery, G. B. (1922). The motion of ellipsoidal particles immersed in a viscous fluid. Proc. R. Soc. Lond. A. 102, 161–179. doi: 10.1098/rspa.1922.0078

CrossRef Full Text | Google Scholar

Johnson, D. W., Dobson, B. P., and Coleman, K. S. (2015). A manufacturing perspective on graphene dispersions. Curr. Opin. Colloid Interface Sci. 20, 367–382. doi: 10.1016/j.cocis.2015.11.004

CrossRef Full Text | Google Scholar

Johnson, K. L., Kendall, K., and Roberts, A. D. (1971). Surface energy and the contact of elastic solids. Proc. R. Soc. Lond. A 20, 301–313. doi: 10.1098/rspa.1971.0141

CrossRef Full Text | Google Scholar

Karagiannidis, P. G., Hodge, S. A., Lombardi, L., Tomarchio, F., Decorde, N., Milana, S., et al. (2017). Microfluidization of graphite and formulation of graphene-based conductive inks. ACS Nano. 11, 2742–2755. doi: 10.1021/acsnano.6b07735

PubMed Abstract | CrossRef Full Text | Google Scholar

Kendall, K. (1975). Thin-film peeling - the elastic term. J. Phys. D Appl. Phys. 8:1449. doi: 10.1088/0022-3727/8/13/005

CrossRef Full Text | Google Scholar

Knieke, C., Berger, A., Voigt, M., Taylor, R. N. K., Röhrl, J., and Peukert, W. (2010). Scalable production of graphene sheets by mechanical delamination. Carbon 48, 3196–3204. doi: 10.1016/j.carbon.2010.05.003

CrossRef Full Text | Google Scholar

Landau, L. D., and Lifshitz, E. M. (1986). Theory of Elasticity. Amsterdam: Elsevier, 7.

Google Scholar

Lawn, B. R. (1993). Fracture of Brittle Solids. Cambridge: Cambridge University Press.

Google Scholar

Lee, C., Wei, X., Kysar, J. W., and Hone, J. (2008). Measurement of the elastic properties and intrinsic strength of monolayer graphene. Science 321, 385–388. doi: 10.1126/science.1157996

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, L. H., Chen, Y., Behan, G., Zhang, H., Petravic, M., and Glushenkov, A. M. (2011). Large-scale mechanical peeling of boron nitride nanosheets by low-energy ball milling. J. Mater. Chem. 21, 11862–11866. doi: 10.1039/c1jm11192b

CrossRef Full Text | Google Scholar

Lin, Y. Y., Hui, C. Y., and Wang, Y. C. (2002). Modeling the failure of an adhesive layer in a peel test. J. Polymer Sci. Part B Polym. Phys. 40, 2277–2291. doi: 10.1002/polb.10289

CrossRef Full Text | Google Scholar

Lindahl, N., Midtvedt, D., Svensson, J., Nerushev, O. A., Lindvall, N., Isacsson, A., et al. (2012). Determination of the bending rigidity of graphene via electrostatic actuation of buckled membranes. Nano Lett. 12, 3526–3531. doi: 10.1021/nl301080v

PubMed Abstract | CrossRef Full Text | Google Scholar

Lister, J. R., Peng, G. G., and Neufeld, J. A. (2013). Viscous control of peeling an elastic sheet by bending and pulling. Phys. Rev. Lett. 111:154501. doi: 10.1103/PhysRevLett.111.154501

PubMed Abstract | CrossRef Full Text | Google Scholar

Obreimoff, J. W. (1930). The splitting strength of mica. Proc. R. Soc. Lond. Ser. A 127:290. doi: 10.1098/rspa.1930.0058

CrossRef Full Text | Google Scholar

Paton, K. R., Anderson, J., Pollard, A. J., and Sainsbury, T. (2017). Production of few-layer graphene by microfluidization. Mater. Res. Exp. 4:025604. doi: 10.1088/2053-1591/aa5b24

CrossRef Full Text | Google Scholar

Paton, K. R., Varrla, E., Backes, C., Smith, R. J., Khan, U., O'Neill, A., Paton, K. R., et al. (2014). Scalable production of large quantities of defect-free few-layer graphene by shear exfoliation in liquids. Nat. Mater. 13:624. doi: 10.1038/nmat3944

PubMed Abstract | CrossRef Full Text | Google Scholar

Pugno, N. M. (2010). The design of self-collapsed super-strong nanotube bundles. J. Mech. Phys. Solids. 58, 1397–1410. doi: 10.1016/j.jmps.2010.05.007

CrossRef Full Text | Google Scholar

Rafiee, M. A., Rafiee, J., Srivastava, I., Wang, Z., Song, H., Yu, Z. Z., et al. (2010). Fracture and fatigue in graphene nanocomposites. Small 6, 179–183. doi: 10.1002/smll.200901480

PubMed Abstract | CrossRef Full Text | Google Scholar

Rieutord, F., Bataillou, B., and Moriceau, H. (2005). Dynamics of a bonding front. Phys. Rev. Lett.. 94:236101. doi: 10.1103/PhysRevLett.94.236101

PubMed Abstract | CrossRef Full Text | Google Scholar

Rohde, F. (1953). Large deflections of a cantilever beam with uniformly distributed load. Q. Appl. Math. 11, 337–338. doi: 10.1090/qam/56438

CrossRef Full Text | Google Scholar

Roman, B. (2013). Fracture path in brittle thin sheets: a unifying review on tearing. Int. J. Fract. 182, 209–237. doi: 10.1007/s10704-013-9869-5

CrossRef Full Text | Google Scholar

Salussolia, G. (2019). Micromechanics of liquid-phase exfoliation of a layered 2D material: a hydrodynamic peeling model. J. Mech. Phys. Solids. 134:103764. doi: 10.1016/j.jmps.2019.103764

CrossRef Full Text | Google Scholar

Santagiuliana, G., Picot, O. T., Crespo, M., Porwal, H., Zhang, H., Li, Y., et al. (2018). Breaking the nanoparticle loading–dispersion dichotomy in polymer nanocomposites with the art of croissant-making. ACS Nano. 12, 9040–9050. doi: 10.1021/acsnano.8b02877

PubMed Abstract | CrossRef Full Text | Google Scholar

Sen, D., Novoselov, K. S., Reis, P. M., and Buehler, M. J. (2010). Tearing graphene sheets from adhesive substrates produces tapered nanoribbons. Small 6, 1108–1116. doi: 10.1002/smll.201000097

PubMed Abstract | CrossRef Full Text | Google Scholar

Shen, J., He, Y., Wu, J., Gao, C., Keyshar, K., Zhang, X., et al. (2015). Liquid phase exfoliation of two-dimensional materials by directly probing and matching surface tension components. Nano Lett. 15, 5449–5454. doi: 10.1021/acs.nanolett.5b01842

PubMed Abstract | CrossRef Full Text | Google Scholar

Shih, C. J., Lin, S., Strano, M. S., and Blankschtein, D. (2010). Understanding the stabilization of liquid-phase-exfoliated graphene in polar solvents: molecular dynamics simulations and kinetic theory of colloid aggregation. J. Am. Chem. Soc. 132, 14638–14648. doi: 10.1021/ja1064284

PubMed Abstract | CrossRef Full Text | Google Scholar

Singh, V., Koch, D. L., Subramanian, G., and Stroock, A. D. (2014). Rotational motion of a thin axisymmetric disk in a low Reynolds number linear flow. Phys. Fluids 26:033303. doi: 10.1063/1.4868520

CrossRef Full Text | Google Scholar

Stoloff, N. S., and Johnston, T. L. (1963). Crack propagation in a liquid metal environment. Acta Metallurg. 11, 251–256. doi: 10.1016/0001-6160(63)90180-9

CrossRef Full Text | Google Scholar

Striolo, A., Michaelides, A., and Joly, L. (2016). The carbon-water interface: modeling challenges and opportunities for the water-energy nexus. Annu. Rev. Chem. Biomol. Eng. 7, 533–556. doi: 10.1146/annurev-chembioeng-080615-034455

PubMed Abstract | CrossRef Full Text | Google Scholar

Tennekes, H., and Lumley, J. L. (1972). A First Course in Turbulence. Massachusetts: MIT press.

Google Scholar

Thomas, J. A., and McGaughey, A. J. H. (2008). Reassessing fast water transport through carbon nanotubes. Nano Lett. 8, 2788–2793. doi: 10.1021/nl8013617

PubMed Abstract | CrossRef Full Text | Google Scholar

Timošenko, S. P. (1940). Strength of Materials: Elementary Theory and Problems. New York, NY: D. Van Nostrand Company.

Google Scholar

Tocci, G., Joly, L., and Michaelides, A. (2014). Friction of water on graphene and hexagonal boron nitride from ab initio methods: very different slippage despite very similar interface structures. Nano Lett. 12, 6872–6877. doi: 10.1021/nl502837d

CrossRef Full Text | Google Scholar

Torrisi, F., and Carey, T. (2018). “Printing 2D materials,” in Flexible Carbon-Based Electronics, eds P. Palermo and V. Samorí (Hoboken, NJ: John Wiley and Sons), 131–184.

Google Scholar

van Engers, C. D., Cousens, N. E., Babenko, V., Britton, J., Zappone, B., Grobert, N., et al. (2017). Direct measurement of the surface energy of graphene. Nano Lett. 6, 3815–3821. doi: 10.1021/acs.nanolett.7b01181

CrossRef Full Text | Google Scholar

Varrla, E., Paton, K. R., Backes, C., Harvey, A., Smith, R. J., McCauley, J., et al. (2014). Turbulence-assisted shear exfoliation of graphene using household detergent and a kitchen blender. Nanoscale 6, 11810–11819. doi: 10.1039/C4NR03560G

PubMed Abstract | CrossRef Full Text | Google Scholar

Voth G. A. (2015). Disks aligned in a turbulent channel. J. Fluid Mech. 772, 1–4. doi: 10.1017/jfm.2015.144

CrossRef Full Text | Google Scholar

Wang, J., Sorescu, D. C., Jeon, S., Belianinov, A., Kalinin, S. V., Baddorf, A. P., et al. (2016). Atomic intercalation to measure adhesion of graphene on graphite. Nat. Commun. 7:13263. doi: 10.1038/ncomms13263

PubMed Abstract | CrossRef Full Text | Google Scholar

Wick, P., Louw-Gaume, A. E., Kucki, M., Krug, H. F., Kostarelos, K., Fadeel, B., et al. (2014). Classification framework for graphene-based materials. Angew. Chem. 53, 7714–7718. doi: 10.1002/anie.201403335

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, L., Ma, T. B., Hu, Y. Z., and Wang, H. (2011). Vanishing stick–slip friction in few-layer graphenes: the thickness effect. Nanotechnology 22:285708. doi: 10.1088/0957-4484/22/28/285708

PubMed Abstract | CrossRef Full Text | Google Scholar

Yi, M., and Shen, Z. (2015). A review on mechanical exfoliation for the scalable production of graphene. J. Mater. Chem. 3, 11700–11715. doi: 10.1039/C5TA00252D

CrossRef Full Text | Google Scholar

Zhao, W., Wu, F., Wu, H., and Chen, G. (2010). Preparation of colloidal dispersions of graphene sheets in organic solvents by using ball milling. J. Nanomater. 2010:528235. doi: 10.1155/2010/528235

CrossRef Full Text | Google Scholar