A numerical survey of nonlinear dynamical responses of discrete pantographic beams

Materials and structures based on pantographic cells exhibit interesting mechanical peculiarities. They have been studied prevalently in the static case, both in linear and nonlinear regime. When the dynamical behavior is considered, available literature is scarce probably for the intrinsic difficulties in the solution of this kind of problems. The aim of this work is to contribute to filling of this gap by addressing the dynamical response of pantographic beams. Starting from a simple spring mechanical model for pantographic beams, the nonlinear equilibrium problem is formulated directly for such a discrete system also considering inertia forces. Successively, the solution of the system of equilibrium equations is sought by means of a stepwise strategy based on a non-standard integration scheme. Here, only harmonic excitations are considered and, for large displacements, frequency-response curves are thoroughly discussed for some significant cases.


Fig. 1
Geometry and discrete mechanical model of a pantographic beam. Within this study, the following choice has been done for the angle α and the total cells number L/ f : α = π/4 and L/ f = 21 regime. 1 Recently, some interesting papers were published analyzing the dynamics of two-, see [17], and three-dimensional, see [18], beams, also in view to applications to robotic arms [19,20].
In the available literature, there exist several integration schemes for linear and nonlinear mechanical problems. Interested readers can find a large number of references and an extensive discussion thereof in Wriggers' book [21,Chap. 6]. Here, the integration scheme proposed by Casciaro in 1975 [22] will be used. Such a method is practically unrecognized, even though it has some very interesting peculiarities.
Synthetically, Casciaro's scheme is based on: (i) a second-order B-spline interpolation for the displacements in the time domain; (ii) the discrete form of the impulse-momentum relationship at each time step; (iii) the optimization of two free parameters to produce the best, in a sense which will be clarified below, response for a linear elastic one-degree-of-freedom system and for an assigned time step length.
In simple words, Casciaro's scheme guarantees for a chosen time step length, numerically stability, besides being able of: (i) filtering of the roundoff errors and (ii) avoiding beat phenomenon for large time steps. We remark that for nonlinear analyses, e.g., for large displacements, there is no guarantee that the response produced by Casciaro's scheme is the best possible response. However, to the Author's knowledge, this is true for any integration scheme reported in the technical literature.
The paper, after this short Introduction, presents in Sect. 2 the intrinsically discrete model for pantographic beams moving in plane. Successively, in Sect. 3, the employed stepwise solution strategy is briefly discussed. Section 4 reports several numerical simulations, each one being thoroughly discussed. Finally, Sect. 5 presents some concluding remarks and future challenges.

Discrete model for pantographic beams
We consider a planar pantographic beam, in brief p-beam, such as that depicted on the top in Fig. 1. The p-beam is made up of two families of extensible Euler-Bernoulli beams connected at their intersecting points by pivots. As already experimented in [16,[23][24][25] for pantographic lattices with orthogonal beams, in [26] for non-orthogonal beams and in [27,28] for pantographic beams, we describe the mechanical behavior of the p-beam by the system of springs depicted on the bottom in Fig. 1.
The interactions consist in extensional and rotational springs to represent in a very simple form the resistance of beams to extension and bending and the torsional stiffness of the pivots. This mechanical description may be improved in a h-refinement way as reported in [29]. Thence, the strain energy of each deformable spring of the p-beam can be written as: 2 E a = 1 2 a 2 , stretching energy, , shear strain energy, being E a , E b , E c the strain energies, respectively, stored in extensional springs, in flexural-rotational springs and shear-rotational springs. The corresponding springs are colored in red and blue for extensional and flexural springs and in black for shear-rotational springs, see again the bottom part of Fig. 1. The associated strain measures for the stretching, cos β for the bending and γ for the shear strain are all defined in terms of the current positions of the p-beam nodes or, equivalently, in terms of nodal displacements. 3 In particular, we assume as strain measures the following: where P i and p i are, respectively, the reference and current positions of the ith node (an analogous convention holds for the nodes j, k and m). The interested reader can find more details in [15]. Since for the ith node p i = P i + u i , the strain energy element contribution depends only on the involved displacement (namely u i and u j for the stretching, u i , u j and u k for the bending and u j , u m and u k for the shear strain). Stiffnesses of the springs involved in the p-beam are denoted by a, b and c, for extensional, flexural and shear-strain springs, respectively. 4 Strain energies for each spring of the considered p-beam completely define its mechanical behavior. We remark that different forms for the bending strain energy E b are possible, see, e.g., [30,31]. For instance, the classic form E b = 1 2b β 2 could equivalently be used where the stiffness parameterb, different from b, is introduced. It is easy to prove that for β's in a small neighborhood of π (i.e. value assumed by β in the reference zero-energy configuration) the two choices for E b are equivalent. The strain energy E of the whole p-beam can be obtained simply by adding the aforementioned contributions. Collecting the nodal displacements, i.e. the Lagrangian variables, in the vector u, we can define the structural reaction s as the gradient of the strain energy respect to the Lagrangian variables Furthermore, we can define the stiffness matrix K as the Hessian of the strain energy or, equivalently, the gradient of the structural reaction The gradient and the Hessian of the strain energy are the tools necessary to tackle the incremental-iterative strategy capable to solve the dynamical equilibrium equations. 5 2 In order to simplify formulas, we omit to utilize a symbol indicating the corresponding spring in space. 3 Nodes are sketched in Fig. 1 as small circles located at the intersection points of the beams. 4 Here, for the sake of simplicity, we do not assume as different the stiffnesses of the two families of oblique beams and, in addition, we consider only one parameter c to describe the shearing stiffness. 5 We remark that both the structural reaction and the stiffness matrix can be assembled using the standard element-by-element procedure, i.e. computing the element structural reaction s e and the element stiffness matrix K e . A remarkable reduction in the required computer memory can be obtained by exploiting the symmetry and the sparsity properties of the stiffness matrix to record compactly its entries.
Besides the strain energy, the study of the dynamical behavior of pantographic beams requires to define the kinetic energy of the system. We distinguish two contributions: that of beams and that of pivots. The first one, for a beam element of length e , is defined as being ρ e the mass density, A e the area of the beam cross section, v the velocity vector written by using the matrix of the blending functions B e and the Lagrangian parameters collected inu e . The matrix B e , using the local coordinate 0 ≤ x ≤ e , assumes the form: In the case of beams with uniform cross section, some straightforward calculations lead to In addition, the mass matrix contribution related to the pth pivot can be derived from its kinetic energy defined as being ρ p , A p , v p and I the density, the cross-section area, the velocity of pth pivot and the 2×2 identity matrix, respectively. 6 The vectoru p collects the Lagrangian parameters which represent the velocity of the pth pivot and is connected to the pivot velocity vector v p = Iu p . Consequently, pth pivot mass matrix contribution is where d p and h are the diameter and the height of the pth cylindrical pivot, respectively. 7 The mass matrix contribution M e of the eth element connecting nodes i and j is given by Eq. (7). This contribution derives from the ordering chosen for vector of nodal velocitiesu e , i.e. x 1 -and x 2 -velocity for the i-node and x 1 -and x 2 -velocity of the j-node, see Fig. 1. Similarly, the pth pivot mass matrix contribution M e is given by Eq. (9)-always with the same ordering rule for the component velocities, i.e. first the x 1 -and then the x 2 -velocity. The mass matrix elementary contributions (7) and (9) are the only quantities necessary to assemble the global mass matrix M.
Kinetic and strain energy, as defined in the immediate foregoing, allow to compute the mass matrix M and the structural reaction s, which are the only quantities strictly necessary to write the system of differential of equations ruling the motion. In the Euler-Lagrange form, see [33], such a system reads 8 where, besides the Lagrangian parameters used to describe the motion, i.e. the displacement vector u, the velocity vectoru and the acceleration vectorü, also the external force vector f appears. 6 Coherently with the hypotheses on the strain energy, the rotational inertia of the pth pivot is neglected. 7 The diameter of the pivots may be different, in general, for each pivot. Indeed, [32] seeks to optimize their diameters. Contrarily, their height h should be considered as constant if we assume a motion in plane. 8 Even if not described in this work, the stepwise strategy which we discuss below can be easily modified to include dissipative forces, see, e.g., [34][35][36] for an insight.

Solution of nonlinear equilibrium equations by a stepwise strategy
As mentioned in Sect. 1, the technical literature is dense of contributions devoted to solve nonlinear equilibrium equations. Some examples are reported in [37][38][39] and references therein. The book [21], updated up to the first years of 2000, provides with a very good introduction on stepwise integration schemes. Starting from a suitable time discretization, possibly non-uniform, these approaches reconstruct the dynamical equilibrium path-i.e. the pair time t and Lagrangian variables collected in the displacement u and in the velocityu vectorsby solving a sequence of initial value problems. Roughly speaking, for each time interval t corresponding to the chosen time discretization, the solution at the final point is computed on the basis of the solution at the beginning of the time step. The choice of the time step length is crucial-especially for nonlinear problemsand, generally speaking, it requires both some preliminary analyses and a noteworthy experience of the analyst.
Actually, there exists in the technical literature a scheme which does not require, at least in principle, a preliminary study. This numerical integration scheme was developed and published by Casciaro [22], notwithstanding that it was communicated on an important journal such as Meccanica, Casciaro's paper is practically unknown to the international scientific community. 9 The keynote idea of Casciaro's numerical integration scheme is based on the use of two free parameters which are determined by optimizing the dynamical response for the a priori chosen time step length. In practice, Casciaro's idea stems from the observation that the numerical solution of a problem in time is affected by errors depending on the chosen time step length. These errors may be triggered by choosing the time step length t as greater than a fraction of the fundamental (minimum) period of the underlying problem (numerical instability). Furthermore, also the roundoff errors, i.e. when the selected time step length t is lower than a specific threshold, can affect significantly the problem solution (instability dues to roundoff errors). Finally, a relatively large time step length t might introduce spurious beatings, i.e. beating due to numerical errors, even for unconditionally stable schemes.
Essentially, Casciaro's algorithm is based on: (i) a parabolic interpolation law for the displacements in the current time interval, using as parameters the values of the displacement and of the velocity at the beginning of the step and of the velocity at the end of the step; 10 (ii) the discrete form of the momentum-impulse relationship.
The displacement vector can be written at the time t within the time step t 1 − t 0 in the form where u 0 ,u 0 andu 1 are the displacement vector, the velocity vector at the beginning of the time step and the velocity vector at the end of the time step, respectively; β 0 and β 1 are two weighting factors, hence β 0 +β 1 = 1, which may be chosen to optimize the integration scheme results. Interpolation (11) when evaluated at the end of the time step, furnishes the vector u 1 collecting the displacements at the end of the time interval: The discrete-form of the momentum-impulse relationship read as: where we notice the momentum change M(u 1 −u 0 ) and the average net impulse, written using the weighting factors α 0 and α 1 , such that α 0 + α 1 = 1, which, once again, can be tuned to optimize the integration scheme.
Remark that owing to the constraint relationships α 0 + α 1 = 1 and β 0 + β 1 = 1, only two parameters among α 0 , α 1 , β 0 , and β 1 are independent, say α and β. These independent parameters can be computed optimally on the basis of the estimates of the first, T 1 , and of the last, T n , natural periods of the considered structure. Casciaro provided in [22] the analytical expressions of the parameters α and β. In formulae, in the case t < T n /2: where τ = 2π t T n . The limit of Eq. (14), In the case t > T n 2 , the optimal values are: Formulae (15) contain two terms. The first one is a hyperbola branch which matches the optimal values of α and β at the end of the interval t < T n /2. The second one is the ratio between two cubic polynomials which tends-when t goes to infinite-to the value 1 2 ; therefore, such a value of α and β which is able to reproduce a quasi-static solution. Figure 2 shows a plot of the optimal values for α e β obtained using the aforementioned expressions.
In the paper [22], Casciaro does not give many details on the proposed optimization scheme. This may be one of the causes implying the lack of success of his method. Actually, a more detailed description of the integration scheme is reported [41], although both in [22] and in [41] there are some typos in the formulas. The more recent work [42] contains some insights into the method used to build the optimal formulas (14) and (15). We have now almost all the tools for the solution of dynamic nonlinear problems involving multidegree-of-freedom systems: a simple, yet effective, mechanical model such as that described in Sect. 2 and an algorithm to tackle the time step solution in the dynamic case such as that described in this section. We only need a solution strategy which takes into account the geometrical nonlinearities of the p-beam response. Hence, we can express the structural reaction s 1 by using the first-order Taylor expansion having used the stiffness matrix K 0 , i.e. the stiffness matrix computed at the beginning of the considered step, see Eq. (4): Considering again Eqs. (12) and (13), we can eliminate the unknown displacement variables, i.e. u 1 by using Eq. (12). The immediately above equations along with Eqs. (12)- (13) are enough to build a Newton-like algorithm for our nonlinear problem.
In a synthetic way, we can start from the definition of the rest computed at the jth iteration and compute the new estimates of the solution in a Newton-like fashioṅ having introduced the iteration matrix H defined as We have thus completed the introduction of the tools of the numerical strategy used to solve the nonlinear equilibrium problem in the considered time step. At this point, some remarks are deemed useful: (i) Newton's method is not able to get around limit points of the iteration matrix. For this reason, in static problems where the iteration matrix is the stiffness matrix, a correction strategy, such as the arc-length method proposed by Riks [43], is capable to overcome this limitation of the Newton's method. In the dynamic case, the iteration matrix contains also the contribution of the mass matrix and, therefore, Riks correction is not necessarily required. (ii) The algorithm depicted above can be easily applied also in the presence of dissipative terms. In this case, if we consider dissipative forces which can be expressed as d = d(u), it is straightforward to show that they appear in the rest of the equilibrium equations as follows: while the iteration matrix H j includes, apart from the mass and the stiffness matrices, also the dissipation matrix C j defined as and can therefore be written as

Numerical simulations
In this section, we begin the exploration of the dynamical behavior of p-beams under harmonic load by means of a parametric study. In synthesis, we are interested to study: (1) the natural frequencies and the corresponding natural modes; (2) the mechanical response for different values of the excitation frequency; (3) the possible presence of jumps in the frequency-response curve indicating dynamical instabilities; (4) the possible presence of secondary and internal resonances.
To this end, we selected four representative datasets characterized by the stiffnesses reported in Table 1. The choice of these datasets considers both experimental data obtained from testing of 3D-printed pantographic specimens and simplified models currently employed in the literature for analyzing pantographic lattices. The spring stiffnesses of the first dataset derive from an identification process described in [44], which was performed on a 3D-printed specimen obtained from selective laser sintering of polyamide powder having a density ρ = 930 kg/m 3 . The specimen is constituted by beams having square cross section with side equal to 1 mm and by cylindrical pivots with height equal to 2 mm and circular cross section with diameter equal to 0.9 mm. Starting from the dataset based on experimental data-dataset 1-we built datasets for three special cases: inextensible beams-dataset 2-inextensible and inflexible beams-dataset 3-and perfect pivotsdataset 4. In datasets 2 and 3, a factor 1000 for increasing the rigidified stiffnesses is considered. Dataset 4 accounts for a special case extensively described in [29]: the pivots are designed in such a way that the corresponding torsional rigidities are practically zero and, therefore, the corresponding stiffness parameter c of each pivot is set to zero. We consider the case of a clamped-clamped p-beam, therefore constraints are set on the nodes A, B, C and D, see

Dataset 1
Making use of the stiffnesses in dataset 1, see Table 1, we first compute the natural frequencies and modes of the p-beam. The calculation is based on the mass matrix M and on the stiffness matrix K computed in the reference configuration, see Sect. 2. Figure 3 reports the first eight natural modes along with the associated natural frequencies. We identify six bending modes-1, 2, 3, 5, 7 and 8-and two extensional modes-4 and 6. In the study of the dynamical behavior of pantographic beams, we limit ourselves to consider harmonic excitations, i.e. excitations representable as p(t) =p sin(2π f t) wherep is a constant vector collecting the generalized forces, f the excitation frequency and t the time. In general, the p-beam displacements depend on the magnitude ofp and on the excitation frequency f .
Fixing the magnitude ofp and varying f , we can compute numerically-using the integration scheme depicted in Sect. 3-the so-called frequency-response curve, see We remark that the computed response amplitude also includes its transient part.
We notice some remarkable peculiarities in these curves. The first concerns the amplitude of the longitudinal displacement- Fig. 4a-which for δ ≈ 1.7 shows a peak. Such a peak may be triggered by mode interaction. Indeed, the considered load, because of symmetry arguments, should entail only deformations corresponding to the first mode, which is, see Fig. 3, a bending mode. The second peculiarity concerns the amplitude of the transversal displacement, which shows a jump for δ ≈ 1.4 − 1.6- Fig. 4b: the amplitude drops of about 20% when the load frequency increases and there is a sharp increment of about 20% when the frequency of the excitation decreases. The third peculiarity concerns the increment of about 8% in the transversal response amplitude, with respect to the trend curve, for δ ≈ 2.7. Most likely, for such a frequency, we have a secondary resonance. Finally, contrarily to the linear case, we notice the absence of a resonance in correspondence of δ = 0. Figure 5a reports the response, intended as the evolution in time of the transversal displacement v M of the p-beam midpoint M, against the time t for the load frequency f = 0.95238 Hz, which corresponds to δ = 1.674. Transversal velocity is reported for the same problem in the phase plane Fig. 5b, while Fig. 5c shows the Fast Fourier Transform of the transversal displacement v M in Fig. 5a. It shows three peaks corresponding to the first three natural frequencies. Figure 6 reports a stroboscopic images of the forced oscillations of the p-beam (displacements are amplified in order to make the image more readable). 12  Having completed the description of the numerical simulation for the symmetrical load, we now consider the case of unsymmetrical load, namely a transversal load acting on the node N . Figure 7 reports the corresponding frequency-response curves, namely the dimensionless amplitude A/A t max of longitudinal- Fig. 7a-and transversal-Fig. 7b-displacements of the p-beam midpoint M versus the detuning parameter In this case, we notice that there is a less significant longitudinal oscillation amplitude within the whole range of considered load frequencies. Deviations from the curve trend are limited. The transversal oscillation amplitude does not present a jump as in the case of the symmetrical load, but there is a remarkable deviation from the trend in the neighborhood of δ ≈ 2.7. Figure 8 shows the response of the p-beam in the unsymmetrical load case for an excitation frequency f = 1.3333 (δ ≈ 2.7). In particular, it reports the evolution of longitudinal, Fig. 5a, and transversal, Fig. 8b, displacement, the evolution of longitudinal, Fig. 8c, and transversal, Fig. 8d, velocity in the phase plane, as well as the frequency-amplitude plots for longitudinal, Fig. 8e, and transversal, Fig. 8f, motions. Figure 9 shows a stroboscopic image of the p-beam oscillations triggered by the unsymmetrical load with frequency f = 1.3333 Hz. 13

Dataset 2
Using the dataset 2, which accounts for the case of inextensible beams, we first compute again natural frequencies and modes. Figure 10 reports the first eight natural modes (red and grey colors show natural modes and the rest configuration, respectively) along with the corresponding natural frequencies. We distinguish again six bending modes-1, 2, 3, 5, 7 and 8-and two extensional modes-4 and 6.  Natural modes are indistinguishable from those deriving from the analysis of the dataset 1. Differences between corresponding frequencies are very limited. For instance, the maximum relative difference between the corresponding frequencies is equal to 1.18%. Numerical simulations analogous to those developed for the dataset 1 do not show significant differences and for this reason are not reported. We remark that the choice of modeling pantographic structures neglecting the extensibility of the beams appears reasonable.

Dataset 3
As well as for datasets 1 and 2, also for the dataset 3 we compute the natural frequencies and modes. Figure 11 reports the first eight natural modes along with the corresponding natural frequencies. In this case, we distinguish four bending modes-2, 4, 6, and 8-and four extensional modes-1, 3, 5, and 7. We notice the differences with respect to datasets 1 and 2. Indeed, we have now four bending modes and four extensional modes but, above all, we notice that the first mode is now an extensional mode. Also natural frequencies present remarkable differences (the first is about 10 times bigger). Figure 12 reports the frequency-response curves for dataset 3 by consideringp = 10. The analysis of this case shows, at least in the range of analyzed frequencies, very flat curves. In particular, the amplitude of the longitudinal oscillations is practically always near to zero while that of the transversal oscillations is practically constant in the whole range of the analyzed frequencies.

Dataset 4
It is the case of a p-beam with perfect pivots. First of all, we compute again the natural frequencies and modes, see Fig. 13 which reports the first eight natural modes along with the corresponding natural frequencies. We distinguish four bending modes-2, 4, 6, and 8-and four extensional modes-1, 3, 5, and 7. We remark that  Figure 14 reports frequency-response curves for longitudinal (a) and transversal (b) motions of the p-beam midpoint in the case of a transversal harmonic load on the p-beam midpoint. As transversal load amplitude, we assumep = 0.02 N. In this case, we observe again that for δ ≈ 4.2 there is a remarkable longitudinal oscillation amplitude. In the analyzed frequency range, dimensionless transversal oscillation amplitudes are practically linear-increasing the load frequency from 0.92 to 1-with respect to the detuning parameter δ.
Analogously to what has been done for the dataset 1, Fig. 15 reports the response of the p-beam midpoint for the symmetrical load with frequency f = 0.5 Hz. More specifically, the plot reports the midpoint transversal displacement evolution, Fig. 15a, the midpoint transversal velocity in the phase plane, Fig. 15b, and the dimensionless amplitude A of the transversal displacement versus the frequency ϕ in Fig. 15c. Figure 16 shows a stroboscopic image of the p-beam oscillations triggered by the symmetrical load with frequency f = 0.5 Hz. 14 Similarly to what has been done for the dataset 1, Fig. 17 reports the frequency-response curves for longitudinal (a) and transversal (b) motions of the p-beam midpoint in the case of a transversal harmonic load on the node N (unsymmetrical load) and forp = 0.005 N. Also in this case, we notice that, in the range of analyzed frequencies, the dimensionless transversal oscillation amplitude is practically linear. What is more remarkable is the longitudinal oscillation amplitude of the p-beam midpoint which varies between 9.5% and 11.5% of the transversal oscillation amplitude of the same point. Figure 18 reports the response of the p-beam midpoint subjected to an unsymmetrical load with frequency f = 0.5 Hz. In particular, Fig. 18a, b report the longitudinal and transversal displacement evolutions, respectively; Fig. 18c, d report the longitudinal and transversal velocities in the phase plane, respectively, and  f report the dimensionless amplitude of the longitudinal and transversal displacements, respectively, versus the frequency ϕ as computed by their Fast Fourier Transform. Figure 19 shows a stroboscopic image of the p-beam oscillation triggered by the unsymmetrical load with frequency f = 0.5 Hz. 15

Concluding remarks and future challenges
In this work, we have dealt with the dynamical response in large displacement regime of pantographic beams subjected to harmonic loads. In particular, we analyzed four different stiffness datasets and computing, among the other things, the corresponding frequency-response curves by means of a stepwise numerical approach. Differences with what is expected in small displacement regime were thoroughly discussed in order to build a solid base foundation for successive explorations. The simple mechanical model and the stepwise analysis strategy used in the present investigation have proven their capability to provide a first set of results. The peculiar results obtained within this work could help all researchers concerned with the dynamics of pantographic structures and materials (see [45]), and with metamaterials in general-see, e.g., the very recent work [46]-to face many open or still unexplored problems.
We believe that future research lines of this work shall include: (i) the enhancement of the mechanical model used for the p-beam modeling such as proposed in [47]; furthermore, the three-dimensional discrete beam formulation proposed in [48] could be very fruitful to consider out-of-plane instability phenomena; (ii) the generalization of the stepwise strategy in order to handle the case of non-conservative forces such as those deriving from damping; (iii) the formulation of continuum, hence synthetic, models capable to forecast the main mechanical behaviors of pantographic beams, see, e.g., [7]; (iv) the extension of the utilized discrete mechanical approach to micropolar media [49], thin plates with surface stresses [50] and strongly anisotropic surface elasticity [51]; (v) the application of the utilized time-integration scheme to continuous problems discretized, for instance, by means of B-spline and NURBS, see, e.g., [52] or for models dealing with granular media, e.g., [53,54].