## Abstract

The corresponding author had the pleasure of attending an excellent symposium titled “Experimental and Theoretical Micro- and Nano-Mechanics: Honoring the Contributions of Prof. Kyung-Suk Kim” organized by Professors Ashraf Bastawros, Wendy Crone, Yanfei Gao, and Ruike (Renee) Zhao as part of the 2022 Society of Engineering Sciences Annual Technical Meeting held in College Station, TX on October 16–19, 2022. The symposium was held in honor of Prof. Kim’s 70th Birthday and celebrated over 40 years of Prof. Kim’s independent research achievements across several areas of Applied Mechanics. The present paper is dedicated to Prof. Kim, a great colleague at Brown University.

We present a new formulation for the multilayer isogeometric Kirchhoff–Love (KL) shells, where the individual layers are assumed to interact through no-penetration and frictional contact. This work is largely motivated by the experiments and analysis presented in Poincloux et al. (2021, “Bending Response of a Book With Internal Friction,” Phys. Rev. Lett., 126(21), p. 218004). We utilize a regularized version of Coulomb’s friction law to model the tangential traction between the contacting shell surfaces. To ensure objectivity (i.e., reference-frame invariance) in the frictional model, we propose two different strategies to extrapolate the velocity vectors of the contact pair at the contact interface: (i) using the underlying KL kinematics of the individual shell layers and (ii) using the Taylor series-based extension from Kamensky et al. (2019, “Peridynamic Modeling of Frictional Contact,” J. Peridyn. Nonlocal Model., 1(2), pp. 107–121). We compare the performance of both approaches through a numerical benchmark example. We then validate our multilayer shell formulation using the “bending response of a book with internal friction” experiments of Poincloux et al. (2021, “Bending Response of a Book With Internal Friction,” Phys. Rev. Lett., 126(21), p. 218004).

## 1 Introduction

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*^{0}-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.

## 2 Frictional Contact Formulation

In what follows, the subscripts (·)_{1} and (·)_{2} 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 $(\u22c5)^$ normalizes a vector with respect to the *l*^{2} norm denoted by $\Vert \u22c5\Vert $. 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 *ξ*_{3} ∈ [−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].

### 2.1 Normal Contact.

**x**

^{3D}and its variation may be computed using the KL kinematics (see Section 2.1 of Ref. [9]) as

**n**is the mid-surface normal vector computed from the parametric derivatives of the mid-surface position vector

**x**

^{2D}. Here, the “+” and “−” signs indicate the upper and lower shell layer surfaces, respectively, and are chosen such that the two mid-surface material points project onto the same spatial location on Γ

_{c}.

*p*(

*d*) is given by

*k*is the contact stiffness parameter,

*r*is the cutoff distance taken as the local shell thickness average, i.e., $r=(h\Gamma 1+h\Gamma 2)/2$, and

*d*is a scalar-valued gap function given by

*d*> 0.

### 2.2 Tangential Contact.

**v**

^{3D}can be computed by taking a Lagrangian time derivative of the position vector in Eq. (2) to obtain

**B**

^{1}and

**B**

^{2}are the auxiliary matrices detailed in Section 2.1 of Ref. [9]. We use Eq. (10) to compute the variation of the velocity as

_{c}may now be defined as

As the analysis in Ref. [25] demonstrates, if the relative velocity is calculated directly between the two shell mid-surface locations as $v1\u21922=v12D\u2212v22D$, it does not transform like a vector under superposed rigid body motions. As a result, this definition violates the objectivity requirements [26]. By defining the relative velocity as in Eq. (11), objectivity is restored because the relative velocity is computed at the same spatial point on the shared contact surface.

#### 2.2.1 Regularized Coulomb’s Friction Law.

**v**

_{1→2})

_{τ}is the tangential part of the relative velocity vector at the contact interface given by

**I**is the second-order identity tensor and $P=I\u2212n^\u2297n^$.

*R*

_{fr}expressing the frictional constitutive law is defined as

*v*is the relative tangential velocity magnitude, and

*c*

_{fr}is the coefficient of friction that depends on the properties of the materials that come into contact. (Table 5.1 in Ref. [22] includes some exemplary values for different material combinations.) The relationship given by Eq. (17) is a regularized version of Coulomb’s law [24,27,28], where $\epsilon fr$ is a regularization parameter corresponding to the speed at which the stick-like behavior transitions to the full slip behavior. Besides being able to model stick-slip behavior, the regularized version of the friction law leads to improved convergence of the Newton–Raphson iterations in the implicit time stepping procedures employed in the present work. Figure 1 shows how the normalized friction magnitude varies with the relative slip speed for a given choice of the constitutive law parameters.

### 2.3 Linearization.

*R*

_{fr}for brevity, the linearization of the traction may be expressed as

*R*

_{fr}in Eq. (19) may be computed as

**C**

_{fr}may be thought of as a tangent damping tensor, which is symmetric and given by

## 3 Numerical Results

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*^{1}-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.

### 3.1 Sliding Panel.

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^{7} N/mm^{2} 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^{3} 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^{3}.

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.

*ρ*is the mass density, Ω is the elastic panel domain, and

**e**

_{1}and

**e**

_{2}are the unit normal vectors in the horizontal (along

*f*

_{h}) and vertical (along

*f*

_{v}) directions, 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.

### 3.2 Bending Response of a Paper Stack With Friction.

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 $\epsilon 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.

## 4 Conclusion

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.

## Acknowledgment

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.

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

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