A continuum model based on Rayleigh dissipation functions to describe a Coulomb-type constitutive law for internal friction in woven fabrics

A continuum model intended to provide predictions for the response of a woven fabric that includes the effects of friction between fibers is proposed. Specifically, we consider a macroscopic formulation in which the fabric weave is composed of two orthogonal families of continuously distributed yarns. The elastic behavior of the planar fabric is characterized by a second-gradient formulation, incorporating the capacity of the fibers to resist a bending deformation. Particular care is devoted to modeling the action of preventing fiber overlapping through a potential energy barrier. The frictional sliding effect of warp threads interwoven with the weft yarns is introduced through a Rayleigh dissipative function that can be appropriately shaped to consider a Coulomb-type law. Spinning friction of yarns belonging to different families also is conceived when a relative rotation between fibers is present to generalize the dissipation phenomenon involved in the considered sheet. Numerical simulations of the proposed model are provided and discussed.

The concept behind pantographic micro-structures assumes that fibers belonging to different families are allowed to have relative displacements and relative rotations (at the micro-level) and that these displacements and rotations are possible with a vanishing or minimal deformation energy. This concept is not entirely new, as it has been exploited in several technological circumstances, particularly in the woven fabrics used as reinforcements in composites. The possibility of these relative displacements and rotations also in pantographic micro-architectures did motivate the formulation of continuum models in which enriched kinematics is introduced. For instance, in [56,57], two kinematical fields are introduced to describe the placement of each of the two families of fibers constituting the micro-architecture.
Instead, in [37], the kinematics is even richer, as two rotation fields are added. These more detailed kinematical descriptions are strongly motivated by the need to describe experimental evidence and the request to use a "computationally manageable" reduced-order model. In fact, at least in principle, the experimental evidence could have been described by introducing a very detailed 3D Cauchy first gradient continuum modeling the microstructure at the length scale of the constituting fibers' diameter, but the computational burden would have been several millions of degrees of freedom. Such a burden, especially when optimization problems need to be solved [23,33,34], would make any analysis almost unattainable.
In the present paper, we formulate a continuum model designed to give a reduced-order description for the response of a woven fabric whose mechanical behavior has some similarities with that of pantographic sheets. The described phenomenology includes the effects of slipping friction and those involved in the contact among fibers. Refraining, for the moment, from any detailed microscopic analysis, we propose directly a macroscopic model aiming to account for the deformation of a fabric whose weave is composed of two orthogonal families of homogeneously distributed yarns. The imagined planar fabric has a complex mechanical behavior: its elasticity properties will be mathematically modeled by introducing a second gradient deformation energy, which will account for the capacity of the fibers to resist bending deformations. Moreover, the interaction between fibers in contact is also assumed to be elastic with an interaction potential that tends quickly to infinity to block possible fibers' interpenetrations. Additional dissipative effects are incorporated in the model in order to account for: (1) the frictional sliding effect of warp threads interwoven with the weft yarns and (2) spinning friction of yarns belonging to different fibers families. The modeling choice consists in introducing a suitable Rayleigh dissipative function. Such a function is chosen to produce a Coulomb-type friction law among fibers. Spinning friction is introduced in order to model those dissipative phenomena occurring when a relative rotation between fibers is present and plays a relevant role in considered deformation phenomena. More generally, modeling internal dissipation in solids is an age-old problem, and for this reason, a huge amount of literature concerning this subject is available (see [11,12,54,64] for helpful insight on the topic). In these contributions, the dissipated energy often participates in damage and plastic phenomena-possibly by producing cracks or enlarging the existing ones. On the other hand, the dissipated energy could also be associated with friction related to reversible internal relative motions. This last aspect is often underestimated or taken into account through some viscous effects that do not fit the experimental responses of the material since it introduces a rate-dependent behavior not present. Only a few studies are available in this regard in the literature about internal Coulomb-type dissipation in solids [36,52] and more specifically in woven fabrics (see, e.g., [31,43]).
The discussed problem is closely related to the adhesion of thin films and fibers [25,44]; therefore, we believe that the applications of the proposed model can also be extended to these kinds of problems.
The proposed model is validated by performing numerical simulations, which produce interesting predictions about the mechanical behavior of the considered fabric. It has to be remarked that the presented methods may be of use also in the design of pantographic metamaterials: in fact, the considered frictional phenomenology has not been studied yet in this last context apart from one preliminary attempt done in [15].

Modeling
To model a woven fabric, a continuum elastic surface S, lying on a rectangular region Ω belonging to a bi-dimensional Euclidean space, is regarded as a double family of straight orthogonal fibers, i.e., the yarns in warp and weft directions, continuously distributed. More precisely, for the plane surface S and ZAMP A continuum model based on Rayleigh dissipation Page 3 of 17 209 at a given time, a motion, herein regarded as planar, is defined by two placement vector fields, denoted by χ 1 and χ 2 , for the two families of fibers. A practical choice for this kind of engineering fabric is to represent all deformation measures in a reference frame aligned with the fiber directions. For this purpose, the unit tangent vectors to the two families of fibers in the reference configuration, denoted by D 1 and D 2 , are introduced. Consequently, to model the elastic behavior of the textile formed by weaving fibers, all elastic constitutive information is represented through a strain energy that depends on the relative fiber displacement, the first, and second gradients of the placement. A strainenergy function of this type [18,56], incorporating the orthotropic response conferred by the initial fiber geometry, is in which K α e , K s and K α b are material parameters related to the extensional, shear, and bending stiffness, respectively, being α = {1, 2} a label that specifies one of the two families of the fibers. The stiffness K c concerns the elastic interaction between fibers belonging to different families. Specifically, the measures of deformation appearing in the stored strain energy are: 1. the stretching of fibers 2. the shearing distortion: namely the change in the angle between fibers belonging to different families; 3. the Lagrangian curvature of fibers: that is the rate of change of the current tangent vector to the fiber with respect to arc length along the same fiber in the reference configuration; 4. the coupling distortion which represents the interaction term between the two families of fibers. This interaction is due to the geometry of the contact between fibers. The mechanism characterizing such an interaction is analogous to a leaf spring (see Fig. 1). It is worth noting that the energy density related to the shear deformation is written in terms of the tangent of the distortion angle γ [49]; this is a natural way of excluding the overlapping of the fibers since the trigonometric function tends only asymptotically to the values ±π/2 for which it would be necessary to provide an infinite amount of energy to the sample.

Model of dissipation
In the present work, the effects of slipping with friction are considered by means of a macroscopic formulation based on the choice of a suitable dissipation functional which depends on the relative sliding velocities between the two interacting families of fibers [20]. For this purpose, a Rayleigh dissipative potential density is introduced as follows: to model a Coulomb-type friction dissipation and is designed by means of the components of the relative velocityχ 1 −χ 2 . Following this approach, the derivative of such a potential density R with respect to the relative sliding velocities in the directions D 1 and D 2 yields the friction stress τ mutually exerted by the fibers during a sliding motion, whose components are: In the literature, Coulomb friction actions are usually modeled by the introduction of the signum function, whose argument is the velocity of the kinematical quantity on which friction forces act. This function, being not Lipschitz continuous, is a source of strong singularities and numerical or chaotic instabilities when appearing in differential equations. Therefore, the signum function is regularized with a hyperbolic tangent modulated with an amplitude ζ-giving the maximum of friction stress which may be exerted-and with a suitably chosen slope given by η -defining the range of velocities where a viscous behavior is present, namely the friction force is an increasing function of the velocity.  Therefore, the virtual work related to the friction has to be calculated as follows:

Spinning friction
Since the nature of the contact is particularly complex in the considered system, in the present paper, we also explore the possibility of adding a dissipation related to the relative rotation of the fibers that are into contact. Indeed, it is reasonably expected that a relative rotation of the small regions belonging to different fibers in touch produces a share of dissipation in the overall motion due to this kind of contact that is not provided by the model given in (9). For this purpose, a spinning friction term [69] is conceived through a Rayleigh dissipative potential as a function of the relative rotation rateγ: in which ζ spin and η spin play an analogous role of the previous sliding friction contribution. Therefore, the virtual work related to the spinning friction becomes Finally, considering the case in which inertial effects can be neglected and to get the equations of motion, the first variation of the strain-energy functional W has to be equated to the dissipative virtual works, i.e., possibly, at first approximation, the last term can be neglected. In what follows, for the sake of simplicity, we consider only kinematic boundary conditions on a portion of the boundary while other parts are free. Therefore, no external work is taken into account here.

Cyclic responses
Numerical simulations based on the energy (2) with dissipation effects due to frictional slip in woven fabrics (9) are performed using a generalized variational principle including Rayleigh dissipative functions ZAMP A continuum model based on Rayleigh dissipation Page 9 of 17 209 Fig. 9. Comparison between relative displacements for the sample with ratio 1:3 for a soft bending stiffness case (top) and a hard bending stiffness case (bottom) and the finite element software COMSOL Multiphysics. A nonlinear second gradient energy density characterizes the adopted model. Therefore, particular care is used in the formulation, especially in the right choice of interpolating functions that must be able to represent a functional space H 2 appropriately (see, e.g., [16,53]). To analyze the effect of dissipation in a woven fabric, cyclic tensile tests have been carried out for samples having diverse sizes. Specifically, we considered a rectangular specimen with three aspect ratios, 1:2, 1:3, and 1:4, having a short edge of 30 cm long. The two families of fibers are considered with the same geometry and material. Table 1 reports the parameters used for the simulations unless diversely specified. Identification procedures based on digital image correlation can be used to obtain these parameters experimentally (see, e.g., [27,28,45,46]). We remark that, in the case of dynamic simulations with Coulomb-friction, stick-slip phenomena can occur and need special numerical treatment as given by appropriate integrators [70,71]. We executed these tests by clamping one edge of the sample and applying a given displacement in the longitudinal direction being suppressed the orthogonal component. The cycles are characterized by increasing amplitude even though the maximum velocity in each of them remains the same (see Fig. 2). Here, the idea is to appraise the dissipation for different amplitudes keeping the applied maximal velocity constant to examine the fabric nonlinear response. Figure 3 shows the force versus displacement plots of the numerical simulations for the different sizes examined under the same imposed displacement. The nonlinear response is more pronounced for the samples that undergo a large displacement compared with the specimen's overall size. As a matter of fact, for the same displacement, the corresponding value of the force is much larger. A second test has been conducted with a different axial displacement-controlled loading (see Fig. 4). Instead of keeping the maximal velocity constant, the cycle frequency is set unvarying to see the role of an increasing velocity together with the displacement amplitude. Figure 5 exhibits some differences between the two displacement-controlled loadings because, in the model of dissipation, the regularization of the signum function through the hyperbolic tangent implies a certain amount of viscous response that is not provided by the signum, which characterizes a pure Coulomb behavior, namely rate-independent. This modeling choice is quite general since, simply changing the slope of the hyperbolic tangent (η), one can vary the share of viscous behavior involved in the motion.
Subsequently, we investigated the role of the stiffnesses on the hysteretic behavior of the fabric. New tensile tests have been carried out with the constant maximal velocity loading in the case of a sample with an aspect ratio of 1:3. In these tests, we considered diminished and increased stiffnesses to illustrate the response of the specimens changing the overall compliance. The stiffnesses used in these tests are: 1) for the soft case K α e = 8.55 × 10 7 N/m and K α b = 24.05 N m; 2) for the stiff case K α e = 1.25 × 10 8 N/m and K α b = 560.8 N m. Figure 6 displays the comparison between a significantly soft sample and a stiff one. The former having a more considerable deformation shows a remarkable nonlinear behavior, but the dissipation involved remains low since the internal stress is rather small. The latter, on the contrary, being stiffer, shows a less pronounced deformation and consequently a more linear response. At the same time, due to a high-stress level, the dissipation is more significant. Finally, we explored the case in which the previously used model is augmented with the contribution (11) related to the spinning friction. Again, for this case, we compared the two axial displacementcontrolled loadings with constant maximal velocity (using the parameters ζ spin = 30 N/m, η spin = 40 s) and constant frequency (using the parameters ζ spin = 80 N/m, η spin = 50 s) for the rectangular sample with an aspect ratio of 1:3. Figure 7 reports the force versus displacement plots. The effect of the spinning friction is evident since the overall dissipation increases compared with the previous cases. The shapes of the hysteretic cycles also change, especially nearby the points of motion reversal. In a detailed description of the dissipative behavior of these materials, one cannot neglect a priori the influence of this phenomenon. Similar comments can be made about the dependence of the imposed displacement history as done before.

Analysis of the deformation
To have a general picture of the fabric behavior, we showcase the equilibrium shapes of the specimens in the tensile tests in correspondence to the maximum displacement applied. In particular, in what follows, we plot through colors the relative displacement (β), the curvatures of the two families of the fibers (κ α ), and the shear distortion (γ). The first plot is significant to see where the friction phenomenon occurs the most; the second plot represents a quantitative assessment of how meaningful the second-gradient contribution is; the third is mainly related to the spinning friction. The model used in the numerical simulations is continuous. Still, we display only a few significant material lines in the plots to follow the deformation pattern of the considered samples better. We remark that the finite element discretization is unrelated to the shown lines. On that matter, we set the size of the mesh in order to be sufficiently accurate in obtaining the results, making a convergence analysis. Figure 8 shows the pattern of the relative displacement involved in the tensile test for different aspect ratios. The qualitative behavior is almost the same, but naturally, the values of the considered field depend on the relative amplitude of the imposed displacement compared with the specimen sizes. Figure 9 exhibits the effect of the stiffnesses of the fibers. The relative displacement between fibers localizes in the same regions. However, the intensities and the equilibrium shapes of the samples under tension differ. Figures 10, 11, and 12 show the curvatures of the fibers in the tensile test. The regions for which the curvatures are not vanishing are significant with respect to the overall sizes of the specimen; therefore, the stored energy in the bending deformation cannot be negligible if the related stiffness is sufficiently high. This aspect allows us to assess the importance of the second-gradient contribution in the considered energy.
For the sake of completeness, we also plot the distortion angle γ, which is significant for the definition of the spinning friction, in Fig. 13. This quantity is almost uniform over a large central zone of the sample.

Conclusions
In the present paper, we propose a weak formulation based on the principle of virtual work to describe the behavior of a planar woven fabric in the framework of continuous deformable surfaces. The model is characterized by an elastic part that takes into account the elongation, and the bending deformations of two orthogonal families of interwoven fibers thought as continuously distributed over a planar reference region. A further contribution in the energy allows us to model the impenetrability of the fibers when they are close to being aligned. The two families of yarns are described by two different independent placement maps, which interact with each other for the effect of the weave but are allowed to possess a relative slipping. This last relative slipping together with the relative rotation between fibers could be a source of dissipation due to friction. In the paper, a suitable Rayleigh functional is introduced to model these phenomena. The effect of the considered friction dissipation is investigated through numerical simulations.
Funding Open access funding provided by Universitá degli Studi dell'Aquila within the CRUI-CARE Agreement.
Open Access. This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.