Viscous peeling of a nanosheet

Abstract

Combining molecular dynamics (MD) and continuum simulations, we study the dynamics of propagation of a peeling front in a system composed of multilayered graphene nanosheets completely immersed in water. Peeling is induced by lifting one of the nanosheet edges with an assigned pulling velocity normal to the flat substrate. Using MD, we compute the pulling force as a function of the pulling velocity, and quantify the viscous resistance to the advancement of the peeling front. We compare the MD results to a 1D continuum model of a sheet loaded with modelled hydrodynamic loads. Our results show that the viscous dependence of the force on the velocity is negligible below a threshold velocity. Above this threshold, the hydrodynamics is mainly controlled by the viscous resistance associated to the flow near the crack opening, while lubrication forces are negligible owing to the large hydrodynamic slip at the liquid-solid boundary. Two dissipative mechanisms are identified: a drag resistance to the upward motion of the edge, and a resistance to the gap opening associated to the curvature of the flow streamlines near the entrance. Surprisingly, the shape of the sheet was found to be approximately independent of the pulling velocity even for the largest velocities considered.

---

1 Introduction

Two-dimensional (2D) nanomaterials, such as graphene, boron nitride, and molybdenum disulfide, are flexible structures of atomic thickness that can bend like a sheet of paper when exposed to sufficiently large forces. When a thin sheet is bound to a substrate by adhesion, external forces can induce peeling of the sheet, a complex phenomenon involving a competition between adhesive, bending, and dissipative forces.^1–8 The peeling of thin structures from rigid or soft substrates has received increasing attention from the soft matter community for its connections with soft wetting⁹ and its many applications in biology and engineering.^10–12 However, the specific properties of 2D materials,¹³ such as low bending rigidity and unusual surface properties, rise questions on the validity of current models, designed for macroscopic sheets, when applied to the description of the peeling of 2D materials.

In many applications, peeling of 2D materials occurs in the presence of a liquid.^3,14–17 For instance, in liquid-phase exfoliation processes for the large-scale production of graphene, colloidal microparticles of graphite are suspended in suitable liquid solvents (e.g. water, NMP), and shear is applied to the fluid–solid mixture until single- or few-layer graphene nanosheets detach from the ‘mother’ graphite microparticles.^14,15 Fluids can also affect peeling of 2D materials in applications that do not involve bulk liquid solution, such as the transfer of 2D materials between substrates, where water is often present owing to condensation from the surrounding air.^18,19 While the mechanics of peeling of 2D materials in vacuum or air has been studied extensively through theory and experiments,^2,20–23 peeling of 2D materials in liquids is a new subject.^24–26 Peeling off 2D materials from a substrate or a stack of other sheets in the presence of a liquid requires initially lifting ‘flaps’,²⁷ which can have nanometric length at the initial stages of the peeling process. The removal by peeling of 2D materials in liquids brings about a new set of scientific challenges, particularly considering that most available theories have been developed to explain the results of macroscopic applications such as adhesive tests or hydraulic fracturing.^1,3–5,8

The presence of a liquid primarily has two effects on the peeling of thin flexible sheets from a flat substrate: the liquid can alter directly the magnitude of the adhesion force between the bonded layers, by modifying e.g. the Hamaker constant;²⁴ and the liquid flow induces viscous forces on the peeled layer, which in turn affect the value of the peeling force.²⁸ For macroscopic sheets, the resistance to the motion of the peeling front originates mainly from lubrication forces.⁴ These forces emerge from the motion of the peeled layer in the direction normal to the substrate, which produces a parabolic flow in the small gap between the sheets. The lubrication forces associated to this flow are extremely sensitive to the boundary condition: even small deviations from the no-slip condition can reduce lubrication forces substantially, and this effect is most marked when the slip length is comparable to the gap size.^29–32 This raises the question as to what mechanisms determine the resistance to motion of the peeling front when the slip length is comparable or larger than system size, i.e. the maximum gap height or the typical crack length.

In this paper, we combine Molecular Dynamics (MD) and continuum modelling to study the peeling of short (few nanometers) graphene layers from a rigid substrate, for the case in which the maximum gap height and maximum crack length are smaller than the slip length characterising the hydrodynamic boundary condition at the solid-liquid interface.³³Fig. 1 illustrates the physical configuration simulated. In our problem, the rigid substrate simulated with MD is composed of a stack of graphene layers, a configuration motivated by applications to liquid-phase exfoliation of graphite microparticles. The peeled layer and the substrate are in contact with a liquid solvent (Fig. 1). As solvent, for this study we choose water because of the quality of published data on graphene–water interaction^29,34,35 and because water in MD displays Newtonian behaviour even at large shear rates^36,37 (in the ESI,† we also show results for NMP as solvent). Peeling is induced by displacing the edge of the flexible graphene layer with an assigned velocity v normal to the substrate. The force on the edge is measured. The objectives of the paper are to quantify the features of the force-velocity curve, and get insights into dissipative mechanisms controlling the dynamics of the peeling front. The slip length λ characterising the surface of pristine graphene in water and many other liquid solvents is rather large.³⁸ Slip lengths of (10 nm) have been reported for water-graphene interface both from experiments^39,40 and ab initio molecular dynamics methods.⁴¹ Therefore our MD simulation results, for which the crack length and maximum gap height are smaller than λ, cannot be explained by classical models based on lubrication theory. The paper therefore analyses alternative sources of viscous dissipation based on the comparison of MD results with different estimates for the viscous contribution to the pulling force.

Fig. 1 Illustration of the physical problem. (a) At time t = 0, a semi-infinite flexible nanosheet is completely bound to a rigid horizontal layer via van der Waals forces. The surrounding fluid does not enter the small gap between the sheets yet. The sheet is inextensible and cannot slide on the horizontal surface. (b) At time t ≥ 0, the edge of the flexible nanosheet is pulled upwards with a velocity v. The horizontal component of the applied force is zero (free horizontal sliding). The magnified view of the peeling front shows that fluid molecules penetrate only until the point where the separation between the nanosheet and the horizontal layer is too small to accommodate a fluid molecule.

Large slip lengths have been measured with many complex, structured fluids such as polymer melts, polymer solutions, colloidal suspensions, and colloidal gels.⁴² Slip is also significant in surfactant-covered solid surfaces when a high shear stress is applied to the fluid,⁴³ surfaces covered by polymer layers,⁴⁴ flow of rarefied gasses in microchannels,⁴⁵ and flow in nano-confined liquid systems.⁴⁶ The results of our investigation may thus have more general implications for soft matter research than our current focus on graphene may suggest. Our work is also relevant to understanding the effect of hydrodynamic forces in adhesion measurements. For example, in the measurement of the adhesion between surfaces in contact with high-viscosity ionic liquids, viscous forces are an important component of the force measured by a Surface Force Apparatus beyond a critical velocity that depends, among other parameters, on the slip length.⁴⁷ The prediction of this velocity is an important practical question.

The structure of the paper is as follows. We first analyse MD simulation results to quantify the dependence of pulling force on velocity. We then develop a one-dimensional non-linear solid mechanics model based on the equation for the elastica to obtain insights into the relation between pulling force, pulling velocity and edge height. We then compare the results between MD and the model in the quasi-static limit, and later consider the velocity dependent case. For the velocity dependent case, we complement the results with finite-element COMSOL simulations for a simplified geometry. In these simulations, the solid–fluid momentum coupling is fully resolved, so we can extract quantities – such as the full profile of the pressure below the peeling sheet – that the one-dimensional elastica model alone cannot provide.

2 Results and discussion

2.1 Molecular dynamics

The peeling simulations of a single graphene nanosheet from a multilayer graphene surface in water are carried out using MD (LAMMPS software⁴⁸). The TIP4P/2005 model⁴⁹ is used for water and the AIREBO force field⁵⁰ is used to model graphene, same as that used in literature to describe graphene-water systems.^26,38,46,51 Additional details of the simulation are given in the ESI.† The initial configuration consists of a stack of 4 graphene layers. The three bottom layers are periodic along the e⃑_x and e⃑_y directions. The top layer is periodic along e⃑_y, but shorter along the e⃑_x direction, thereby creating two edges. The graphene layers are in contact with water, and a piston wall⁵² is used to enclose the fluid in the e⃑_z direction and impose a pressure p₀ = 1 atm (Fig. 2a). The length of simulation box is equal to 7.2 nm in e⃑_x, 2.5 nm in e⃑_y and 8 nm in e⃑_z, and the distance between the piston and the top graphene layer is approximately 4 nm.

Fig. 2 Analysis of Molecular Dynamics (MD) simulations. (a) Side view of the MD system with water in red and white, graphene in pink, and the rigid wall in gray. (b) Close-up view of the MD system for 4 values of the peeled height (v = 1 m s⁻¹). (c) Peeling force as a function of edge height extracted from MD for different peeling velocities; v = 1 m s⁻¹ (green disks), v = 10 m s⁻¹ (red squares), and v = 50 m s⁻¹ (blue lozenges). The colored areas correspond to the standard deviation, and the black full line to the static case (v = 0). σ = 3.4 Å is the equilibrium sheet height. (d) Dynamic peeling force minus the static peeling force (F − F(v = 0)) as a function of v. Red disks are F_plt, i.e., F calculated at h₀ = 2 nm, and the blue squares are F_max.

The bottom three layers are maintained fixed, and the top layer is free to move. A velocity of magnitude v is applied along e⃑_z to one of the edges of the top layer. Velocities in the range of 1 to 100 m s⁻¹ are used. Lower values for v are difficult to reach due to computational limitation. For values of v ≥ 100 m s⁻¹, cavitation is observed (discussed further in Section 2.4). The other edge is maintained near its original position along e⃑_x by a harmonic potential. We make sure that restraining the motion of the top layer along e⃑_x has no impact on the final results by performing the same simulation in absence of the harmonic potential. We evaluate the sheet height along e⃑_y, h(x), as the vertical distance between the centers of the carbon atoms of the deformable sheet and the sheet immediately under it. The equilibrium sheet height, σ, is measured to be ≈3.4 Å when the sheets are completely adhered. Starting from the equilibrium height (σ), we measure the distance h₀ between the peeled edge and the fixed layer, and the force F (per unit width) along e⃑_z resisting peeling. We term F the peeling force and h₀ the edge height. In addition to dynamic simulations, we also perform static simulations by imposing a constant value for h₀ for 80 ps followed by an acquisition period of 120 ps.

The peeling force F was extracted from MD as a function of h₀. We find that as h₀ increases initially, F reaches a maximum, F_max (Fig. 2c). As h₀ increases further, F decreases to a plateau value F_plt. A similar trend was recently reported by Ouyang et al.⁵³ in MD simulations of peeling a graphene nanoribbon off a hexagonal boron nitride monolayer. Qualitatively, the trend of F vs. h₀ in Fig. 2c can be divided into three regions. In region I, F increases to a maximum F_max; in region II, F slowly decreases as h₀ increases; and in region III, F reaches a plateau value F_plt. We observe this qualitative trend for 3 values of the peeling velocity (v = 1,10 and 50 m s⁻¹.).

For v = 1 m s⁻¹, F is approximately equal to the force F⁰ obtained from the steady-state simulation (v = 0). For v = 10 m s⁻¹ and 50 m s⁻¹, the dependence on peeling velocity is evident (Fig. 2c). To ascertain that the increase of the force with the velocity is originating from the viscous dissipation within the fluid, we performed a simulation in absence of fluid with v = 100 m s⁻¹. This simulation showed no evident increase of the force with respect to the static case in vacuum (ESI,† Fig. S2). Both the characteristic maximum (F_max) and plateau (F_plt) values of the force increase with v. For v < 50 m s⁻¹, the difference F–F⁰ between the total peeling force and the peeling force at steady state increases approximately linearly with v (Fig. 2d). We observe that F_plt depends more strongly on v than F_max. For larger velocities, the rate of increase of F decreases. From Fig. 2d, we can estimate the order of magnitude of an effective friction coefficient ξ_eff ∼ (F − F_v=0)/v using F_max and F_plt as two extremes. We find 2.6 mPa s ≤ ξ_eff ≤ 7.6 mPa s. We notice that η_eff has the same order of magnitude of the viscosity of water (η = 0.855 mPa s for the TIP4P/2005 water model used in this study⁵⁴).

We compare snapshots of the system in the 3 regions of Fig. 2c. In region I the fluid does not intercalate between the adhered layers and no fluid molecules are present inside the crack (Fig. 2b). The slope of the sheet is small in this regime. In region II the water molecules start entering the crack. Here, h₀ is greater than ≈5 Å and the slope of the sheet increases with h₀, as can be seen in Fig. 2b. In region III, the sheet near the free edge is nearly vertical.

In the simulations, the fluid is in contact with the solid at all times for v < 50 m s⁻¹: the fluid molecules penetrate the crack until there is enough spacing between the carbon atoms to accommodate a fluid molecule, i.e., until a threshold sheet height of ≈5 Å. This motion is driven by pressure as the Péclet number (Pe) is larger than 1 for v > 10 m s⁻¹ and the fluid is practically incompressible (the Mach number for v = 100 m s⁻¹ is ). Here Pe ≈ v[small script l]/D_w, where D_w ≈ 2.3 × 10⁻⁹ m² s⁻¹ is the diffusion coefficient of water⁵⁵ and [small script l] ≈ 3 Å corresponds to one molecular diameter.

To interpret the results in Fig. 2, we build a continuum model inspired by the MD system. The model complements the MD results by analysing the balance of the adhesive, bending, and viscous forces for rate-dependent peeling.

2.2 A continuum model

We consider a continuum model for the peeling of an elastic sheet from a stationary rigid surface. The sheet has length L, thickness d, and bending stiffness B. The deformable sheet is bound to the stationary surface by an adhesion force related to an adhesion potential ϕ. The left edge of the sheet is pulled upwards with an assigned velocity v, while being allowed to move freely in the horizontal direction. This latter condition is enforced by setting the horizontal components of bending and tensile stresses at the edge to be equal (Fig. 3). We assume that the sheet is inextensible. Therefore, while moving upwards the edge moves in the e⃑_x direction so that the inextensibility condition is satisfied at all times. An incompressible fluid of viscosity μ fills the gap until a minimal threshold gap height is reached (see below). Through its motion, the fluid exerts a tangential hydrodynamic traction f on the deformable sheet, and determines a pressure difference ΔP between the bottom and top surface of the sheet.

Fig. 3 Schematic of the mathematical model. The deformable sheet is being pulled with a constant velocity v at the left edge. The blue shaded region represents the fluid of viscosity μ.

The shape of the sheet is described by a function h(s, t). When mapping to the MD simulations, h is the vertical distance between the centers of the atoms composing the deformable and stationary graphene sheets. Here s is the curvilinear coordinate line along the top sheet and t is the time. The corresponding Cartesian coordinates are x(s, t), h(s, t). Assuming that the sheet is inextensible, the angle of inclination of a point on the top sheet with the x-axis, θ(s, t), can be linked to h(s, t) and x(s, t) via

h_s = sin θ, x_s = cos θ, (1)

where ()_s represents the derivative with respect to s.

The equations of equilibrium governing the shape of the deformable sheet can be obtained from a force balance on an element ds (Fig. 3),⁵⁶ yielding

Bθ_sss + Tθ_s − ϕ_hcos θ = ΔP, (2)

and

Bθ_ssθ_s − T_s + ϕ_hsin θ = f. (3)

Here T is the axial tension in the sheet, ϕ_h is the vertical adhesion force per unit area, ΔP = P⁻ − P⁺ is the difference between the normal forces per unit area exerted by the fluid below and above the sheet, and f is the corresponding tangential force per unit area (directed towards the crack tip). Without loss of generality, the pressure P⁺ is assumed to be zero. When the sheet is pulled upwards, the fluid exerts a downward hydrodynamic tension on the sheet, hence ΔP<0. The potential ϕ(h) models the adhesion of the deformable graphene layer with the stationary layer. The potential ϕ(h) is taken as a standard 4-10 Lennard-Jones potential between two thin plates, consisting of an attractive term and a repulsion term:^57–59. Here A is the depth of the potential well and σ is the inter-layer equilibrium separation. Following previous work on graphene-graphene interactions,⁵⁸ we set σ = 3.4 Å, a value that is also consistent with our MD simulations.

For given ΔP and f, eqn (1)–(3) are solved to find θ, θ_s, θ_ss, T and h. Five boundary conditions are required: at the right edge s = L, θ_s = θ_ss = 0 (zero moment and zero vertical force);⁶⁰ at the left edge s = 0, θ_s = 0 (zero moment), h_t = v (kinematic condition) and [Bθ_sssin θ − Tsin θ]_s=0 = 0 (zero horizontal force). The upwards force acting on the left edge is

F = [Bθ_sscos θ + Tsin θ]_s=0 (4)

The adhesive and viscous contributions to F are calculated as

The finite van der Waals radius of carbon atoms makes the effective channel height available for the fluid molecules smaller than the channel height h defined from the center of the carbon atoms.^34,61,62 Accordingly, hydrodynamic boundary conditions are applied at a distance d/2 from the center of each sheet. The effective height of the nanochannel is thus h* = h − d, with d = 2σ_C ≃ 3.4 Å where σ_C denotes the van der Waals radius of a carbon atom. The finite range of the carbon–water interaction leads to a threshold for h below which no fluid molecule fits in the nanochannel as discussed in Section 2.1. We choose this threshold height, h_i = 5 Å for water–graphene,⁵¹ a value consistent with observation from MD profiles (Fig. 2b).

The threshold sheet height results into an interface position s_i beyond which there is no fluid (Fig. 3). Therefore we set ΔP = 0 for s > s_i. Due to the fluid flow, ΔP ≠ 0 left of this interface. Thus, there is a integrable discontinuity in the pressure across the interface at s_i.⁶³ The physics of this moving interface resembles that of the ‘dry cracking’ problem analysed by Lister et al. (2019), but is different from that of the ‘fluid lag’ concept,^64,65 which requires a vapour region ahead of the interface. To account for the pressure discontinuity, we divide the s-axis into two domains, to the left and to the right of the interface s_i. The solutions for the two domains are coupled to each other by enforcing continuity of θ, κ, κ_s, T and h across s_i.⁶⁰

Applying mass conservation and accounting for incompressibility we obtain

where ū is the height-averaged fluid velocity. Assuming that the liquid front position moves with a velocity equal to the average fluid velocity at that point, we arrive at the kinematic condition⁶³where the velocity is evaluated at positions just to the left of s_i.

The eqn (1)–(3), (7), (8) are solved together with the corresponding boundary conditions using the Boundary Value Problem (BVP) solver of MATLAB.⁶⁶ As initial condition, we use the stationary solution corresponding to a small assigned edge displacement h(s = 0) = h_i. The initial shape of the sheet is obtained by running the iterative solver with v = 0 starting from the assigned shape

until convergence. The regions to the left and to right of s_i are discretized with two uniform meshes. We assume that the interface always moves in the positive x-direction as the sheet is peeled. Therefore, at each step of the iteration, the interface is moved from the previous position by a small amount equal to the grid spacing by using an implicit time-marching method. We repeat the above steps until F reaches a constant value. For the range of parameters we consider, this usually happens when the edge slope reaches π/2.

2.3 Comparison MD-continuum in the quasi-steady case

In this section, we consider quasi-steady peeling where viscous forces are absent (ΔP = f = 0) and therefore the forces on the sheet are rate independent. In the following we focus particularly on the prediction of the quasi-steady plateau force F⁰_plt and the maximum force F⁰_max obtained by solving eqn (1)–(3). The computed values of the force are compared with MD data to calibrate values of A and B.

In the limit h₀ ≫ σ, the left edge has quasivertical orientation (θ (s = 0) ≈ π/2), the curvature of the sheet near the crack tip is approximately independent of h₀⁶ and the bending energy is constant with respect to h₀. The energy balance thus involves only the external work done by F⁰ and the adhesion work: ^6,14,67

F⁰_plt ≈ A. (10)

By calculating the work of adhesion required to separate two parallel sheets in vacuum using MD²⁵ we obtain A ≃ 0.29 N m⁻¹, a value close to the force plateau in vacuum ≃0.28 N m⁻¹ (Fig. 4a).

Fig. 4 Quasi-steady peeling. (a) Dependence of F⁰ on h₀, comparing the 1D continuum model (lines) and the MD data (solid symbols). The plot shows continuum simulations for 3 combinations of adhesion energy (A) and bending rigidity (B) and, MD simulations in vacuum and water as solvent. (b) F⁰vs. h₀ calculated from eqn (12) and (15) compared with numerical solution for A = 0.28 N m⁻¹ and B = 1 eV. (c) Dependence of h_x(0) on 1/a₀ for h₀ = 1.26σ calculated numerically and from eqn (14). (d) F⁰_maxvs. and the dotted line is the plot for y = 1.18x. The top and bottom insets show the dependence of F⁰_max (N m⁻¹) on A (N m⁻¹) and B (eV) in log–log scale.

To obtain an approximate analytical solution for F⁰_max, we exploit the fact that the maximum value of the force is attained when slope of the sheet is small. For small slopes, the shape of the deformable sheet is governed by:^4,5,68

Bh_xxxx − ϕ_h = ΔP. (11)

The edge force F⁰ ≃ Bh_xxx can be obtained by multiplying eqn (11) by h_x, integrating with respect to x and then using the boundary condition κ(0) = 0. The same result can be obtained by using the principle of virtual work.⁶⁹ For ΔP = 0, we getInserting values of h_x(0) in the above equation from numerical resolution of the continuum model, we get an excellent agreement with the complete numerical solution of F⁰ (Fig. 4b). We now express h_x(0) as a function of h₀, so that we can maximise F⁰ with respect to h₀ to obtain F⁰_max. Linearising ϕ_h in eqn (11) about the equilibrium separation giveswhich has the solution h − σ = (h₀ − σ) e^−x/a₀cos(x/a₀), where is the cohesion length. ThusPlugging this relation in eqn (12) yieldsUsing the expression for ϕ(h) (in Section 3), the above equation gives the same trend as the complete numerical solution (Fig. 4b) and has a fair agreement for the value of h₀ for which F⁰ is maximum (h₀ ≃ 1.26σ). For σ/a₀ < 2, a condition which often holds in the context of thin film peeling experiments, the numerical solution confirms the 1/a dependency of h_x(0) in eqn (14) at h₀ = 1.26σ (see Fig. 4c). Maximizing F⁰ with respect to h₀ giveswhere the constant of proportionality is . We confirm the power-law relation of A and B with F⁰_max in eqn (16) by numerically computing F⁰_max for A in range [0.1,10] N m⁻¹ (keeping fixed B = 1 eV and σ = 0.34 nm) and for B in range [0.1,10] eV (keeping fixed A = 0.28 N m⁻¹ and σ = 0.34 nm) (Fig. 4d insets). As a further confirmation, F⁰_max is computed by varying σ in range [0.1,0.5] nm and A, B in the aforementioned range. Expressing F⁰_max as a function of , the data collapses onto a straight line (Fig. 4d). We estimate the proportionality constant to be 1.18, close to the value of 0.88 obtained by exact minimisation of the small-slope expression. In summary, the maximum force is given approximately byRef. 27,70 showed that for small opening angle, , where [small gamma, Greek, dot above] is the critical shear rate for exfoliation. Using F⁰_max∼ μ[small gamma, Greek, dot above] a and a∼ a_o, we arrive at the same scaling as in eqn17. This suggests that for a_o > 0.34σ, which turns out to be the case for the parameters B, A and σ characterising our simulations, the force F⁰_max is larger than F⁰_plateau by a factor . In peeling of macroscopic elastic films B is much higher than in graphene.^4,71–74 Our predictions for the maximum and plateau forces are based on a continuum framework, and therefore hold for macroscopic systems provided that a₀ ≫ σ.

The values of F⁰_max and F⁰_plateau depend on A and B. We can therefore use the quasi-static estimates of peeling force to calibrate these two parameters. Using eqn (10) for peeling in vacuum we obtain A ≃ 0.28 N m⁻¹, a value that is very close to the graphene-graphene adhesion energy reported in previous MD simulations and experiments.^75,76 For peeling in water, we obtain A ≃ 0.25 N m⁻¹, comparable to the value in vacuum.^24,25 Previous experiments and MD simulations estimated the adhesion energy of graphene in water to be of the same order magnitude as found here.^3,25,77,78 Knowing the adhesion energy, eqn (17) can be used to calculate the bending stiffness. The value we obtain is B = 2.4 eV. However for h₀ > 0.8 nm fitting the continuum peeling force with MD using BFGS method⁷⁹ yields B = 1.0 eV. This suggests that the bending stiffness varies between 1.0–2.4 eV depending on h₀. Previous works have also reported a variation of bending rigidity in a similar range, depending on the extent of deformation of the sheet.^80,81

We compare sheet profiles in the continuum model with those obtained from MD. The quasi-static evolution of sheet shape in MD shows good agreement with that obtained by solving eqn (2) and (3) with the fitted MD values of B and A for vacuum (Fig. 5a). The sheet profiles in Fig. 5b and c confirm our observation on the variance of B on h₀: B = 2.4 eV shows a slightly better fit for h₀ = 0.355 nm while B = 1.0 eV is a considerably better fit for h₀ = 1.75 nm. In the cohesive region, the oscillatory pattern visible in Fig. 5b is a consequence of the competition between the bending and adhesive forces.⁸² In the following we use B = 1.0 eV and A = 0.25 N m⁻¹ when discussing the values of F for graphene–water system.

Fig. 5 Comparison of quasi-steady sheet profiles from the continuum model (lines) and MD simulations (disks) in vacuum. (a) Sheet profiles for 4 different values of h₀, (0.355, 0.52, 1.04, 1.75) nm. The parameters A = 0.28 N m⁻¹ and B = 1 eV for the continuum profiles. (b and c) The MD profiles for h₀ = 0.355 nm and h₀ = 1.75 nm respectively, compared with continuum profiles for two values of B. Here, A = 0.28 N m⁻¹.

2.4 Analysis of the velocity-dependent case

In this section, we analyze the effect of rate-dependent peeling on the sheet's shape and peeling force. To get insights on the dependence of the shape of the sheet on the peeling velocity, we first write eqn (2) and (3) in dimensionless form. With the horizontal length of the channel scaling with and estimating the characteristic curvature as ,^56,83,84 the characteristic slope near the crack tip is . We rewrite eqn (2) and (3) using the following dimensionless variables (represented with hat symbols):The dimensionless governing equations arewhere and ĥ = h/(ε a₀). The dimensionless sheet thickness is d̂ = d/(εa₀) and the dimensionless pulling velocity is v̂ = μv/(2A) (With the normalisation eqn (18), P̂ and f̂ are the ratios of the pressure and friction terms to the first terms in (19) and (20), respectively. For ε ≪ 1, the small slope limit of the balance equations is recovered.). Characteristic values from our MD simulations are μ = 10⁻³ Pa s and A ∼ 0.1 N m⁻¹, hence v̂ turns out to be in the range [0.005–0.5] for v ∈ [1–100] m s⁻¹. In eqn (19) and (20), ΔP̂ is an increasing function of v̂, so there is a lower range of value of v̂ for which ΔP̂ ≪ 1 and f̂ ≪ 1, i.e. the quasi-static case.

Examining the effect of v on the sheet shape extracted from MD (Fig. 6a), we noticed that the profile of the sheet is practically independent of v for the range of velocities we considered. To interpret this observation, we solve the coupled eqn (19) and (20) with a simple prescription for ΔP̂ in which ΔP̂ = −n for ĥ greater than the interface height ĥ_i and zero otherwise:

Here we fix ε = 0.5, f̂ = 0 and n = 0 or 1. If the slope is not negligible, an increase in peeling force with ΔP̂ must be accompanied by an increase in tension at the left boundary (eqn (4)). From eqn (20), in the case when [small straight theta, Greek, circumflex] and its first two derivatives are independent of ΔP̂, T̂_ŝ remains nearly constant, and therefore the tension curve must uniformly shift upwards with the increase in tension at the boundary. Considering that for n = 1, there is a significant increase in F̂, we find that T̂ indeed shifts upward with n for ĥ_i = 3 (Fig. 6b). However, the increase for ĥ_i = 2.3 is notably non-uniform implying a significant change in shape. A smaller value of ĥ_i means that the fluid wets further inside the sheet at the peeling front. Therefore, the unvarying shape of the sheet with v can be explained by the absence of fluid in the curved part of the sheet near the front (see Fig. 2a).

Fig. 6 Rate dependence of sheet profile and tension. (a) Sheet profiles from MD for h₀ = 1.75 nm at v = 0 (squares) and v = 50 m s⁻¹ (disks). Dotted line shows solution of continuum model (A = 0.25 N m⁻¹, B = 1.0 eV) with ΔP = f = 0. (b) Non-dimensional tension along the sheet for the pressure profile in terms of n defined in eqn (21) with two values of ĥ_i and ĥ₀ = 4.

The fact that the shape of the sheet is approximately independent of the peeling velocity even in presence of water means that we can estimate the bending and adhesive forces from the quasi-static results, and use these estimates when the velocity is not small. As a corollary, we can decompose unambiguously the peeling force F into a non-dissipative contribution, given by eqn (10) for large angles and eqn (12) for small angles, and a velocity-dependent viscous contribution.

As customary in elasto-hydrodynamic problems, we investigate the relation between pressure-drop and axial velocity in the lubrication limit of small slopes.^5,85–87 From the continuum model, we analyse the velocity dependent case for small deflections (i.e., h₀ < 1 nm) as the 1-D approximations for pressure are valid only for small slopes. In this limit, for stationary channels presenting Navier-slip boundary conditions at the walls, the depth-averaged fluid velocity is related to the pressure gradient via⁸⁸

which recovers the no-slip expression for λ ≪ h*/6. Here the height of the channel is taken to be h* = h − d, to account for the thickness of the sheet (see Section 2.2). Eqn (22) is in the form used, e.g., to model the flow of thin liquid films in the weak-slip regime.^89–91 Momentum balance in the flow direction giveswhere f = μ(du⃑/dy)·n⃑ is the tangential friction and n⃑ is the unit vector normal to the boundary. We solve eqn (2) and (3) with the prescription of pressure and friction described above^4,5 and p|_s=0 = 0 (Fig. 7). From the shape of the sheet, we calculate the edge force using eqn (4).

Fig. 7 Results of lubrication model described by eqn (22) and (23). (a) Comparison of lubrication pressure profile along the sheet (h₀ = 0.8 nm) from the continuum model for different combinations of λ and d. (b) F^vis (solid lines) and F^adh (dashed lines) as a function of λ using eqn (22). The peeling velocities, v = 50 m s⁻¹ (red) and v = 10 m s⁻¹ (green), h₀ = 0.8 nm and B, A, σ, d are for graphene-water system. The inset is an enlarged version of the same plot showing the difference in F^adh values.

If the slip length was zero, the crack propagation would result into a diverging viscous stress near the crack tip,^85,92 as in the moving contact line problem.^93,94 In theory of moving contact line,^93,95 introduction of slip regularises the stress divergence. The slip lengths are comparatively large in our MD simulations (λ ≃ 60 nm ≫ a₀ ≃ 0.3 nm)³⁸ and significantly reduce the magnitude of the viscous stress near the crack tip (Fig. 7a). The choice of d has no noticeable effect on lubrication pressure as λ ≫ h− d.

For large slip lengths (6λ/h* ≫ 1), the viscous component of F due to lubrication is much smaller than the adhesive component even at relatively large peeling velocities (Fig. 7b). As an instance, for λ = 10 nm and v = 50 m s⁻¹ the viscous component F^vis is about one order of magnitude smaller than F^adh. If the shape of the sheet is invariant with respect to v, the lubrication pressure for λ ≫ h* should scale, according to (22), approximately as

Here, we use the continuity equation for small slopes, which yields . Therefore, F_vis ∝ 1/λ for sufficiently large λ, which our continuum simulation confirms (Fig. 7b). The criterion to neglect viscous stresses due to lubrication (P̂ ≪ 1) yields the following “large-slip” condition:

Using eqn (25), we find that the viscous component of F arising from lubrication can be neglected for v ≪ 10⁵ m s⁻¹ in our MD simulation. Thus, lubrication forces cannot explain the dependence of peeling force on velocity seen in the MD results of Fig. 2d.

In the MD simulations of confined liquid in nano-channels,^96–98 a large hydrodynamic slippage is known to generate a more uniform flow profile (the usual parabolic flow expected for pressured driven flow in no-slip channels reduces to a uniform flow as the slip length increases^29,99). Similarly, in the problem considered in the current paper the effect of large slip is to make the velocity field near the crack tip approximately uniform (Fig. 8). The streamlines (in the frame of reference of the stationary layer) near the crack-tip are parallel to the lower boundary. The velocity field does satisfy both the no-penetration and the tangential slip boundary conditions at the liquid–solid boundaries. A small e⃑_z component of the fluid velocity can only be seen close-to the edge, above the sheet and in front of the entrance of the flap.

Fig. 8 Fluid velocity vector field computed from molecular dynamics simulation. The blue lines denote velocity streamlines. The data is averaged over 8 simulations at v = 50 m s⁻¹.

Calculating pressure, viscous shear and normal stresses from MD is challenging due to thermal noise and the difficulty of applying volume-average to smallgap regions. Therefore, we have performed continuum simulations using a finite element software (COMSOL) of a peeled sheet of finite thickness for 4 values of h₀ moving at different peeling velocities (Fig. 9a). In these simulations, we solve the incompressible Stokes equations

∇p = μ∇²u; ∇·u = 0 (26)

where p is the pressure field and u the velocity field, with free slip boundary condition at all surfaces, corresponding to λ = ∞. The boundaries, ABCD in Fig. 9a, are prescribed based on the solution of eqn (2) and (3) for ΔP = f = 0 (see Fig. 6a). The boundaries are placed at a distance d/2 from the centre-line positions of the moving and stationary sheets, given by y = h(x) and y = 0 respectively. To avoid singularity, point B and C are replaced by rounded corner of radius r = 0.5 Å. To avoid finite size effect due to the reservoir, the height (FG) and width (GH) of the computational domain is chosen to be much larger than any other dimension. The interface (segment DE in Fig. 9a) of height h_i− d, located at position s = s_i, moves with velocity u_i in the e⃑_x direction, where u_i = ds_i/dt is the average fluid velocity at s_i. The reservoir walls are modelled as zero pressure gradient outlets. We perform simulations for h₀ = 1, 2, and 5 nm.

Fig. 9 Analysis of COMSOL simulations of the dynamically peeled sheet. (a) Schematic of the simulation for h₀ = 1 nm. The region with fluid is shaded with grey. The reservoir boundaries are shown not to scale. (b) Representative pressure and velocity field. (c) Absolute pressure just below the sheet (along CD) extracted from the simulation for different peeling heights for v = 50 m s⁻¹. (d) Simulation schematic illustrating the configuration of sheet for h₀ = 5 nm and the corresponding velocity streamlines. The surface near the edge of the sheet (corresponding to y > 4 nm) is highlighted in green.

Fig. 9b shows the fluid velocity vectors and the corresponding pressure distribution for v = 50 m s⁻¹ and h₀ = 1 nm. The pressure in between the gap is predominantly negative, with a large negative pressure in the region surrounding the edge. The absolute pressure along the arc CD has a sharp minimum near the edge and then increases smoothly as the crack tip is approached (Fig. 9c). A convective motion of fluid around the edge of the sheet is visible in Fig. 9b and 8. The integral effect of the viscous and pressure stresses associated to this convective motion give rise to a drag resistance, which we term “edge drag” force. To characterise the magnitude of this force, we have developed two independent estimates. One is based on the COMSOL simulations discussed previously. We consider the case of h₀ = 5 nm shown in Fig. 9d for different values of v as the viscous effect of the edge can be isolated in this configuration. To estimate the viscous stress at the edge, we integrate the y-component of stress along the surface of the sheet above y = 4 nm (highlighted in green in Fig. 9d). The estimated viscous stress increases approximately linearly with v and is comparable in magnitude with A for v > 10 m s⁻¹ (Fig. 10a). As a second estimate, using MD simulations we have computed the vertical force on a flat rigid vertical nanosheet of length ∼ 1 nm moving vertically with an assigned velocity v.‡ In both cases we get values which are very close, in magnitude and trend, to the value of F_max − F⁰_max. Both estimates give values of the edge force comparable in order of magnitude to the quasi-steady adhesion forces for v > 10 m s⁻¹. Therefore, the sharp minimum in (Fig. 9c) can be attributed to the motion of fluid displaced as the edge moves upward.

Fig. 10 (a) Viscous resisting force for a single nanosheet moving upwards with velocity v (red disks); F_max − F⁰_max from MD (blue squares); integral of the y-component of the stress near the edge of the sheet from COMSOL integrated along the surface highlighted in green in Fig. 9d (line). Here, A = 0.25 N m⁻¹ is the adhesion energy measured in water. (b) Estimate of F^ent from eqn (27) for h₀ = 0.6, 0.8 and 1 nm. Red disks are F − F⁰, i.e., the total increase in peeling force in MD at h₀ = 1 nm.

The “edge drag”' force does not capture all the viscous forces as it does not explain the increase in F − F⁰ occurring when fluid starts entering the crack (region II of Fig. 2c). The fact that the pressure does not suddenly increase to zero and saturates to a significant value for h₀ = 1 nm (Fig. 9c) suggests that another source of pressure drop in our system can be due to the motion of fluid entering the gap between the sheets. To estimate this contribution, we refer to previous results on pressure-driven flow in two-dimensional channels which share similarities to the entrance flow below the sheet and near the edge in our configuration. Hasimoto (1958)¹⁰⁰ studied analytically the Stokes-flow hydrodynamic resistance of an infinitely thin plate presenting a slot of height 2h. This configuration is relevant to our case when h₀ is small and the channel walls are thus nearly parallel. The pressure drop along the channel was characterised in terms of the hydrodynamic resistance R_H, i.e. the ratio of the magnitude of the pressure drop to the volumetric flow rate Q ∼ ūh. Hasimoto's solution gives

Hasimoto's formula is not an accurate representation of hydrodynamic resistance for smoothly-converging channels, so the accuracy of eqn (27) decreases for increasing h₀. Numerical studies of pressure-driven flow in smoothly converging axis-symmetric entrances with fixed walls show that R_H depends on the ratio of λ to the minimal channel height a, the ratio of the channel length and a, as well as the specific shape of the channel.^101–104 For completeness, a COMSOL analysis of the entrance resistance R_H for pressure-driven flow in smoothly-converging stationary 2D channels (as opposed to axi-symmetric) in the case of infinite slip length, for different values of the radius of curvature of the entrance, is presented in the ESI.† The analysis show that Hasimoto's solution gives values of R_H only 20% smaller than those provided by COMSOL when the radius of curvature of the entrance is comparable to, or smaller than, the minimum channel height (see Fig. S4 in ESI†). Therefore, we use eqn (27) as an approximation for small values of h₀. We notice that a viscous resistance due to converging streamlines is, physically, akin to the resistance due to extensional viscous stresses described in the lubrication formulation of ref. 90, 105 and 106, except that in our case the gradients are large because the fluid layer is not slender.

We compute ΔP, appearing in eqn (2), using eqn (27) with the available channel height in place of h, and evaluating Q from the average velocity defined in eqn (7). We calculate the component of the force associated with entrance flow, F^ent, from the total force by removing the adhesive component (eqn (6)). We find that F^ent is comparable to, but smaller than, the value F − F⁰ extracted from MD at h₀ = 1 nm (Fig. 10b). This suggest that entrance flow is an important contribution to the viscous resistance to the peeling front motion.

As F_max lies in region I of Fig. 2c, we infer that the viscous contribution to F_max must arise from the motion of the fluid mostly close to the edge, because there is practically no fluid in between the sheets in region I. The peeling force in region II & III will instead receive contributions from both the edge and the entrance resistance. Indeed, the difference between F − F⁰ and the theoretical values for F_ent in Fig. 10b is comparable to the viscous force due to the edge drag in Fig. 10a. So, the sum of these two contributions is close the total viscous resistance (not exactly equal, expectedly, due to the use of simplified models). Using the COMSOL estimate for the edge drag force and eqn (27) for the entrance pressure drop in our continuum model, we obtain an increase in F with v for v > 1 m s⁻¹. This increase in F is qualitatively similar to that suggested by the MD data in Fig. 2c. To build an accurate model, one needs to separate the edge and entrance contributions to the viscous resistance in MD. However, a systematic way to separate these two contributions was unfortunately not found since both contributions are, to a first approximation, linear in the velocity and their magnitudes are also comparable.

The MD results in Fig. 10a and b show approximately a linear trend with v, except at v = 100 m s⁻¹. For this particular value the MD data for the force is lower than expected by extrapolation of the force-velocity curve at smaller velocity. Fig. 11 illustrates that this velocity value corresponds to the formation of a vacuum pocket (cavitation) near the crack tip. (At this velocity value, capillary forces may play an important role.) Cavitation at large peeling velocity have been reported^10,107 in the case of peeling of a scotch tape. In molecular dynamics, the threshold negative pressure for cavitation in pure water is reported to be ∼0.2 GPa,^108,109 which is of the same order as the pressure near the crack measured from our COMSOL calculation.

Fig. 11 Snapshot of MD simulation near the crack tip at different peeling velocities for h₀ = 1.6 nm.

2.5 Threshold velocity

In the initial stage of peeling, using eqn (13) the slope at the entrance is . Balancing the terms in the fluid continuity equation gives . Thus,To obtain this expression we have here used the scaling suggested by Hasimoto's solution (27). Using a practical numerical threshold of 0.01 in place of ≪ 1, the condition ΔP̂ ≪ 1 translates to a practical threshold ofto be compared with the lubrication estimate (25). The entrance-flow estimate (29) gives v_th ∼ 1 m s⁻¹ in our problem of peeling of mono-layer graphene in water.

The parametric dependence of the threshold velocity using the entrance pressure is different from the one obtained using the large-slip lubrication pressure. This is because the typical length-scale characterising the velocity gradients in large-slip lubrication is , while it is in entrance flow. Fig. 12, obtained with the 1-D continuum model, confirms the ∼B^−1/2 and ∼B^−1/4 dependencies of the threshold velocity in the large-slip lubrication case and in the entrance flow case, respectively. In this test we define v_th as the velocity for which F increases by 1% with respect to the quasi-steady value. In the literature on silicon wafer bonding in air, the velocity of crack has been reported to have the same scaling in B and A as in eqn (29).⁷¹ An inverse 3/2 dependence with sheet thickness has also been previously reported.¹¹⁰

Fig. 12 Threshold peeling velocity. (a) v_th for entrance model (cross symbol) and lubrication model with λ = 60 nm (plus symbol) as a function of B, keeping A = 0.25 N m⁻¹, σ = d = 0.34 nm, h₀ = 1 nm and h_i = 0.5 nm. (b) F vs. v for different values of B for entrance model, keeping other parameters same as (a). The value v_th (cross symbol) is defined as the point where F = 1.01F⁰.

Eqn (29) predicts an increase in threshold velocity with adhesion energy. The adhesion force increases with A and the peeling velocity at which viscous forces start becoming comparable to adhesion forces also increases. The inverse dependence on σ can be explained similarly: for the same value of , adhesion forces increase with decreasing σ. Eqn (29) predicts a weak inverse dependence on bending rigidity. A more rigid sheet leads to smaller deformations, larger crack lengths and larger fluid velocity at the entrance. Hence the total viscous force on the top sheet is increased. This increment in viscous force must be larger than the increment in the quasi-steady peeling force (with B) at small slopes. Our continuum model indeed confirms that F, unlike v_th, increases with B (Fig. 12b).

Comparing the scaling for lubrication pressure and entrance pressure (eqn (24) and eqn (28) respectively), we get the condition λ/a₀ ≫ 1 to neglect the effect of the lubrication pressure in comparison with the entrance pressure for small peeling angles.

2.6 Comparison with the literature

For the static peeling, we define a horizontal length scale a₀ from the balance between adhesive and bending stresses near the crack tip. Previous studies^28,111 suggest that the stresses at a length near the peeling front control the shape of the peeled flap. In our case, a₀ ∼ 0.3 nm, comparable to the size of one fluid molecule, and comparable to the length of the curved part of sheet near the peeling front. As fluid is absent in the curved part for s > s_i, the shape of the peeling front in the neighbourhood of this region is approximately quasi-static, with characteristic outer curvature .¹¹² Away from the front, the shape is practically linear owing to the zero moment condition at the edge. Correspondingly, the shape of the sheet is approximately independent of viscous effects. Our case is contrary to that of an elastic polymer sheet with a pre-wetting layer of liquid underneath, where the viscous effects near the peeling front control the sheet's shape.^1,4,63

In previous studies on viscous peeling, the fluid is either present only inside the crack^1,60,63 or confined within a blister;^4,28,68 in these studies lubrication dissipation due to a Poiseuille flow in the gap was the primary source of viscous dissipation. In contrast, our configuration is completely immersed in the fluid. Therefore, viscous dissipation also comes from the motion of the edge relative to the surrounding fluid²⁷ and suction of fluid from the outer “reservoir” as the opening angle increases.¹⁰¹ When the slip is large compared to the system's dimensions, these forms of dissipation dominate over the usual lubrication dissipation.

3 Summary and outlook

We have carried out MD simulations of peeling of a deformable graphene layer in the presence of liquid water. The MD data is rationalised by a 1-D continuum model based on equilibrium equations for large deflections of a non-linear elastica. These solid mechanics equations have been used to model, e.g., inextensible elastic rods^113,114 and fibres,¹¹⁵¹¹⁶ also for large deflections. Finite-element COMSOL simulations are also carried out to investigate certain flow features and to better characterise the hydrodynamic load terms present in MD. The effect of identified viscous stresses are estimated, for small deflections, using the 1-D model.

In the quasi-static regime, we have been able to capture with the 1-D model the dependence of the peeling force on the edge height h₀, provided that two values were used for the bending rigidity, one for small deformations, and another for large deformations. This feature is likely due to the dependence of the mechanical response of the crystal structure of graphene to the tension applied to it.¹¹⁷ The quasi-static force profile as a function of h₀ displays a maximum for intermediate values of h₀, and then reaches a plateau. The plateau force in MD was found to scale linearly with the adhesion energy A, independent of the bending rigidity B of the sheet. The maximum quasi-steady peeling force instead was found to scale proportionally to B^1/4 and A^3/4. In the velocity dependent cases, the analysis, corroborated by finite-element COMSOL simulations, suggests that entrance flow effects (associated to the curvature of the streamlines at the entrance of the nanochannel), and the vertical drag force exerted by the fluid on the edge of the sheet as it moves upwards, are dominant contributions to the viscous resistance to the motion of the sheet.

For our range of parameters, the critical velocity below which the deformation is quasi-static is about 1 m s⁻¹. An estimate for the critical velocity beyond which entrance flow effects make the peeling force velocity dependent was developed (eqn (29)). This estimate is valid provided that the slip length is much larger than the length scale a₀. The dependence of the shape of the sheet on v was found to be practically negligible even in the regime where the peeling force displayed a significant dependence on v.

For very large peeling velocities, the MD simulations reveal the formation of a cavity near the crack tip. The formation of this cavity was shown to correlate with a decline in the trend of peeling force vs. velocity curve. Snapshots from the MD simulations (Fig. 11) show quite evidently that this decline is simply due to the reduction in contact area between the liquid and the solid, so it would be erroneous to associate the reduction in viscous resistance to a variation in slip or viscosity at a critical velocity (cavitation for instance has been used to explain the appearance of sudden slip in atomically smooth mica surface coated with surfactants or self-assembled monolayers⁴³).

In the initial stages of peeling of a graphene layer, the linear dimension of the peeled flap can be comparable to the hydrodynamic slip length characterising the liquid-solid interaction (the slip length is a few tens of nanomaters for many graphene solvent combinations³⁸). For these small system sizes, the peeling force and the dynamics of the sheets will be controlled by mechanisms similar to the ones discussed in this paper, so our results can be used, for example, to build analytical models to predict thresholds for graphene exfoliation in shear mixing.^{14,27,70,75,118}

A possible future development of the current work is the inclusion of thermal fluctuations in the elastic, continuum model. Fluctuations for instance could play an important role in peeling initiation (methods to include the effect of thermal fluctuations in peeling models are available^119–121). Another relevant avenue of research would be to consider all-atom or coarse-grained MD simulations simulations accounting for longer nanosheets that we were able to simulate here. This could clarify whether a transition from edge- and entrance-dominated viscous dissipation to lubrication-dominated viscous dissipation occurs at a critical value of the crack length.

The current work may have implications for other soft matter systems. Peeling processes have been shown to play an important role in a variety of soft matter problems, from the adhesion of cells¹²² and biological membranes,¹²³ to the mechanics of pressure-sensitive adhesives.¹²⁴ Much of the work done on peeling and thin-film lubrication focuses on phenomena occurring near the crack tip,⁷¹ or in the bulk fluid away from the entrance.¹⁰⁵ The current work instead highlight the importance of entrance effects, whether associated to the entrance flow or the downward drag on the edge. Entrance effects will be particularly important for soft matter systems such as polymer melts⁴² that can display slip length sometimes exceeding hundreds of microns.¹²⁵ Slip lengths tend to increase as the surface roughness decreases, thus significant slip effects are expected for molecularly smooth surfaces, graphene being a primary example but certainly not the only one.¹²⁶ For both no-slip and slip sheets peeled from a substrate and initially completely bound to it, entrance effects will be more pronounced during the initial stages of peeling, when the interfacial crack is short and the hydrodynamics of the entrance dominates the flow everywhere below the peeled sheet.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was funded by the European Research Council, under the European Union's Horizon 2020 Research and Innovation programme (project FLEXNANOFLOW, grant agreement no. 715475). This research utilised Queen Mary's Apocrita HPC facility, supported by QMUL Research-IT.¹²⁷

References

1. A. D. McEwan, Rheol. Acta, 1966, 5, 205–211.
2. E. Memet, F. Hilitski, Z. Dogic and L. Mahadevan, Soft Matter, 2021, 17, 2704–2710.
3. M. Z. Miskin, C. Sun, I. Cohen, W. R. Dichtel and P. L. McEuen, Nano Lett., 2018, 18, 449–454.
4. J. R. Lister, G. G. Peng and J. A. Neufeld, Phys. Rev. Lett., 2013, 111, 1–5.
5. A. E. Hosoi and L. Mahadevan, Phys. Rev. Lett., 2004, 93, 8–11.
6. B. Roman, Int. J. Fracture, 2013, 182, 209–237.
7. A. Juel, D. Pihler-Puzović and M. Heil, Annu. Rev. Fluid Mech., 2018, 50, 691–714.
8. C. Creton and M. Ciccotti, Rep. Prog. Phys., 2016, 79, 046601.
9. B. Andreotti and J. H. Snoeijer, Annu. Rev. Fluid Mech., 2020, 52, 285–308.
10. R. Villey, C. Creton, P.-P. Cortet, M.-J. Dalbe, T. Jet, B. Saintyves, S. Santucci, L. Vanel, D. J. Yarusso and M. Ciccotti, Soft Matter, 2015, 11, 3480–3491.
11. A. Ibarra, B. Roman and F. Melo, Soft Matter, 2016, 12, 5979–5985.
12. D. A. Dillard, B. Mukherjee, P. Karnal, R. C. Batra and J. Frechette, Soft Matter, 2018, 14, 3669–3683.
13. Q. Cao, X. Geng, H. Wang, P. Wang, A. Liu, Y. Lan and Q. Peng, Crystals, 2018, 8, 357.
14. J. Stafford, N. Uzo, U. Farooq, S. Favero, S. Wang, H.-H. Chen, A. L’Hermitte, C. Petit and O. K. Matar, 2D Mater., 2021, 8, 025029.
15. X. Chen, R. A. Boulos, J. F. Dobson and C. L. Raston, Nanoscale, 2013, 5, 498–502.
16. G. F. Schneider, V. E. Calado, H. Zandbergen, L. M. Vandersypen and C. Dekker, Nano Lett., 2010, 10, 1912–1916.
17. J. Seo, C. Kim, B. S. Ma, T.-I. Lee, J. H. Bong, J. G. Oh, B. J. Cho and T.-S. Kim, Adv. Funct. Mater., 2018, 28, 1707102.
18. B. Li, A. V. Klekachev, M. Cantoro, C. Huyghebaert, A. Stesmans, I. Asselberghs, S. De Gendt and S. De Feyter, Nanoscale, 2013, 5, 9640–9644.
19. J. Wu, L. Xie, Y. Li, H. Wang, Y. Ouyang, J. Guo and H. Dai, J. Am. Chem. Soc., 2011, 133, 19668–19671.
20. J. Annett and G. L. Cross, Nature, 2016, 535, 271–275.
21. M. Ishikawa, M. Ichikawa, H. Okamoto, N. Itamura, N. Sasaki and K. Miura, Appl. Phys. Express, 2012, 5, 065102.
22. R. C. Sinclair, J. L. Suter and P. V. Coveney, Adv. Mater., 2018, 30, 1705791.
23. C. Ke, M. Zheng, G. Zhou, W. Cui, N. Pugno and R. N. Miles, Small, 2010, 6, 438–445.
24. É. Bordes, J. Szala-Bilnik and A. A. Pádua, Faraday Discuss., 2018, 206, 61–75.
25. C. J. Shih, S. Lin, M. S. Strano and D. Blankschtein, J. Am. Chem. Soc., 2010, 132, 14638–14648.
26. V. Sresht, A. A. Pádua and D. Blankschtein, ACS Nano, 2015, 9, 8255–8268.
27. G. Salussolia, E. Barbieri, N. M. Pugno and L. Botto, J. Mech. Phys. Solids, 2020, 134, 103764.
28. T. V. Ball and J. A. Neufeld, Phys. Rev. Fluids, 2018, 3, 074101.
29. S. K. Kannam, B. Todd, J. S. Hansen and P. J. Daivis, J. Chem. Phys., 2011, 135, 016313.
30. F. Calabrò, K. Lee and D. Mattia, Appl. Mathematics Lett., 2013, 26, 991–994.
31. E. M. Kotsalis, J. H. Walther and P. Koumoutsakos, Int. J. Multiphase Flow, 2004, 30, 995–1010.
32. M. Neek-Amal, A. Lohrasebi, M. Mousaei, F. Shayeganfar, B. Radha and F. Peeters, Appl. Phys. Lett., 2018, 113, 083101.
33. F.-C. Wang and Y.-P. Zhao, Soft Matter, 2011, 7, 8628–8634.
34. M. Neek-Amal, F. M. Peeters, I. V. Grigorieva and A. K. Geim, ACS Nano, 2016, 10, 3685–3692.
35. H. W. Horn, W. C. Swope, J. W. Pitera, J. D. Madura, T. J. Dick, G. L. Hura and T. Head-Gordon, J. Chem. Phys., 2004, 120, 9665–9678.
36. H. Zhang, H. Ye, Y. Zheng and Z. Zhang, Microfluid. Nanofluid., 2011, 10, 403–414.
37. T. Werder, J. H. Walther, R. Jaffe, T. Halicioglu and P. Koumoutsakos, J. Phys. Chem. B, 2003, 107, 1345–1352.
38. C. Kamal, S. Gravelle and L. Botto, Nat. Commun., 2020, 11, 1–10.
39. A. Maali, T. Cohen-Bouhacina and H. Kellay, Appl. Phys. Lett., 2008, 92, 053101.
40. D. Ortiz-Young, H.-C. Chiu, S. Kim, K. Votchovsky and E. Riedo, Nat. Commun., 2013, 4, 1–6.
41. G. Tocci, L. Joly and A. Michaelides, Nano Lett., 2014, 14, 6872–6877.
42. S. G. Hatzikiriakos, Soft Matter, 2015, 11, 7851–7856.
43. Y. Zhu and S. Granick, Phys. Rev. Lett., 2001, 87, 096105.
44. Y. Zhu and S. Granick, Macromolecules, 2002, 35, 4658–4663.
45. V. Hemadri, V. V. Varade, A. Agrawal and U. Bhandarkar, Phys. Fluids, 2017, 29, 032002.
46. K. Falk, F. Sedlmeier, L. Joly, R. R. Netz and L. Bocquet, Nano Lett., 2010, 10, 4067–4073.
47. R. Lhermerout and S. Perkin, Phys. Rev. Fluids, 2018, 3, 1–13.
48. S. Plimpton, J. Comput. Phys., 1995, 117, 1–19.
49. J. L. Abascal and C. Vega, J. Chem. Phys., 2005, 123, 234505.
50. S. J. Stuart, A. B. Tutein and J. A. Harrison, J. Chem. Phys., 2000, 112, 6472–6486.
51. S. Gravelle, L. Joly, C. Ybert and L. Bocquet, J. Chem. Phys., 2014, 141, 18C526.
52. S. Marchio, S. Meloni, A. Giacomello, C. Valeriani and C. Casciola, J. Chem. Phys., 2018, 148, 064706.
53. W. Ouyang, O. Hod and M. Urbakh, ACS Appl. Mater. Interfaces, 2021, 13, 43533–43539.
54. M. A. González and J. L. Abascal, J. Chem. Phys., 2010, 132, 096101.
55. M. Holz, S. R. Heil and A. Sacco, Phys. Chem. Chem. Phys., 2000, 2, 4740–4742.
56. L. D. Landau and E. M. Lifšic, Theory of elasticity: volume 7, Elsevier, 1986, vol. 7.
57. J. E. Jones, Proc. R. Soc. London, Ser. A, 1924, 106, 463–477.
58. B. T. Kelly, Physics of graphite, Applied Science, 1981, p. 477.
59. R. Okamoto, K. Yamasaki and N. Sasaki, Mater. Chem. Front., 2018, 2, 2098–2103.
60. A. Ghatak, L. Mahadevan and M. K. Chaudhury, Langmuir, 2005, 21, 1277–1281.
61. H. Mosaddeghi, S. Alavi, M. Kowsari and B. Najafi, J. Chem. Phys., 2012, 137, 184703.
62. S. Gravelle, C. Ybert, L. Bocquet and L. Joly, Phys. Rev. E, 2016, 93, 1–7.
63. D. Pihler-Puzović, A. Juel, G. G. Peng, J. R. Lister and M. Heil, J. Fluid Mech., 2015, 784, 487–511.
64. Z.-Q. Wang and E. Detournay, J. Appl. Mech., 2018, 85, 041010.
65. J. R. Lister, J. Fluid Mech., 1990, 210, 263–280.
66. J. Kierzenka and L. F. Shampine, ACM Transactions on Mathematical Software (TOMS), 2001, 27, 299–316.
67. K. Kendall, J. Adhesion, 1975, 7, 55–72.
68. Y. N. Young and H. A. Stone, Phys. Rev. Fluids, 2017, 2, 1–21.
69. J. Obreimoff, Proc. R. Soc. London, Ser. A, 1930, 127, 290–297.
70. L. Botto, Front. Mater., 2019, 6, 302.
71. F. Rieutord, B. Bataillou and H. Moriceau, Phys. Rev. Lett., 2005, 94, 2–5.
72. A. Ghatak and M. K. Chaudhury, Langmuir, 2003, 19, 2621–2631.
73. A. Ghatak, L. Mahadevan, J. Y. Chung, M. K. Chaudhury and V. Shenoy, Proc. R. Soc. London, Ser. A, 2004, 460, 2725–2735.
74. J. N. Israelachvili, Intermolecular and surface forces, Academic press, 2015.
75. S. Gravelle, C. Kamal and L. Botto, J. Chem. Phys., 2020, 152, 104701.
76. J. Wang, D. C. Sorescu, S. Jeon, A. Belianinov, S. V. Kalinin, A. P. Baddorf and P. Maksymovych, Nat. Commun., 2016, 7, 1–7.
77. C. D. van Engers, N. E. Cousens, V. Babenko, J. Britton, B. Zappone, N. Grobert and S. Perkin, Nano Lett., 2017, 17, 3815–3821.
78. K. Novoselov, O. A. Mishchenko, O. A. Carvalho and A. C. Neto, Science, 2016, 353, aac9439.
79. J. Nocedal and S. J. Wright, Nonlinear Equations, Numerical Optimization, Springer New York, 2006, pp. 270–302.
80. Q. Wang, Phys. Lett. A, 2010, 374, 1180–1183.
81. J. W. Kang and S. Lee, Comput. Mater. Sci., 2013, 74, 107–113.
82. D. Dillard, J. Appl. Mech., 1989, 56, 382–386.
83. N. Glassmaker and C. Hui, J. Appl. Phys., 2004, 96, 3429–3434.
84. T. J. Wagner and D. Vella, Soft Matter, 2013, 9, 1025–1030.
85. J. R. Lister, D. J. Skinner and T. M. Large, J. Fluid Mech., 2019, 868, 119–140.
86. C. Dhong and J. Fréchette, J. Appl. Phys., 2017, 121, 044906.
87. D. Pihler-Puzović, G. G. Peng, J. R. Lister, M. Heil and A. Juel, J. Fluid Mech., 2018, 849, 163–191.
88. B. Cross, C. Barraud, C. Picard, L. Léger, F. Restagno and É. Charlaix, Phys. Rev. Fluids, 2018, 3, 062001.
89. A. Münch, J. Phys.: Condens. Matter, 2005, 17, S309.
90. A. Münch, B. Wagner and T. P. Witelski, J. Eng. Math., 2005, 53, 359–383.
91. F. Brochard-Wyart, P.-G. De Gennes, H. Hervert and C. Redon, Langmuir, 1994, 10, 1566–1572.
92. J. Desroches, E. Detournay, B. Lenoach, P. Papanastasiou, J. R.-A. Pearson, M. Thiercelin and A. Cheng, Proc. R. Soc. London, Ser. A, 1994, 447, 39–48.
93. T. S. Chan, C. Kamal, J. H. Snoeijer, J. E. Sprittles and J. Eggers, J. Fluid Mech., 2020, 900, A8.
94. C. Huh and L. Scriven, J. Colloid Interface Sci., 1971, 35, 85–101.
95. Y. Sui, H. Ding and P. D. Spelt, Annu. Rev. Fluid Mech., 2014, 46, 97–119.
96. J.-L. Barrat and L. Bocquet, Phys. Rev. Lett., 1999, 82, 4671.
97. C. Bakli and S. Chakraborty, Nanoscale, 2019, 11, 11254–11261.
98. L. Bocquet and E. Charlaix, Chem. Soc. Rev., 2010, 39, 1073–1095.
99. E. Lauga, M. Brenner and H. Stone, Microfluidics: The no-slip boundary condition, Springer Handbook of Experimental Fluid Mechanics, Springer Berlin Heidelberg, 2007, pp. 1219–1240.
100. H. Hasimoto, J. Phys. Soc. Jpn., 1958, 13, 633–639.
101. S. Gravelle, L. Joly, F. Detcheverry, C. Ybert, C. Cottin-Bizonne and L. Bocquet, Proc. Natl. Acad. Sci. U. S. A., 2013, 110, 16367–16372.
102. C. Belin, L. Joly and F. Detcheverry, Phys. Rev. Fluids, 2016, 1, 054103.
103. J. Happel and H. Brenner, Low Reynolds number hydrodynamics: with special applications to particulate media, Springer Science & Business Media, 2012, vol. 1.
104. K. H. Jensen, A. X. Valente and H. A. Stone, Phys. Fluids, 2014, 26, 052004.
105. O. Bäumchen, L. Marquant, R. Blossey, A. Münch, B. Wagner and K. Jacobs, Phys. Rev. Lett., 2014, 113, 014501.
106. J. Eggers, Rev. Mod. Phys., 1997, 69, 865.
107. S. T. Thoroddsen, H. Nguyen, K. Takehara and T. Etoh, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 82, 046107.
108. J. W. Biddle, R. S. Singh, E. M. Sparano, F. Ricci, M. A. González, C. Valeriani, J. L. Abascal, P. G. Debenedetti, M. A. Anisimov and F. Caupin, J. Chem. Phys., 2017, 146, 034502.
109. M. A. González, C. Valeriani, F. Caupin and J. L. Abascal, J. Chem. Phys., 2016, 145, 054505.
110. S. Bengtsson, K. Ljungberg and J. Vedde, Appl. Phys. Lett., 1996, 69, 3381–3383.
111. G. G. Peng and J. R. Lister, J. Fluid Mech., 2020, 905, A30.
112. B. Roman and J. Bico, J. Phys.: Condens. Matter, 2010, 22, 493101.
113. R. Frisch-Fay, Flexible bars, Butterworths, 1962, p. 220.
114. B. Audoly, Introduction to the Elasticity of Rods, Fluid-Structure Interactions in Low-Reynolds-Number Flows, The Royal Society of Chemistry, 2015, pp. 1–24.
115. A. Lindner and M. Shelley, Elastic Fibers in Flows, Fluid-Structure Interactions in Low-Reynolds-Number Flows, The Royal Society of Chemistry, 2016, pp. 168–192.
116. M. Jabbarzadeh and H. C. Fu, J. Comput. Phys., 2020, 418, 109643.
117. X. Shi, B. Peng, N. M. Pugno and H. Gao, Appl. Phys. Lett., 2012, 100, 191913.
118. Z. Li, R. J. Young, C. Backes, W. Zhao, X. Zhang, A. A. Zhukov, E. Tillotson, A. P. Conlan, F. Ding and S. J. Haigh, et al. , ACS Nano, 2020, 14, 10976–10985.
119. J. Qian, J. Lin, G.-K. Xu, Y. Lin and H. Gao, J. Mech. Phys. Solids, 2017, 101, 197–208.
120. C. Zhao, J. E. Sprittles and D. A. Lockerby, J. Fluid Mech., 2019, 861, R3.
121. Y. Zhang, J. E. Sprittles and D. A. Lockerby, Phys. Rev. E, 2019, 100, 023108.
122. E. Décavé, D. Garrivier, Y. Bréchet, F. Bruckert and B. Fourcade, Phys. Rev. Lett., 2002, 89, 108101.
123. P-G. de Gennes, Adhesion of soft objects, The Physics of Complex Systems (New Advances and Perspectives), IOS Press, 2004, pp. 3–15.
124. H. Perrin, A. Eddi, S. Karpitschka, J. H. Snoeijer and B. Andreotti, Soft Matter, 2019, 15, 770–778.
125. L. Léger, H. Hervet, G. Massey and E. Durliat, J. Phys.: Condens. Matter, 1997, 9, 7719.
126. S. Granick, Y. Zhu and H. Lee, Nat. Mater., 2003, 2, 221–227.
127. T. King, S. Butcher and L. Zalewski, Apocrita - High Performance Computing Cluster for Queen Mary University of London, 2017.

---

Footnotes

† Electronic supplementary information (ESI) available: Technical details of MD simulations, MD results for peeling in vacuum and NMP, and a COMSOL analysis of the entrance pressure drop for circular entrances with slip. See DOI: https://doi.org/10.1039/d1sm01743h

‡ As this sheet has two identical ends, the drag force at each end of the nanosheet (red circles in Fig. 10a) is calculated as half of the required pulling force.

---