1 Introduction

An important modeling issue in the design of metamaterials micro-architecture concerns the need of describing those phenomena occurring when parts of the afore-mentioned microstructure, initially distant, enter in contact because of deformation. This kind of phenomena occurs at micro-level especially in lightweight metamaterials when their microstructure occupies a small percentage of volume fraction and includes parts undergoing large relative displacements. Starting from the basics outlined in [1,2,3], we want to present a first approach that intends to give some research perspectives in the study of contact phenomena occurring in the micro-architecture of fibrous metamaterials. To this end, we will follow in the presentation of the results already available in the literature a criterion of simplicity so that in the following we can introduce more easily the new concepts that are discussed here for the first time.

For seek of simplicity, we will start by addressing the problem of modeling contact interactions in a particular class of complex fibrous metamaterials, namely the pantographic architectured metamaterial. As it will be clear in the following, the topological arrangement of fibers in a metamaterial is of fundamental importance. Interesting results on this topic can be found in [4, 5].

First of all, before we can define what a pantographic microstructure or metamaterial is (and there is a remarkable length-scale difference between these concepts), it is necessary to clarify the terminology used in the present paper. These terminologies are necessary to underline the role of the different scales that play a relevant role in the modeling process and that account for the hierarchical structural ordering that are the main reason of most of the problems that we will present (as in any problem of metamaterial synthesis, the hierarchy plays a fundamental role).

Fig. 1
figure 1

Geometry of a pantographic structure (a). Top view (b) and side view (c) of a 3D-printed real sample

Full size image

The first question that needs to be answered if we want to deal effectively with the complexity of microstructural elements as the hinges, also known as pivots in the literature, is: which is the microstructure of a pantographic metamaterial? We can answer this question by describing the geometrical properties of the object or by showing a real sample like the one in Fig. 1b and c. From a geometrical point of view, as it is also clear observing Fig. 1a, a pantographic structure is constituted by two fiber layers that are on two parallel planes held together by mechanical hinges. The directions of the fibers, in the initial configuration, are orthogonal, and the hinges are located in correspondence of the points of intersection between the fiber directions. But to fully characterize the topic, it is better to adopt the typical approach of those who want to conceive a new metamaterial. This approach consists in specifying a priori the equations governing the metamaterial and from these to find backwards the microstructure that, once homogenized, produces a continuum governed by these equations. In this context, the real question to be answered to describe the pantographic metamaterial is: what are the mechanical properties of this metamaterial and which equations govern it?

To answer this question, we can refer to [3]. In fact, in [3] the basic characteristics of the pantographic metamaterial are outlined. We can distinguish two basic aspects: the primary request is technological and implies, producing a secondary request, a property related to mathematical modeling. The primary demand is to obtain a material that can be stretched without any energy expense (or with a very low expense of it). The mathematical implication, as we will see below, is that the governing equations can be obtained considering a higher gradient continuum, specifically a second gradient one.

In this context, the role of the pivots, we will see, is merely to ensure that the two fiber layers that make up the pantographic metamaterial can exhibit relative rotations at their points of contact. Section 2 deals with this problem. The role of the pivots can be generalized, and one can use them, for example, as actuators of deformation in case they are involved in flexoelectric or piezoelectric processes. This has been attempted in [6, 7].

2 The role of perfect pivots in pantographic micro-architectures

The earliest conception of a pantographic structure as a micro-model coupled with a second gradient macroscopic continuum model occurs in [1, 2]. The key insight consisted in finding the micro-model that would lead, in the sense of some kind of homogenization procedure, to the simplest model characterized by a second gradient continuum. In fact, it was shown there how looks like a microstructure for which elongations require no strain energy and gradient of elongations are associated with a strain energy.

The above-mentioned procedure can be qualified as a multiscale approach with levels organized according to a hierarchical classification. Mechanics provides several examples of multiscale approaches, mainly designed to define the interrelation between macro-models and micro-models. Highly important examples include those due to Maxwell [8] and Saint-Venant [9,10,11]. In the domain of multiscale approaches, a high efficiency approach involves asymptotic identification. When the micro- and macro-models have been postulated a priori, a kinematic matching among them can be conjectured and, successively, equivalence between the powers in the corresponding motions (micro and macro) is applied. If the kinematic matching has been sufficiently judicious, the range of applicability of macro-model covers a large class of micro-dynamic processes. In practice, the kinematical matching determines the range of applicability of calculated continuum limit. A pristine discussion of this point is due to Piola [12,13,14,15,16] so that, sometimes, this micro-macro identification is called “Piola identification.”

Using Piola identification, we can estimate the constitutive parameters of the macro-model equations with respect to the parameters of the fundamental micro-architecture elements. Besides, before presenting some peculiar features of the model used in [1, 2], we would like to underline once more the main characteristic of Piola identification. In it, at first we postulate the macro-model as a second gradient continuum, and only then we look for a possible micro-model that produces the specified macro-model through homogenization. Such macro-theory-driven strategy has a very strong potential: we do not search for random microstructures trying to identify one that is adequate for our objectives, but instead, keeping the macro-theory in mind, we aim to design the appropriate microstructure in order to comply with the theoretical expectations.

When pantographic structures have been introduced in the conception of new metamaterials, second gradient models were already available in the literature [17,18,19,20,21,22,23,24,25,26,27,28,29,30]. In this regard, we can refer to the theories of elastica studied by Euler, Bernoulli and Navier as the very first example of a second gradient model. The model presented by Euler is a 1D prototype. Almost a century later, it was proposed the first model of microstructured continuum by the Cosserat brothers [31]. As shown in [32], when a suitable constraint is introduced, Cosserat continua can be regarded as a generalization of an incomplete second gradient continuum model.

It is, however, noteworthy that the genesis of the higher gradient 3D models—and also of the peridynamic generalized continua—can be attributed to Piola [12, 13] or maybe to earlier authors (Lagrange, D’Alembert).

With Germain [22,23,24, 33], we refer to second gradient continua those media whose Lagrangian volumetric deformation energy depends both on the first and second gradient of the placement field (of course, these dependencies must be independent from the observer). We refer to an incomplete linear second gradient material if its strain energy only depends on \(\nabla u\) and \(\nabla \omega (u)\), where \(\omega (u)\) is the skew-symmetric part of the gradient of the displacement, \(\nabla u\), \(\omega (u)=\nabla u-\varepsilon (u)\) with \(\varepsilon (u)\) the symmetric part of \(\nabla u\). Complete 2D and 3D second gradient models can be found in the capillarity models or even in the damage and plasticity theory [34,35,36,37,38,39,40,41,42].

2.1 Some microstructure can be described by a second gradient model

In the sixties and seventies, the problem of a microstructured medium was approached by Mindlin, Toupin and Germain in several papers. In their research, they discussed different perspectives of higher order models to treat materials with a microstructure, and pointed out how the mere existence of the microstructure, in some cases, can induce higher order terms in the equilibrium equations for the material under investigation. Unlike classical homogenization techniques, in this case equations are obtained that contain terms dependent on second or higher order derivatives of displacement, thus inducing higher gradient theories.

According to Germain’s work [22,23,24], we can demonstrate how the kinematics due to the presence of the microstructure generates a second gradient continuum at macroscopic level. In the classical representation, a continuum is composed of a continuous distribution of particles, geometrically represented by a point M and characterized by a velocity field, defined by its \(U_i\) components. However, in a theory that takes into account the presence of a microstructure, each particle represented by a point M must be characterized by a more sophisticated kinematics. So, how can we describe the presence of the microstructure from a kinematic point of view? For this purpose, it is necessary to consider the continuum at microscopic level: each particle must be considered as a continuum of small extension itself. This small continuum replacing the adimensional point is addressed by Germain as \(\mathcal {P}(M)\). Germain [22,23,24] demonstrates in detail how the previous hypothesis implies that the velocity field \(U_i\) associated with the continuum \(\mathcal {P}(M)\) and the field \(\varphi _{ij}\) of relative velocity gradients (resulting in a second-order tensor) must be considered. This ultimate result makes it necessary to introduce the second relative velocity gradient \(v_{ijk}\), which is a third-order tensor. Therefore, in [22,23,24] it is clearly proven that the presence of a microstructure can be naturally described by including higher order terms in the considered continuum theory. This is the starting point in the study of pantographic structures. A certain microstructure is chosen to have a second gradient continuum as simple as possible and then homogenization techniques are used to determine an appropriate continuum model.

A clear and elegant comparative analysis of the advantages obtained using the principle of virtual work (or the principle of minimum action) can be found in Hellinger [43, 44].

2.2 The pantographic microstructure

The pantographic architecture, see a planar 2D array design in Fig. 1a and an example of its practical realization in Fig. 1b and c, has been first introduced in the context of homogenized generalized media in [45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60]. In order to achieve the associated homogenized macro-model, it is required to examine a structure consisting of n pantographic elements and to analyze the behavior of the limit material when n tends to infinity.

Certain formal asymptotic expansion techniques, formerly employed in [3, 61,62,63,64], are consistently explored in [65,66,67,68] for establishing the effective properties of periodically assembled linear elastic rods-based structures. Noticeably, when the flexural stiffness of homogeneous elastic bars forming the pantographic arrays and torsional stiffnesses of interconnecting cylinders are lower than the beams elongation rigidity [3], one achieves attractive macro-models. In searching for macro-models for micro-architectural metamaterials, there is generally a complex evaluation to be carried out, i.e., we have to establish if the energy associated with the second gradient terms is negligible with respect to the first gradient term. Longly, it has been assumed that the second gradient energetic terms were always irrelevant [16]. This assumption was demonstrated to be inadequate at the very beginning of the modern continuum mechanics by Gabrio Piola [12, 13], albeit some scholars have longly ignored the opinion of Gabrio Piola in this context while praising some of his results [51].

In order to prove that the synthesis of second gradient metamaterials with non-negligible second gradient energetic terms is not only physically conceivable, but can be approached mathematically, it has been proven that [1, 2]:

  1. i.

    pantographic architectures provide for the synthesis of Casal-type beams [69, 70];

  2. ii.

    it is possible to synthesize second gradient plates with two families of pantographic substructures: i.e., plates having strain energy depending on the in-plane second gradients of displacement [71];

  3. iii.

    when perfect hinges are provided, the macro-strain energy of short beam pantographic structures does not incorporate at all the shear first gradient terms, and in the pantographic fabrics the so-called floppy-modes can be observed: i.e., homogeneous deformations that correspond to a null strain energy.

2.3 A second gradient continuum proto-model for the pantographic metamaterial

Consider, now, a reference configuration of a pantographic lattice (see Fig. 1a). As it should be now clear, it is constituted by two families of mutually orthogonal fibers which intersect one by one in a bidimensional array of regularly distributed material points. These points are modeling the hinges or pivots we have already introduced. Let be \((D_1,D_2 )\) an orthogonal basis for the reference configuration. For simplicity, \((D_1,D_2 )\) can be taken along the directions of the two fiber families which are, in the reference configuration, orthogonal. So, we will refer to the reference configurations of the two families of fibers by using the two unitary vectors \(D_{\alpha }\) with \(\alpha =1,2\).

We now consider a 2D continuum whose reference shape is given by a rectangular domain \(\varOmega = \left[ 0,\mathrm {L}\right] \times \left[ 0,\ell \right] \subset \mathbb {R}^2\), where \(\mathrm {L}\) and \(\ell \) can be identified as the lengths of the sides of the pantographic structure in Fig. 1a and very often, it is assumed that \(\mathrm {L}=3\ell \). If we want to study only planar motions, then the current shape of the rectangle is mathematically described by a regular placement function \(\varvec{\chi }: \varOmega \rightarrow \mathbb {R}^2\). In [3], it is shown how, by employing the so-called Piola’s Ansatz and assuming that \(\varvec{\chi }( \cdot )\) is at least twice differentiable, one obtains the following macroscopic homogenized strain energy for the pantographic metamaterial

$$\begin{aligned} \mathcal {U} (\varvec{\chi } ( \cdot ))= & {} \int _{\varOmega }\sum _{\alpha }\frac{\mathrm {K}_{e}^{\alpha }}{2} (||\varvec{F} D_{\alpha }|| -1)^{2}\; \mathrm{d} \varOmega \nonumber \\&+ \int _{\varOmega } \sum _{\alpha }\frac{\mathrm {K}_{b}^{\alpha }}{2}\left[ \frac{\nabla \varvec{F} \vert D_{\alpha } \otimes D_{\alpha } \cdot \nabla \varvec{F}\vert D_{\alpha } \otimes D_{\alpha }}{||\varvec{F}D_{\alpha }||^{2}}\right. \nonumber \\&- \left. \left( \frac{\varvec{F} D_{\alpha }}{||\varvec{F} D_{\alpha }||} \cdot \frac{ \nabla \varvec{F}\vert D_{\alpha } \otimes D_{\alpha }}{||\varvec{F} D_{\alpha }||}\right) ^{2}\right] \; \mathrm{d} \varOmega \end{aligned}$$
(1)

where \(\varvec{F}\) denotes the deformation gradient \(\nabla \varvec{\chi }\) and \(\left( \nabla \varvec{F}\vert D_{\alpha } \otimes D_{\alpha }\right) ^{\beta }=F^{\beta }_{\alpha ,\alpha }=\chi ^{\beta }_{,\alpha \alpha }\). The first integral in Eq. 1 models the fiber elongation energy, while the second integral represents the fiber bending energy. We want to remark that the bending energy is written in terms of the gradient of the deformation tensor \(\nabla \varvec{F}\) corresponding to a second gradient of the placement function \(\nabla ^2\varvec{\chi }\). This is coherent with what we had previously stated, as the second gradient term is representative of the bending of the fibers.

2.4 A mechanical diode

As we have remarked, the continuum energy in Eq. 1 takes into account only bending and elongation of fibers: no connective microstructural elements deformation is considered. If this is the case, one can recall as an analogy an electric circuit component: the diode.

In many applications of technological relevance, in fact, the voltage-current response of an ideal diode, under static conditions, can be approximated by a piecewise linear function. In this approximation, the current can be considered null if the voltage between anode and cathode is less than or equal to a certain threshold value \(V_T\) ; if, on the contrary, the voltage is higher, the diode can be approximated to a voltage generator, whose current is imposed by the circuit to which it is related. In mechanics, pantographic metamaterial exhibits a behavior formally identical to the one exhibited by the diode in electrical circuits (this analogy has been studied in detail in [72, 73]). As already explained in previous sections, the macroscopic deformation of this metamaterial results, at microscopic level, in the deformation of its basic constituents, i.e., fibers and pivots. If the pivots are perfect and their torsional stiffness is therefore equal to zero or negligible, so the overall deformation of the structure results in the sole deformation of the fibers.

Fig. 2
figure 2

Particular of a bias extension test. In (a), two bent fibers are evidenced in red. In (b), the same two red fibers enter in contact (color figure online)

Full size image

A bias extension test of a pantographic structure (with perfect pivots or not, see Fig. 2) consists in two successive stages: (i) the bending of the fibers (coupled to an extremely negligible elongation of the fibers), and then (ii) an elongation of fibers after the fibers in the central part of the structure have come into contact. Moreover, in the second stage of deformation, when the fibers are in contact, one should take into account also this fiber contact. We will deal with the modeling of the fiber contact interaction in the following sections.

When the connective microstructural elements are not perfect hinges, then Eq. 1 has to be completed by adding a shear first gradient term. We will deal with this further term in the next Section.

3 Deformation mechanisms of the connecting microstructural elements

When the microstructural connections (we have called them ‘hinges’ or ‘pivots’) cannot be considered as perfect hinges, one has to deal with their deformation mechanisms and it is necessary to add two further terms to the strain energy in Eq. 1 .

As it is shown in Fig. 3, if we model the microstructural connections as elastic cylinders, one must take into account the torsion and the shear of the cylinders. The torsion of the cylinders is, from a macroscopic point of view, related to the (macro) shear deformation, and it will provide for a new energy term. At this point, it is fundamental to underline that the shape of this term depends heavily on the microscopic mechanics of the cylinder used to mimic the hinges. So, we can have a perfect hinge, an elastic hinge, a plastic hinge and all the intermediate possibilities.

Fig. 3
figure 3

Deformations of a cylinder used to model the pivot: it can be sheared (resulting in a relative sliding of two fibers) and it can be torqued (resulting in a relative rotation between two fibers)

Full size image

A simple form of this energy term, which is adjustable with successful results in identification procedures with both elastic and plastic hinges, has been assumed (as shown in [3]) to depend on the angle between the interconnected fibers from the pivot to a certain power \(\beta \): the latter parameter can be obtained from the fit of the experimental data and in general can also depend on the type of material by which the structure is manufactured. So, the proposed form for the shear energy is

$$\begin{aligned} \mathfrak {U}_s (\varvec{\chi } ( \cdot ))= & {} \int _{\varOmega }\frac{\mathrm {K}_{s}}{2} \left| \arccos \left( \frac{\varvec{F} D_{1}}{||\varvec{F} D_{1}||} \cdot \frac{\varvec{F} D_{2}}{||\varvec{F} D_{2}||}\right) -\frac{\pi }{2}\right| ^{\beta } \; \mathrm{d}\varOmega \end{aligned}$$
(2)

where \(\mathrm {K}_{s}\) is the shear stiffness and can be related to the torsional stiffness of the cylinders modeling the connective microstructural elements and the function \( \arccos \left( \frac{\varvec{F} D_{1}}{||\varvec{F} D_{1}||} \cdot \frac{\varvec{F} D_{2}}{||\varvec{F} D_{2}||}\right) -\frac{\pi }{2}\) represents the variation of angle between two fibers in correspondence of a chosen pivot.

Obviously, the shear energy can also be modeled in other ways. For example, a model that captures more accurately the phenomenology of a torqued cylinder is due to Ogden [74]. A version of Ogden’s shear energy adapted to pantographic structures is the one shown in [75]

$$\begin{aligned} \mathfrak {U}_{Ogd} (\varvec{\chi } ( \cdot ))= & {} \mathrm {K}_{s_1}\left[ \left( 1+\left( \frac{\vartheta }{\vartheta _0}\right) ^2\right) ^{\beta }-1\right] \nonumber \\&-\mathrm {K}_{s_2}\left( \log (\mathcal {J}-\mathcal {J}_0)-(\log (1-\mathcal {J}_0)-\mathcal {J}+1)\right) \end{aligned}$$
(3)

where \(\vartheta \) is the variation of the angle between the fibers in correspondence of the interconnecting pivot and \(\mathcal {J}\) is related to the deteminant of \(\varvec{F}\), while \(\vartheta _0\) and \(\mathcal {J}_0\) are the values of these entities in some particular points of the deformation history. These values can be obtained by fitting the experimental data.

Some problems related to the shear energy are discussed in [76]. An attempt to find a further form for shear energy, able to well reproduce the observed phenomenology, is being performed recently.

The other deformation mechanism which has to be taken into account for properly modeling the connective elements is related to the (micro) shear deformation of the cylinders (see Fig. 3). From a macroscopic point of view, it is responsible for the relative displacement between the fibers of the two families which are connected by the considered pivot. From a microscopic point of view, it is pretty clear that this mechanism is strongly related to the ratio between the height and the radius of the cylinder. The slender the cylinder, the more visible this (micro) shear effect will be. The resulting energy term will be analyzed in the next Section.

4 Interactions between fibers due to relative displacements

4.1 Sliding of fibers of one family with respect to the other one

In Sect. 2, it has been presented the strain energy for the pantographic structure and its principal implications. The main assumption consists in the fact that the two families of fibers are described by using a single placement function \(\varvec{\chi }\).

Fig. 4
figure 4

Sliding of two fibers. The pivot is not modeled as a rigid connection (as it was in previous models), but as a linear spring

Full size image

As it has been shortly discussed in the previous Section, for properly taking into account the whole deformation mechanisms of the connective microstructural elements, it is necessary to introduce and analyze the possibility to allow to the fibers also to slide one with respect to the other. To this aim, we have to modify the homogenized strain energy for the system as presented in Eq. 1 by using two independent placement functions \(\varvec{\chi }^{\alpha }\) (each for a fiber family) and introducing an interaction term in the energy [77]

$$\begin{aligned} \mathfrak {U}_{fc}(\varvec{\chi } ( \cdot ))= \int _{\varOmega }\frac{\mathrm {K}_{fc}}{2} ||\varvec{\chi }^{1} -\varvec{\chi }^{2}||^{2} \; \mathrm{d}\varOmega \end{aligned}$$
(4)

We could imagine also a different interaction energy between the fiber families, for example by adding some suitably conceived extensional springs linking the pivots of different fibers. In this case, we would have a more complex interaction energy, depending on the first gradient of displacement

$$\begin{aligned} \mathfrak {U}^{(\nabla )}_{fc}(\varvec{\chi } ( \cdot ))= \int _{\varOmega }\frac{\mathrm {K}^{(\nabla )}_{fc}}{2} || \nabla \varvec{\chi }^{1} - \nabla \varvec{\chi }^{2}||^{2} \; \mathrm{d}\varOmega \end{aligned}$$
(5)

The energy term that we have now presented does not complete the discussion about the interaction between the fibers. In fact, we are completely neglecting the interaction between the fibers of the same family (at this level, nothing prevents the fibers to interpenetrate each other when we deform the structure). We will have to introduce then a new term which takes into account this kind of interaction and our simplifying assumption is to write it using only the two placement fields \(\varvec{\chi }^{\alpha }\).

4.2 Contact interaction between fibers of the same family

We want now to account for the interactions between fibers of the same family. As it is well explained in Fig. 5, depending on the direction of the relative motion between two parallel fibers, we expect to have two possible interactions: frictional interaction (Fig. 5a) and contact interaction (Fig. 5b). The friction interaction can be mathematically modeled by exploiting a Rayleigh potential, and we will deal with it at the end of this Section. The contact interaction can be dealt with as it is described in the following.

Fig. 5
figure 5

Frictional interaction (a) and contact interaction (b) between two fibers of the same family

Full size image

Without lack of generality, consider now just one of the two families of fibers (as shown in Fig. 6) and we will neglect in what follows the family index \(\alpha \) in the placement function. In the contact energy term, we have to introduce the micro-level constraint of impenetrability between two fibers of the same family. Let consider two adjacent fibers, for example in the family oriented along \(D_1\) in the reference configuration.

Fig. 6
figure 6

Kinematics of parallel fibers. Reference configuration (left) and current configuration (right) where the interdistance between the fibers has to be computed

Full size image

By referring to Fig. 6, we calculate the change in distance between the fibers of the same family. In the deformed configuration, the distance between two adjacent fibers, once chosen the point \(\chi \left( \bar{x}_1,\bar{x}_2\right) \), can be defined as the minimum for \(x_1\) of the squared norm

$$\begin{aligned} ||\chi \left( x_1,\bar{x}_2+\varepsilon \right) -\chi \left( \bar{x}_1,\bar{x}_2\right) ||^2 \end{aligned}$$
(6)

The problem consists in solving the following equation

$$\begin{aligned} \frac{\partial }{\partial x_1}\left[ \left( \chi \left( x_1,\bar{x}_2+\varepsilon \right) -\chi \left( \bar{x}_1,\bar{x}_2\right) \right) \cdot \left( \chi \left( x_1,\bar{x}_2+\varepsilon \right) -\chi \left( \bar{x}_1,\bar{x}_2\right) \right) \right] =0 \end{aligned}$$
(7)

Having called \(\varvec{v}^1=\chi \left( x_1,\bar{x}_2+\varepsilon \right) -\chi \left( \bar{x}_1,\bar{x}_2\right) \), we can rewrite Eq. 7 as

$$\begin{aligned} \frac{\partial }{\partial x_1}\chi \left( x_1,\bar{x}_2+\varepsilon \right) \cdot \varvec{v}^1=0 \end{aligned}$$
(8)

or, using the definition of \(\varvec{F}\), as

$$\begin{aligned} \varvec{F}\left( x_1,\bar{x}_2+\varepsilon \right) D_1\cdot \varvec{v}^1=0 \end{aligned}$$
(9)

For solving this equation, we must expand the function \(\chi \left( x_1,\bar{x}_2+\varepsilon \right) \)

$$\begin{aligned} \chi \left( x_1,\bar{x}_2+\varepsilon \right)\simeq & {} \chi \left( \bar{x}_1,\bar{x}_2\right) +(x_1-\bar{x}_1)\left. \frac{\partial \chi }{\partial x_1}\right| _{\left( \bar{x}_1,\bar{x}_2\right) }+\varepsilon \left. \frac{\partial \chi }{\partial x_2}\right| _{\left( \bar{x}_1,\bar{x}_2\right) } \end{aligned}$$
(10)

so

$$\begin{aligned} \varvec{v}^1\simeq \bar{\varvec{F}}\left( \begin{array}{c} x_1-\bar{x}_1\\ \varepsilon \end{array}\right) \end{aligned}$$
(11)

where \(\bar{\varvec{F}}=\varvec{F}(\bar{x}_1,\bar{x}_2)\). As done for \(\chi \left( x_1,\bar{x}_2+\varepsilon \right) \), we perform the same kind of expansion for \(\varvec{F}\left( x_1,\bar{x}_2+\varepsilon \right) \) obtaining

$$\begin{aligned} \varvec{F}\left( x_1,\bar{x}_2+\varepsilon \right) =\bar{\varvec{F}}+\nabla \bar{\varvec{F}}\left( \begin{array}{c} x_1-\bar{x}_1\\ \varepsilon \end{array}\right) \end{aligned}$$
(12)

Finally, we can rewrite Eq. 7 as

$$\begin{aligned} \left( \bar{\varvec{F}}D_1+(x_1-\bar{x}_1)\partial _1\bar{\varvec{F}}D_1+\varepsilon \partial _2\bar{\varvec{F}}D_2\right) \cdot \left( (x_1-\bar{x}_1)\bar{\varvec{F}}D_1+\varepsilon \bar{\varvec{F}}D_2\right) =0 \end{aligned}$$
(13)

Now, by considering that \(x_1-\bar{x}_1\) is of the same order as \(\varepsilon \) and neglecting terms of the order \(\varepsilon ^2\) we obtain

$$\begin{aligned} x_1-\bar{x}_1=-\varepsilon \frac{\bar{\varvec{F}}D_1\cdot \bar{\varvec{F}}D_2}{||\bar{\varvec{F}}D_1||^2} \end{aligned}$$
(14)

By substituting this result in Eq. 11, we obtain an estimation for \(\varvec{v}^1\) (and recalling the family index \(\alpha \))

$$\begin{aligned} \frac{\varvec{v}^1}{\varepsilon }\simeq= & {} \bar{\varvec{F}^1}\left( \begin{array}{c} -\frac{\bar{\varvec{F}^1}D_1\cdot \bar{\varvec{F}^1}D_2}{||\bar{\varvec{F}^1}D_1||^2}\\ 1 \end{array}\right) \end{aligned}$$
(15)

A similar discussion can be carried out for the other family of fibers. Similarly to what we have done before, by defining \(\varvec{v}^2=\chi \left( \bar{x}_1+\varepsilon ,x_2\right) -\chi \left( \bar{x}_1,\bar{x}_2\right) \) the final result is

$$\begin{aligned} \frac{\varvec{v}^2}{\varepsilon }\simeq= & {} \bar{\varvec{F}^2}\left( \begin{array}{c} 1\\ -\frac{\bar{\varvec{F}^2}D_1\cdot \bar{\varvec{F}^2}D_2}{||\bar{\varvec{F}^2}D_2||^2} \end{array}\right) \end{aligned}$$
(16)

From a phenomenological point of view, we expect the contact interaction energy term to behave as in Fig. 7, where \(\delta \) is the interdistance between the centerlines of two adjacent fibers (see Fig. 6) and \(\varepsilon \delta _0\) is the contact interdistance between two fibers. This interdistance is a parameter that has to be settled depending on the geometrical features of the considered pantographic structure and it can be simply computed as the thickness of the fiber.

Fig. 7
figure 7

Theoretical behavior of the contact interaction energy versus the interdistance between two fibers

Full size image

A mathematical expression of the behavior in Fig. 7 can be obtained by considering a piecewise continuous function, such as

$$\begin{aligned} u^{\alpha }\left( \chi (\cdot )\right) =\left\{ \begin{array}{cr} \frac{\varepsilon \delta _0-||\varvec{v}^{\alpha }||}{||\varvec{v}^{\alpha }||^2} &{} \text {if}\,\, ||\varvec{v}^{\alpha }||\le \varepsilon \delta _0\\ 0 &{} \text {if}\,\, ||\varvec{v}^{\alpha }||>\varepsilon \delta _0 \end{array}\right. \end{aligned}$$
(17)

The energy terms accounting for contact of fibers can be then written as

$$\begin{aligned} \mathfrak {U}_{\phi _{\alpha }}(\varvec{\chi } ( \cdot ))= \int _{\varOmega }\frac{\mathrm {K}_{\phi _{\alpha }}}{2} \;f\left( u^{\alpha } \right) \; \mathrm{d}\varOmega \quad \text {for}\,\, \alpha =1,2 \end{aligned}$$
(18)

where \(f\left( u^{\alpha } \right) \) is a generic monotonically increasing function of \(u^{\alpha }\).

4.3 Complete form of the strain energy

Actually, we can formally write the complete strain energy of the pantographic structure. It is composed by Eq. 1 and by the new terms we wrote to take into account the effects of interaction between the fibers: (i) of different families (sliding mechanism) and (ii) of the same family (contact). So, we have

$$\begin{aligned} \mathcal {U} (\varvec{\chi } ( \cdot ))= & {} \int _{\varOmega }\sum _{\alpha }\frac{\mathrm {K}_{e}^{\alpha }}{2} (||\varvec{F}^{\alpha } D_{\alpha }|| -1)^{2}\; \mathrm{d} \varOmega \\&+ \int _{\varOmega } \sum _{\alpha }\frac{\mathrm {K}_{b}^{\alpha }}{2}\genfrac[{}{}{}{ \nabla \varvec{F}^{\alpha } \vert D_{\alpha } \otimes D_{\alpha } \cdot \nabla \varvec{F}^{\alpha }\vert D_{\alpha } \otimes D_{\alpha }}{||\varvec{F}^{\alpha } D_{\alpha }||^{2}} \\&- \left. \left( \frac{\varvec{F}^{\alpha } D_{\alpha }}{||\varvec{F}^{\alpha } D_{\alpha }||} \cdot \frac{ \nabla \varvec{F}^{\alpha }\vert D_{\alpha } \otimes D_{\alpha }}{||\varvec{F}^{\alpha } D_{\alpha }||}\right) ^{2}\right] \; \mathrm{d} \varOmega \\&+ \int _{\varOmega }\frac{\mathrm {K}_{s}}{2} \left| \arccos \left( \frac{\varvec{F}^{1} D_{1}}{||\varvec{F}^{1} D_{1}||} \cdot \frac{\varvec{F}^{2} D_{2}}{||\varvec{F}^{2} D_{2}||}\right) -\frac{\pi }{2}\right| ^{\beta } \; \mathrm{d}\varOmega \\&+ \int _{\varOmega }\frac{\mathrm {K}_{fc}}{2} ||\varvec{\chi }^{1} -\varvec{\chi }^{2}||^{2} \; \mathrm{d}\varOmega \\&+ \int _{\varOmega }\sum _{\alpha }\frac{\mathrm {K}_{\phi _{\alpha }}}{2} \;f\left( u^{\alpha } \right) \; \mathrm{d}\varOmega \end{aligned}$$

Finally, we want to deal with dissipation. In fact, as it is pictorially explained in Fig. 5a, one of the possible relative motions which can arise between two parallel fibers consists in a parallel shifting. This shifting motion produces a friction, and it can be modeled by using a Rayleigh dissipation potential [78]. The discussion on how we can include dissipation in the modeling is presented in the next Section.

5 Introducing dissipation: the Hamilton–Rayleigh Principle

As it should be now clear, the approach to Mechanics employed in this paper is variational. The origins of this formulation can be found in Lagrange [79, 80] and Piola [12, 13], and it has been extended to continuum mechanics by Paul Germain and others (we have previously discussed this point). A big question, which is yet debated, consists in the possibility of including dissipation in the variational framework. A possibility is to employ the so-called Hamilton-Rayleigh Principle [48]. It can be summarized as follows:

  1. i.

    we separate the problem in two parts, one conservative and the other dissipative, and we describe only the conservative part using an action functional;

  2. ii.

    we calculate the first variation of the action functional;

  3. iii.

    we add suitable linear functionals to model dissipative phenomena (Rayleigh dissipation functionals);

  4. iv.

    the Hamilton-Rayleigh Principle can be formulated by summing the Lagrangian and the dissipative virtual works and postulating this sum to be vanishing for every admissible variation of motion. This is an energy balance; the share energy not conserved is dissipated indeed.

Before discussing the mathematical formulation of the above explained procedure, we want to remark that it is debated about the real nature of dissipative systems (see, for example, [81]). In fact, it has been conjectured that a dissipative system can always be considered a system that approximates or simplifies a conservative one. This conjecture is supported by the fact that one can obtain a non-Lagrangian system from a Lagrangian one in at least two ways: (i) neglecting some terms in the Lagrangian action functional and therefore obtaining an approximation that is not Lagrangian; (ii) considering that a non-Lagrangian system is in effect a subsystem of a Lagrangian system with many more degrees of freedom. In the latter case, the dissipative phenomena would be only apparent: in fact, the energy which is “lost” by dissipation would simply transit to the hidden degrees of freedom and could even return to the visible ones, but only after the Poincaré time, which can be very long.

In the following, we give the main points of the variational approach to dissipation and we eventually apply these concepts to the case studied in the present work for including dissipation in the study of pantographic structures.

5.1 Least action principle and Euler–Lagrange equations

Consider a Lagrangian action functional

$$\begin{aligned} \mathcal {A}\left( \varphi _{\alpha }(\cdot )\right) =\int _T\mathcal {L}\left( X^{A},\varphi _{\alpha },\frac{\partial \varphi _{\alpha }}{\partial X^{A}}\right) \,\mathrm{d}V \end{aligned}$$
(19)

where: \(\varphi _{\alpha }(X^{A})\) is any set of n tensor fields on \(\mathbb {R}^m\) (\(\alpha =1,\dots ,n\) and \(A=1,\dots , m\)); \(T\subset \mathbb {R}^m\); dV is the Lebesgue measure and the coordinates \(X^{A}\) are used to refer to the generic point X in T. The least action principle states that a real motion for the considered mechanical system minimizes the Lagrangian action functional 19. Moreover, such a motion has to be searched for among the motions for which the first variation of the action vanishes. For searching for an extremal of the action potential, it is necessary to compute its first variation. This computation brings to

$$\begin{aligned} \delta \mathcal {A}\left( \delta \varphi _{\alpha }(\cdot )\right)= & {} \int _{t_0}^{t_f}\int _V\left( \left( \frac{\partial \mathcal {L}}{\partial \varphi _{\alpha }}-\frac{\partial }{\partial t}\left( \frac{\partial \mathcal {L}}{\partial (\partial \varphi _{\alpha }/\partial t)}\right) \right) \delta \varphi _{\alpha }\right. \nonumber \\&+\left. \sum _{A=1}^{3}\frac{\partial \mathcal {L}}{\partial (\partial \varphi _{\alpha }/\partial X^{A})}\frac{\partial \delta \varphi _{\alpha }}{\partial X^A}\right) \,dV \end{aligned}$$
(20)

where the following quantity is generally accounted for as virtual work

$$\begin{aligned} \delta \mathcal {W}(\delta \varphi _{\alpha };t)= & {} \int _V\left( \left( \frac{\partial \mathcal {L}}{\partial \varphi _{\alpha }}-\frac{\partial }{\partial t}\left( \frac{\partial \mathcal {L}}{\partial (\partial \varphi _{\alpha }/\partial t)}\right) \right) \delta \varphi _{\alpha }\right. \nonumber \\&+\left. \sum _{A=1}^{3}\frac{\partial \mathcal {L}}{\partial (\partial \varphi _{\alpha }/\partial X^{A})}\frac{\partial \delta \varphi _{\alpha }}{\partial X^A}\right) \,\mathrm{d}V \end{aligned}$$
(21)

so that the stationarity condition of the action functional can be written as

$$\begin{aligned} \delta \mathcal {A}\left( \delta \varphi _{\alpha }(\cdot )\right) =\int _{t_0}^{t_f}\delta \mathcal {W}(\delta \varphi _{\alpha };t)\,\mathrm{d}t=0 \end{aligned}$$
(22)

The main result of this principle consists in a system of differential equations (and some boundary conditions) called Euler–Lagrange equations (here in a multi-dimensional space of dimension m, which can be 4 for the space-time)

$$\begin{aligned}&\frac{\partial \mathcal {L}}{\partial \varphi _{\alpha }}-\sum _{A=1}^{m}\frac{\partial }{\partial X^{A}}\left( \frac{\partial \mathcal {L}}{\partial (\partial \varphi _{\alpha }/\partial X^{A})}\right) -\frac{\partial }{\partial t}\left( \frac{\partial \mathcal {L}}{\partial (\partial \varphi _{\alpha }/\partial X^{A})}\right) =0\, ,\quad \forall X\in T \end{aligned}$$
(23)
$$\begin{aligned}&\sum _{A=1}^{m}\frac{\partial \mathcal {L}}{\partial (\partial \varphi _{\alpha }/\partial X^{A})}N_{A}=0\, ,\quad \forall X\in \partial T/\partial _d T \end{aligned}$$
(24)

where \(\partial T/\partial _d T\) is the difference between \(\partial T\) and \(\partial _d T\) (\(\partial _d T\subset \partial T\) is the part of the boundary of T where the functions constituting the action functional have an assigned value) and \(N_{A}\) is the extern unit normal to \(\partial T/\partial _d T\).

We remark that the least action principle, when formulated for action functionals admitting first differentials, can be considered as a particular form of the principle of virtual work (see [48, 78] for all the details). Moreover, it can be shown that the principle of least action implies the principle of virtual work.

5.2 The principle of virtual work

As it has been recalled in the previous subsection, the principle of virtual work states that the motion of a system is characterized by imposing that

$$\begin{aligned} \left( \forall t\in \left[ t_0,t_f\right] \right) \left( \forall \varphi _{\alpha }\left( X^1,X^2,X^3\right) \right) \left( \delta \mathcal {W}(\varphi _{\alpha };t)=0\right) \end{aligned}$$
(25)

Following Lagrange and Piola and being inspired by some particular cases, we assume in general that the virtual work functional can be separated in three distinct functionals: inertial work (\(\mathcal {W}^{iner}(\varphi _{\alpha };t)\)); internal work (\(\mathcal {W}^{int}(\varphi _{\alpha };t)\)); external work (\(\mathcal {W}^{ext}(\varphi _{\alpha };t)\)). In the case of Lagrangian systems, i.e., when the expression of virtual work derives from the action functional expressed in terms of a Lagrangian density \(\mathcal {L}\), identifications in terms of the Lagrangian density for the virtual work functionals can be found (see [48] for more details). Using the above introduced inertial, external and internal virtual work functionals, the principle of virtual work stated in Eq. 25 can be rewritten as

$$\begin{aligned}&\left( \forall t\in \left[ t_0,t_f\right] \right) \left( \forall \varphi _{\alpha }\left( X^1,X^2,X^3\right) \right) \left( \delta \mathcal {W}^{TOT}(\varphi _{\alpha };t)\right. \nonumber \\&\quad \left. =\delta \mathcal {W}^{iner}(\varphi _{\alpha };t)-\delta \mathcal {W}^{int}(\varphi _{\alpha };t)+\delta \mathcal {W}^{ext}(\varphi _{\alpha };t)=0\right) \end{aligned}$$
(26)

5.3 Hamilton–Rayleigh principle

The main point of this approach consists in reformulating the Principle of Virtual Work in terms of an Action functional and a Rayleigh dissipation functional. Rayleigh dissipation functional \(\mathcal {R}\) is defined as a linear functional on the set of velocities, not on the set of motions as the action \(\mathcal {A}\). The Rayleigh dissipation functional is a map that associates, at every instant t, a positive real number with the ordered set of fields \(\left( \varphi _{\alpha }, \frac{\partial \varphi _{\alpha }}{\partial X^A}, \frac{\partial \varphi _{\alpha }}{\partial t}, \frac{\partial }{\partial X^A}\left( \frac{\partial \varphi _{\alpha }}{\partial t}\right) ; t\right) \)

$$\begin{aligned} \mathcal {R}\left( \varphi _{\alpha }(\cdot ;t), \frac{\partial \varphi _{\alpha }(\cdot ;t)}{\partial X^A}, \frac{\partial \varphi _{\alpha }(\cdot ;t)}{\partial t},\frac{\partial }{\partial X^A}\left( \frac{\partial \varphi _{\alpha }(\cdot ;t)}{\partial t}\right) ; t\right) \in \mathbb {R}^+ \end{aligned}$$
(27)

where \(X^A\) are Lagrangian coordinates. A general formulation of Hamilton–Rayleigh principle can then be written reformulating the principle of virtual work of Eq. 25

$$\begin{aligned} \left( \forall t\in \left[ t_0,t_f\right] \right) \left( \forall \varphi _{\alpha }\left( X^1,X^2,X^3\right) \right) \left( \delta \mathcal {W}(\varphi _{\alpha };t)=\delta \mathcal {R}(\varphi _{\alpha };t)\right) \end{aligned}$$
(28)

So the principle of virtual work reformulated following Hamilton–Rayleigh in the case of pantographic metamaterial has the form

$$\begin{aligned} \delta \mathcal {W}\left( \delta \chi ^{\alpha }\right) =\frac{\partial \mathcal {R}}{\partial \nabla \dot{\chi }^{\alpha }}\delta \chi ^{\alpha } \end{aligned}$$
(29)

As the Rayleigh dissipation functional is defined on the set of velocities, for computing it we have to consider the velocities of the fibers in the family direction. From a mathematical point of view, we will have

$$\begin{aligned} \mathcal {R}=\sum _{\alpha =1}^2\frac{\mu ^{\alpha }}{2} ||\nabla \dot{\varvec{\chi }}^{\alpha }|\varvec{F}^{\alpha }D_{\alpha }||^2 \end{aligned}$$
(30)

where \(\mu ^{\alpha }\) is a suitably chosen coefficient which has to be calibrated on the physical system.

6 Conclusion

We have proposed an extension to the strain energy for the pantographic metamaterial previously introduced in [3]. The constitutive micro-structural elements have been studied in detail for fully understanding their contribute to the global macroscopic deformation behavior of the structure. Specifically, we have examined the mechanics of pivots (as it has been already presented in [76, 77]) and the interactions between the constituent fibers. In fact, in the available literature about pantographic structures the fibers of the two families have been always considered as non-interacting. We have shown, instead, that the fibers interact between them at least in three different ways: (i) indirectly, through the microstructural connections (the pivots) which could allow for relative sliding between the two fiber families; (ii) directly, as the fiber of one family can touch each other and (iii) can scroll introducing dissipation.

From a mathematical point of view, these effects have been modeled firstly by introducing two placement fields for the two fiber families and adding a coupling term to strain energy and, secondly, by adding two further terms taking into account the interdistance between parallel fibers and a Rayleigh dissipation potential (for modeling the friction).

The present work proposes, as said, a modification to the continuum model of pantographic structures. A further study must be carried out for providing verifications. More precisely, a first step will consist in implementing the new energy terms in the numerical codes available in the study of such structures. Preliminary results for the sliding mechanisms have been obtained in [77] using the continuum model and in [82, 83] employing a “mesoscopic” description in which the pantographic structure is modeled as an assembly of Euler–Bernoulli nonlinear beams. For the complexity of the topic, it will be necessary to rely on powerful numerical tools. Many examples can be found in the literature [84,85,86,87,88,89,90,91,92,93,94,95,96], and some numerical strategies have been already successfully employed in the case of pantographic structures [5, 97,98,99,100,101,102,103,104].

A second stage will concern the experimental validation. Many experiments have been performed for testing the pantographic metamaterial [105,106,107,108,109,110,111]. Clearly, the model validation will be the more reliable the more accurate the experimental data and their analysis will be. A recent data analysis procedure that is very promising for its accuracy and benefits is Digital Image Correlation. Through this technique (whose basic skills can be found in [112,113,114]), it is possible to obtain an estimate of the displacement and deformation fields. Moreover, this technique can be applied in the calibration of the constitutive parameters of the model to be tested. Some results already available in the field of pantographic structures can be found in [101, 115]. From the interaction between experiments and theory, mediated by numerical simulations, it is possible to design a pantographic structure with optimal mechanical properties. This has been performed in a preliminary way in [116].