Equilibrium of Two-Dimensional Cycloidal Pantographic Metamaterials in Three-Dimensional Deformations

A particular pantographic sheet, modeled as a two-dimensional elastic continuum consisting of an orthogonal lattice of continuously distributed fibers with a cycloidal texture, is introduced and investigated. These fibers conceived as embedded beams on the surface are allowed to be deformed in a three-dimensional space and are endowed with resistance to stretching, shearing, bending, and twisting. A finite element analysis directly derived from a variational formulation was performed for some explanatory tests to illustrate the behavior of the newly introduced material. Specifically, we considered tests on: (1) bias extension; (2) compressive; (3) shear; and (4) torsion. The numerical results are discussed to some extent. Finally, attention is drawn to a comparison with other kinds of orthogonal lattices, namely straight, parabolic, and oscillatory, to show the differences in the behavior of the samples due to the diverse arrangements of the fibers.


Introduction
In recent decades, metamaterials has attracted the attention of the scientific community considerably.They are artificially constructed materials characterized by an interior architecture that is designed with the aim of producing an optimized combination of responses to specific external excitation.One may conceive many possible metamaterials depending on the particular properties to be obtained.For example, metamaterials may be mechanical, acoustic, electromagnetic, optic, and so on.Typically, metamaterials show qualitative exotic behaviors due to their substructure rather than to the composition of the material of which they are made [1][2][3][4].
In this paper, we investigate a particular type of mechanical metamaterial [5][6][7]: a pantographic sheet [8][9][10][11][12][13][14][15].This particular material has been introduced as the first synthetic example of material exhibiting a strain-energy function involving second gradient contributions relevant at a macroscopic scale of observation [16][17][18].This main feature is inherited by the presence of "fibers" that, having a resistance to being bent, require the presence of spatial second-order derivatives of the displacement to properly describe their energy.
It is worth noticing that the pantographic structure may be used as a reinforcement of composite material, e.g., a two-phase continuum with a "pantographic-inspired" architecture, which is a 3D fiber network, surrounded by a matrix, capable of large deformations [19][20][21][22].
From a geometric point of view, a pantographic surface resembles very closely an elastic network structure.Therefore, many of the achievements in modeling these systems can be fruitfully employed also for the pantographic sheets (see, e.g., [23][24][25][26][27][28][29]).According to these results, the pantographic structure is identified with a surface composed of two families of continuously distributed fibers in the framework of continuum theory.Indeed, this kind of perspective is also justified by homogenization procedures applied to the discrete system of fibers [30][31][32][33][34][35][36].Besides, since the fibers are modeled in accordance with the beam theory, which entails the description of their cross-section rotation, the pantographic sheets may be regarded as micropolar continua.Thus, results from this research topic are available too (see, e.g., [37][38][39]).
For the sake of completeness, we have to mention that pantographic structures, as microstructured materials, may also show impressive dynamical features related to the complexity of the interior architecture [40][41][42][43][44][45].Besides, in this context of engineered materials, it is also worth mentioning the study of dynamically self-organized "intelligent" materials [46] or functionally graded materials [47].
The paper is organized as follows.Section 2 describes the pantographic system with initially cycloidal fibers.Section 3 is devoted to giving some numerical examples of the behavior of this arrangement of fibers.Finally, Section 4 reports a comparison with other kinds of fiber dispositions and discusses the differences.

Pantographic Sheets with Initially Cycloidal Fibers
The pantographic sheet investigated is characterized by an orthogonal lattice of fibers initially curved and arranged parallel to two families of cycloids.The sample under study is rectangular and cut suitably in order to have a fiber disposition symmetric with respect to the longitudinal axis of the specimen and with a higher density along the most extended lateral edges.This arrangement of the fibers has been specially conceived with the purpose of employing the superior capability of curved fibers to bear loads (fibers tend to behave like arches) and to have reinforcement in long boundaries (see Figure 1).Therefore, the sample has a notably anisotropic behavior with an extremely soft core (due to the pantographic effect [16,17,[48][49][50][51]) and a stiffer part at the borders.
The two families of curves with one parameter, namely ϕ and ψ, in an initial rectangular domain B of size 2R × 2πR are represented in Cartesian coordinates, for the first family, as and, for the second family, as where α and β are constants that identify the particular curve belonging to the family.The deformation of the pantographic surface is defined by the map r(X 1 , X 2 ) : B ⊂ R 2 → R 3 , which takes the material particle initially placed in X = (X 1 , X 2 ) ∈ B and yields the current position of that particle in a three-dimensional Euclidean space.
By assuming a Cartesian coordinate system {e i } with e 1 and e 2 parallel to the edges of the sample, and introducing the components of displacement u i in that basis, the explicit expression of the placement map becomes where Greek indexes range from 1 to 2 while Latin ones from 1 to 3. The differentiation with respect to the Cartesian coordinates of r is represented by the following notation Recalling the definition of the deformation gradient F = ∇r, we have The Cauchy-Green deformation tensor is thus given by and the strain tensor E can be evaluated, by components, as follows where δ αβ is the Kronecker delta.The considered problem also needs to explicate the second gradient of the deformation, ∇F = ∇∇r, namely ∇F = F ,α ⊗ e α , because it is able to describe the curvatures and the twist of the fibers [52,53].Indeed, we recall here that the considered continuous surface can be thought of as an infinite distribution of curved fibers arranged in the domain, as mentioned above (see also [54]), assumed to be material curves with no relative slipping and tied together at their points of intersection.
According to Steigmann and dell'Isola [53], a suitable fiber decomposition is introduced by means of a couple of orthogonal unit vector fields {L(X), M(X)}, which determines the fiber directions in the plane of B for the reference configuration.Besides, the fibers directions in the current configuration, i.e., after a deformation, are given by the unit vectors {l(X), m(X)} as follows where λ and µ are the stretches of fibers.Knowing the current directions of fibers {l, m}, which define the tangent plane to the deformed pantographic sheet at the material point initially in X, it is possible to evaluate the shear distortion angle γ, i.e., the change of angle passing from {L(X), M(X)} to {l, m}, by The tangent vector t ϕ to each curve of the family in Equation ( 1) is obtained by Using Equation (1), with simple algebraic calculations, we obtain and, exploiting the properties of the trigonometric functions, Choosing a root of Equation ( 12), the parameter ϕ can be expressed in terms of Cartesian coordinates as follows therefore, the tangent vector to any curve of the family in Equation ( 1) at the point X, in Cartesian coordinates, becomes and its norm is Eventually, the unit vectors L(X) and M(X) are In Equations ( 16), we take advantage of the fact that the two families of curves in Equations ( 1) and ( 2) are orthogonal in the intersection point.Alternatively, M(X) may be evaluated with the same strategy used for L(X) applied to the curves in Equation (2).
Using the fiber decomposition previously defined, the gradient of deformation is thereby the Cauchy-Green deformation tensor takes the expression: and the second gradient of the deformation becomes [53] ∇∇r where and η l and η m being the geodesic curvatures of the deformed fibers, namely the part of the fiber curvatures in the tangent plane of the deformed sheet, and φ l and φ m the so-called Tchebychev curvatures.Besides, define the orthogonal directions of the fibers embedded on the pantographic sheet in the current configuration, whereas in which κ l and κ m are the normal curvatures of the deformed fibers, and τ represents the twist of these fibers.These last three measures of deformations vanish if the deformed equilibrium shape of the sample remains plain.As a result, they are related to the fiber curvatures and the torsion activated by an out-of-plane flexure of the pantographic surface.
To specify the behavior of the considered bi-dimensional continuum, according to the authors of [53,55], we assume a strain-energy density incorporating the curvilinear orthotropic symmetry due to the initial fiber geometry such that where the first term depends on the variables (see, e.g., [53,55]): which represent, respectively, the fiber extension along L and M directions, as well as the area dilation.
In particular, the first term of the energy density, W I , is assumed to be where Y L , Y M , and G LM are positive material constants.The second gradient energy term, W I I , is postulated as follows and k T are material parameters.In particular, to evaluate the terms of this second part of the energy, we employ the representation [53]: and the projections on the tangent plane to the deformed surface as well as on the normal direction n.

Numerical Examples of Equilibrium Shapes
Under the framework of the weak formulation, Finite Element simulations were performed using the model previously described and proposed in [53].The energy density in Equation ( 24) was directly used in the software COMSOL Multiphysics R because it very quickly allowed the definition of a custom form of strain-energy function.Afterwards, the program evaluated the related weak form of the equilibrium equations together with its finite-element implementation.To deal with the particular problem at hand, which entailed the presence of second derivatives of the unknown kinematical variables sought for, we used Argyris finite elements [63,64], which are C 1 /H 2 -conforming.Indeed, the related shape functions were particular Hermite quintic polynomials defined over a triangular element with 21 degrees of freedom.This element provides continuity for the values of the shape functions as well as for all their first and second partial derivatives at the vertices and for the integrals of the normal derivatives along edges.
The maximum size of mesh elements was set to ≈ 10 −1 R after a convergence study [65], which involved a total of ≈ 4.6 × 10 4 degrees of freedom for the whole sample because each adopted element has 21 degrees of freedom.
For the sake of simplicity, the equilibrium mechanical problem to be solved was rewritten in a non-dimensional form by normalizing the strain-energy density in Equation ( 24) with a reference stiffness Y 0 : and using a reference radius, R (the one used to define the cycloids), to normalize lengths.Herein, a superimposed tilde denotes non-dimensional quantities.Therefore, the material parameters appearing in Equations ( 26) and ( 27) were normalized in the following way: Specifically, the values of these parameters employed in the numerical analysis presented below were: ỸL = 100, ỸM = 100, GLM = 0.2, The considered examples were selected to illustrate the main behavior of the cycloidal arrangement of the fibers.Explicitly, we considered a bias extension test, a compressive test, a shear test, and a torsion test.All tests were conceived to show the occurrence of large deformations.We emphasize that the deformed equilibrium configurations displayed in all figures of this paper are given at the same scale as the geometry of the pantographic sheet, i.e., there is no magnification of the deformation.Moreover, for the sake of understandability, in all the equilibrium shapes, we also plot with solid black lines the material curves representing the fibers in the current configuration.
The paper aims to provide an initial characterization of the considered pantographic sheet.To fulfill this purpose, the commercial code used is sufficiently adequate.However, we are well-aware that, to have an efficient numerical implementation for complex systems such as those analyzed, ad hoc codes must be implemented optimized for each particular case investigated.As examples of numerical formulations, reference may be made to the works in [66][67][68][69] related to a Hencky discrete approach, that in [70] for Ritz's method, those in [71][72][73][74][75][76] for shell or plate structures, [77][78][79][80] for rods, and those in [81][82][83][84][85] regarding applications to composites and lattices.
Herein, we considered only the study of static cases, neglecting inertial contributions.However, we plan to systematically investigate bifurcation scenarios and, primarily, the dynamics of deformation.Naturally, this requires time-dependent weak-formulation techniques based on the introduction of kinetic energy which may depend on the particular microstructure of the pantographic sheet.

Bias Extension Test
In the bias extension case, the cycloidal pantographic sheet was clamped to the left short end, and a displacement parallel to the longitudinal direction of the sample was prescribed on the opposite edge.In this first case, only Dirichlet conditions on displacement were considered.Therefore, we applied no double forces or couples.Figure 2 reports the equilibrium shape for an imposed displacement ũ1 of 1.5 and ũ2 being zero on the same edge.The lack of symmetry of the cycloidal lattice with respect to a vertical central axis reflects in the equilibrium shape, which exhibits the same feature.In addition to the distribution of the energy densities W I and W I I , the figure shows also some material lines in the current configuration which initially lay along cycloids as in Figure 1.In Figure 2, it is possible to see that both the energy densities are concentrated in a narrow layer at the long edges, as is expected, given the greater density of the fibers in that region and the more compliant disposition of those in the center.
Figure 3 shows the contributions of the elastic energy versus the imposed displacement.We split the energy due to the first gradient in a piece accounting for the elongation of the fibers, namely the integral over the domain B of the first two terms of Equation (26), Ψ e , and in a piece due to the shear deformation, i.e., the integral of the last term of Equation ( 26), Ψ s .Moreover, the part related to the second gradient responsible for the geodesic bending in this in-plane motion is also reported, i.e., Ψ b , as well as the total energy stored in the deformation Ψ tot .The most significant contribution is related to fiber elongation.Similarly, the accumulation of the deformation along the most extended edges is visible also in Figure 4 that reports the distribution of the first gradient measures of deformation, namely ε L , ε M , and γ.The shear distortion γ is an exception because, for a certain amount, it is also spread in the central region.

Bias Compression Test
In a compression test, the analyzed surface with embedded fibers shows a buckling behavior.In detail, we considered a zero displacement on the short left side as in the previous case, while the imposed displacement in the opposite edge decreased until ũ1 reached the value of −2.Besides, here, we also added the boundary conditions related to the normal derivative of the out-of-plane displacement, which was set to be zero on both the short edges, i.e., ũ3,1 =0.
Figure 5 shows the buckled shape (Figure 5a) at the final imposed displacement.Figure 5b exhibits the second gradient contribution of the energy density since it is the most appreciable in such a case.It is also emphasized that here a lack of symmetry is perceivable if we observe the energy distribution and the material lines in the current configuration.

Shear Test
The equilibrium shape obtained in a shear test is presented in Figure 6.In this case, we considered the same boundary conditions of the bias test but the displacement imposed on the right short side was now ũ2 equal to 2 while ũ1 was fixed to zero on the same edge.In addition, in this case, both the energy density of the first and second gradient tend to be concentrated close to the long edges.The energy contributions depicted in Figure 7, differently from the bias extension test, show that the stretching and the geodesic bending energies are almost the same until the displacement equals 1.2.Similar considerations can be done in this case as for the bias extension test, for the first gradient measures of deformation, as can be detected by Figure 8.

Torsion Test
The last example is a torsion test where the displacements on the left short side were fixed, while the displacements of the short right side were evaluated to obtain a rigid rotation of an angle θ with respect to the central longitudinal axis of the sample.Moreover, the derivative with respect to the X 1 coordinate of the out-of-plane displacement was also put to zero on the short sides.In Figure 9 are reported two equilibrium shapes for a torsion angle of 1.3π and 1.5π in Figure 9a,b, respectively.The former shape resembles a helix, while the latter equilibrium configuration is a buckled one.Indeed, in Figure 9b, the deformed surface presents a particular folding near the long borders deviating from the helix-type shape.Figure 9c shows with the colors the distribution of the second gradient energy density for the buckled shape.The presence of the buckling is evident if we see the plot of the longitudinal reaction on the short left side versus the torsion angle θ in Figure 10.After an angle of 1.4π, a drop in the reaction R X 1 due to the buckling behavior is quite evident.The presence of the buckling is also detectable looking at the energy contributions varying the rotation angle θ, as indicated by the diagrams in Figure 11, in correspondence of the suddenly change of the slope of the energy plots.Moreover, we observe that, after the occurrence of the buckling, the flexural/twisting energy increases notably while the extension energy decreases.This fact is conforming to the above-mentioned folding at the sides.

Comparison among the Pantographic Sheets with Orthogonal Lattices
In this section, we compare the results obtained for the cycloidal lattice with those related to linear and curvilinear lattices, especially for parabolic [55] and oscillatory ones [65].The bias extension test described in Section 3.1 was numerically performed for all the cases in the same conditions adopted for the cycloidal sample, namely the size of the domain, the material parameters, and the boundary conditions.
In Figure 12, we plot the particular arrangement of fibers conceived for the parabolic pantographic sheet in the considered domain.Figure 13 shows the main characteristics of this pattern of fibers resulting from the bias extension test.In detail, due to the straight transversal fibers at the center of the specimen, the sample may experience a buckling in an extension test [86].In such a test, the central transversal fibers are subjected to a significant compression indeed.Figure 13b,c displays the buckled shapes for different initial imperfection.If the deformation remains in the plane of the sample, we have a geodesic buckling, as reported in Figure 13b; otherwise, the buckled configuration assumes a typical spoon-shaped form (see Figure 13c and [87] for experimental evidence).Figure 13a shows the equilibrium configuration without any buckling.According to Giorgio et al. [65], the oscillatory network of fibers is described by the two families of curves, as represented below in Cartesian coordinates: where the amplitude A and the frequency of oscillation ω determine the behavior of the pantographic surface; in fact, the value of ω defines the number and the position of the regions in which the fibers tend to be denser, as shown in Figure 14, where a transverse accumulation of the fibers is evident.Of course, α and β select the specific curve of each family.In this case, the distribution of the fibers therefore entails a sort of transverse stiffeners which can be properly designed by choosing the values of A and ω.In what follows, we set A = 0.5 R and ω = 5/(2 R); Scheme of an oscillatory orthogonal network of fibers.Solid blue lines are the graphs of Equation ( 32), while solid red lines refer to Equation (33).
Figure 15 shows the different contributions of the elastic energy versus the imposed displacement in the bias extension test.The oscillatory pantographic sheet is stiffer than the cycloidal one: for the same displacement imposed, it stores twice as much deformation energy.Moreover, the shear energy is almost negligible.
Figure 16 displays by colors, for an imposed displacement ũ1 = 1.65, the distributions of energy density related to the stretching (Figure 16a), the shear (Figure 16b), and the second gradient energy W I I contribution (Figure 16c).Figure 16 also shows the equilibrium configuration in the considered test.As it is expected, in correspondence of the "stiffeners" due to the pattern of the fibers, the transverse deformation of the sample is less pronounced.
Finally, all the longitudinal reactions R X1 related to the different kinds of fiber arrangement are plotted in Figure 17 as a function of the imposed displacement for the sake of comparison.In the case of parabolic fibers, not only is the case without buckling plotted but also the two diagrams of the two detected bucklings are shown.From them, it is clear that the out-of-plane buckling is the most advantageous from an energetic point of view; this is the reason for which it is easily obtainable by experiments.The sample with straight fibers is the most compliant because of the pantographic effect, which is dominant.On the other hand, the oscillatory lattice corresponds to the stiffest one because of the multiple extra rigid regions of the pattern.The cycloidal lattice has an intermediate behavior due to the simultaneous presence of a pantographic effect in the central zone and the particular orientation of the fibers along the more extended edges.Indeed, this particular fiber arrangement contributes to the longitudinal strength of the sample.We further remark that the cycloidal network has a load-bearing capacity almost similar to the parabolic lattice in the case of a spoon-shaped buckling.The case of a parabolic lattice without buckling is a metastable configuration.Thus, it is not attractive in many applications.A little more enticing may be the parabolic net with in-plane buckling in those cases in which there are extra constraints that compel the deformed configuration to remain plain.

Figure 1 .
Figure 1.Scheme of a cycloidal orthogonal network of fibers.Solid blue lines are the graphs of Equation (1), while solid red lines refer to Equation (2).

Figure 2 .
Figure 2. Distribution of the strain-energy density for the bias extension test: (a) W I contribution; and (b) W I I contribution.

Figure 4 .
Figure 4. Distribution of the deformation measures of first gradient for the bias extension test: (a) ε L ; (b) ε M ; and (c) γ.

Figure 5 .Figure 6 .
Figure 5. Bias compression test: (a) buckled shape, colors indicate the out-of-plane displacement; and (b) second gradient energy contribution, W I I .

Figure 9 .Figure 10 .Figure 11 .
Figure 9. Torsion test: (a) equilibrium shape without buckling, where colors indicate the out-of-plane displacement; (b) buckled shape, where colors indicate the out-of-plane displacement; and (c) buckled shape, where colors indicate the second gradient energy W I I .

Figure 12 .Figure 13 .
Figure 12.Scheme of a parabolic orthogonal network of fibers.