Multilayer Shells Interacting Through Friction

Contact-impact treatment of shell-to-shell interactions with friction has a broad spectrum of applications. One traditional engineering scenario that involves frictional contact is car crash, where the structural members self-contact, buckle, and slide past each other, exhibiting finite deformations and dissipating impact energy. On the interface modeling side, more recently, the mechanical response of layered architected structures has gained much attention due to their numerous applications in biology [1,2], thin films [3,4], and multilayer graphene [5–7]. In all these applications, interlayer friction plays an essential role in defining the global response of multilayer structural systems. Although analytical and experimental methods to study the response of multilayer structures comprising of thin elastic rods and plates interacting through friction were recently proposed [8,29], there is still a demand for accurate and efficient computational methodologies of general purpose and applicability in this area. Motivated by this demand, we develop a computational framework to model the response of multilayer structures governed by interlayer friction. Our approach models each layer as a large-deformation thin shell and defines frictional contact interaction between the layers. At the level of each layer, we employ an updated Lagrangian, isogeometric Kirchhoff–Love (KL) shell formulation, building on the authors’ recent work in Ref. [9].

We note that over the past decade, isogeometric analysis (IGA) [9–14] has matured to be a computationally efficient, accurate methodology for the modeling of thin shells and deployed in numerous applications like fluid–structure interaction [15], failure analysis of carbon-fiber reinforced composites [9], curvature-resisting material surfaces [16], bioprosthetic heart-valve analysis [17], and car crash [14]. IGA has also proven to be advantageous in the modeling of shell-to-shell contact [18] over C⁰-continuous Lagrange polynomial discretizations. The ability to exactly represent shell mid-surfaces through smooth B-splines directly translates into better accuracy and smoothness in representing the contact tractions, especially for frictional contact with sliding [19,20].

Because we aim to simulate structural systems involving many layers, we develop and employ a simple and efficient penalty-based approach to enforce normal and tangential contact constraints. Although more sophisticated contact-constraint enforcement techniques such as Lagrange multiplier or mortar methods may be used, they typically suffer from a high computational cost due to the extra unknowns employed [21,22], and often require complex implementation [23].

The remainder of this paper is structured as follows. In Sec. 2, we formulate the normal and tangential parts of the contact model, discuss the two options of extrapolating the velocity vector at the contact interface, present a regularized Coulomb’s friction law [24], and show the linearization of the frictional contact terms that is needed for computer implementation of implicit time stepping. In Sec. 3, we present the numerical results focusing on experimental validation. In Sec. 4, we draw conclusions and discuss future research directions.

In what follows, the subscripts (·)₁ and (·)₂ indicate a quantity on sides 1 and 2, respectively, of a shared contact interface denoted by (Γ_c). The superscript (·)^3D denotes an extrapolation of a shell mid-surface quantity, labeled (·)^2D, to the contact interface via the KL kinematics [9,12,13]. We denote by h_Γ the thickness of a single shell layer. The operator $(⋅)^$ normalizes a vector with respect to the l² norm denoted by $‖⋅‖$⁠. The notation δ(·) stands for the variation with respect to the mid-surface unknowns, ξ_i=1,2 indicates the shell in-plane parametric coordinate directions, and ξ₃ ∈ [−1, +1] corresponds to the shell through-thickness parametric direction. In what follows, we focus exclusively on the contact formulation. The detailed formulation of the isogeometric KL shell that describes the behavior of the individual layers may be found in Ref. [9].

We compute two examples: (i) a relatively simple numerical benchmark problem from Ref. [25] where we assess the performance of the KL-based versus Taylor expansion-based velocity extrapolation and (ii) a more involved three-point bending test of a paper stack, where we are able to compare our results with the experimental measurements of Ref. [29]. For all the calculations we use C¹-continuous quadratic non-uniform rational B-spline (NURBS) and employ three-point Gaussian quadrature in the plane and two-point Gaussian quadrature through the thickness in each layer.

Following Ref. [25], we simulate a sliding action of an elastic panel on a rigid flat surface. The problem configuration is depicted in Fig. 2. The rectangular rigid surface measures 10 mm × 2 mm and has a constant thickness of 0.0625 mm. The deformable elastic square panel has an edge length of 1 mm and a constant thickness of 0.125 mm. The elastic panel is assumed to be made of an isotropic material with Young’s modulus E = 2.5 × 10⁷ N/mm² and Poisson’s ratio ν = 0.25. The flat surfaces come into contact after applying a uniformly distributed vertical body force per unit volume of f_v = 1.0 N/mm³ to the elastic panel. After static equilibrium is achieved, we induce friction by sliding the elastic panel on the rigid surface using a constant, uniformly distributed horizontal body force per unit volume of f_h = 1.0 N/mm³.

We compute the problem using two meshes, denoted by M1 and M2, and two different values of the friction coefficient c_fr = 0.2 and c_fr = 0.5. We also compare the performance of the KL kinematics-based and Taylor expansion-based velocity extrapolation. The total number of elements for M1 and M2 are 100 × 20 + 10 × 10 = 2,100 and 200 × 40 + 20 × 20 = 8,400, respectively.

In Figs. 3–6, we compare the time history of the effective friction coefficient $cfreff$ to the constant value c_fr employed in the regularized Coulomb’s friction law (Eq. (17)). For both values of c_fr, the Taylor series-based extension results in $cfreff$ fluctuate around the analytical value. The fluctuations arise due to the fact that locally the relative friction velocity may have a small out-of-plane component. The fluctuation levels reduce appreciably with mesh refinement. These results are consistent with the original reference on the Taylor expansion-based approach in Ref. [25]. In contrast to the Taylor expansion-based case, the KL extension gives an accurate prediction of the effective friction coefficient and shows no fluctuations, even on the coarse mesh. With this result, in the next section, we simulate a more complicated scenario using only the KL kinematics-based extension version of the formulation.

Inspired by experiments from Ref. [29], we present a series of three-point bending test calculations of a 40-layer paper stack. A rigid indenter that drives the transverse deflection is placed halfway between the supporting rollers, which have a radius of 6.8 mm and are positioned 65 mm from the center the stack. Each elastic plate measures 220 mm in the longitudinal direction and 30 mm in the width direction and has a constant thickness of 0.286 mm. (See Fig. 7 for the problem configuration.) The stack layers are made of polyethylene terephthalate. The material parameters for each layer are: Young’s modulus E = 2.4 GPa and Poisson’s ratio = ν = 0.44. In the experiments, the supporting rigid rollers are coated with a thin film of vinyl polysiloxane to mimic frictionless contact between the bottom layer of the stack and the rollers. To provide dry-friction conditions and avoid interlayer adhesion, the surface roughness of each layer is increased by applying sandpaper.

For the discretization of each layer in the stack, we use 353 elements in the longitudinal direction and a single element in the width direction. In a similar fashion, to discretize the frictionless roller supports, we employ 804 elements in the circumferential direction with a single element in the width direction. The problem mesh consists of 42 NURBS patches with 15,728 elements in total. In the computations, we assume that the layers are airtight. That is, there are no initial gaps present between the layers in the stack.

The stack is displaced at the midline of the top layer in the transverse direction up to 19.5 mm and then unloaded to the undeformed configuration. The stack is displaced at a constant speed of 1 mm/s with a short ramp-up and ramp-down of the velocity at the start, end, and reversal of the loading. We employ several values of the coefficient of friction to demonstrate the influence of interlayer friction on the bending response of the stack.

We carry out several computations of the 40-layer stack and compare the predicted force–displacement response with the experimental results in Fig. 8. (See Fig. S1(b) in the supplementary material of the experimental reference [29].) We first perform a computation without interlayer friction (i.e., we only enforce the no-interpenetration condition and assume the frictional contact contributions are identically zero). We compare the results with those for a single-layer computation and a force response amplified by a factor 40. The two force-deflection curves are nearly indistinguishable and provide a lower-bound for the stack bending stiffness, which scales with the number of layers for the frictionless case. (See also the analysis and discussion in Ref. [29].)

We then carry out two additional simulations, one employing c_fr = 0.50 and another c_fr = 0.75. The same regularization parameter $εfr=0.001mm/s$ is employed in both cases. In Figs. 9(a)–9(c), we visualize the horizontal component of the tangential traction due to friction at the three representative stages of the deformation process. The figure indicates that throughout the deformation, the tractions cluster below the indenter and above the roller supports. In addition, the tractions appear to be smoothly distributed along the thickness of the stack despite the discontinuity in the kinematics between the layers.

Introducing friction increases the bending stiffness during the loading part of the curve. This is because friction in the stack creates additional forces that counteract the downward force produced by the indenter (see Fig. 9(a)). During unloading, however, the sign of the traction vector due to friction reverses and creates a contribution that reduces the indenter force (see Figs. 9(b) and 9(c)). This creates a hysteresis behavior that is very well captured in the experiments as well as in our simulations. We also note that the computational results are self-consistent in that the loading part of the force–displacement curve is above, and the unloading part is below, the data for the frictionless case.

Compared to the experimental findings in Ref. [29], we have a fairly accurate prediction of the peak load and hysteresis curve. For c_fr = 0.75, the peak load in the computations is only underestimated by $3.5%$⁠. On the other hand, for c_fr = 0.5, we obtain a better estimation of the hysteresis compared to the experiments. The discrepancies in the peak load and the amount of dissipated energy may be attributable to the inherent minor differences between the modeling assumptions and physical reality. For instance, in the early stages of the experiment for up to 1 mm of deflection, the reaction forces reported appear to be near zero. This may be explained by the presence of small initial gaps between the layers in the experimental setup. Because this is not reflected in the modeling, the computations produce non-zero reaction forces as soon as the indenter comes into contact with the top layer.

During the unloading stage, as the interlayer friction reverses sign, we also observe a slight separation and dynamic “chattering” of layers in the bottom of the stack, which is consistent with the experimental observations reported in Ref. [30]. This finer-scale dynamic phenomenon can only be captured using a multilayer modeling approach developed in this work.

We developed and validated a new multilayer isogeometric KL shell formulation, where the individual layers are assumed to interact through the enforcement of no-penetration and frictional contact conditions. The latter makes use of a regularized Coulomb’s law, which is able to model stick-slip behavior and help improve the nonlinear convergence of the implicit time integration procedure. We compared two strategies for extrapolating the velocity vector in the frictional contact formulation and found that the approach based on directly using the KL shell kinematics delivers the best results. Finally, we carried out a set of computations for a 40-layer stack subjected to bending and obtained very good results in being able to reproduce the load–displacement curves and hysteresis behavior reported in the experiments. We also note that the distribution of the traction due to friction is smooth across the thickness of the stack, which is attributable to the accurate approximation of the KL shell kinematics by the IGA approach.

The present work clearly demonstrates the ability of the proposed approach to handle frictional contact scenarios involving a large number of layers with finite sliding. The formulation presented opens the door for problems where shell-to-shell frictional contact plays an important role, e.g., extreme-event mitigation using architected structures [31] or heart-valve analysis [32]. In addition, the present approach may be further extended to model the behavior of multilayer structural systems where both friction and van der Waals interactions play an important role.

The authors were supported through the ONR Grant No. N00014-21-1-2670. We gratefully acknowledge the high-performance computing resources provided by the Center for Computation and Visualization (CCV) of Brown University.

There are no conflicts of interest.

The authors attest that all data for this study are included in the paper.