Parameter identification of a second-gradient model for the description of pantographic structures in dynamic regime

Pantographic structures are examples of metamaterials with such a microstructure that higher-gradient terms’ role is increased in the mechanical response. In this work, we aim for validating parameters of a reduced-order model for a pantographic structure. Experimental tests are carried out by applying forced oscillation to 3D-printed specimens for a range of frequencies. A second-gradient coarse-grained nonlinear model is utilized for obtaining a homogenized 2D description of the pantographic structure. By inverse analysis and through an automatized optimization algorithm, the parameters of the model are identified for the corresponding pantographic structure. By comparing the displacement plots, the performance of the model and the identified parameters are assessed for dynamic regime. Qualitative and quantitative analyses for different frequency ranges are performed. A good agreement is present far away from the eigenfrequencies. The discrepancies near the eigenfrequencies are a possible indication of the significance of higher-order inertia in the model.


Introduction
Metamaterials are referred to as a group of engineered materials that possess peculiar properties, which are not commonly observed in regularly used materials [1][2][3][4]. The properties that we observe at the length scale of the sample (macro-scale) are inherited from a designed substructure of the material at a lower scale. Thus, these multi-scale materials are also called micro-architectured materials [5][6][7]. Metamaterials reveal specific behaviors, particularly under dynamic loads. These unusual mechanical behaviors at the macro-scale are due to the interactions of the underlying microstructure of the metamaterial. For example, we may mention the so-called locally resonant acoustic metamaterials [8], or the lattice fiber materials [9] as well as inclusions' effect due to their moduli mismatch [10]. The role of material symmetry considering static and dynamic properties of materials is discussed in [11].
In principle, it is possible to describe the mechanical behavior of the metamaterials with a complex substructure by means of the standard first-gradient theory of elasticity. The computational cost owing to a mesh with such a small length scale is excessively high [12,13]. It is well known that homogenization [14] adds additional terms to the energy definition effected by the microstructure. A homogenized continuum of this multi-scale material is simulated by a far lower computational cost. From a numerical point of view, this homogenized model is in favor. Yet, these additional parameters have to be determined [15][16][17][18].
For generating an adequate homogenized continuum, there exist different techniques. An approach for modeling metamaterials is employing the generalized continuum mechanics [19][20][21][22]. Generalization of the classical continuum mechanics makes it possible to model the mechanical properties of metamaterials, accurately [23][24][25][26][27]. One possible generalization is the strain-gradient theory, considering the second gradient of the displacement field as an independent kinematical descriptor of the system [28,29]. Consequently, the dynamics described by such a model will also depend on additional parameters. The parameter values are related to the properties of the microstructure [30,31]. The validity range of the strain-gradient elasticity is shown to be sufficiently large in terms of frequency and wavelength [32]. We emphasize that the generalization adds parameters to the deformation energy as well as to the kinetic energy. Particularly, the higher-order inertial terms are not yet understood from an experimental perspective.
An example of metamaterial is the so-called pantographic structure characterized by such a substructure increasing the significance of the aforementioned additional parameters [33][34][35][36][37]. The pantographic structure is capable of undergoing large tensile deformation while remaining in the elastic regime [38], and it has a very favorable strength-to-weight ratio [39]. In this structure, two orthogonal groups of parallel beams are mutually connected by means of pivots. The pivots can be of two types: either perfect pivots (hinges) [40,41], or elastic cylinders [42][43][44][45]. Herein, we utilize pivots in the form of deformable cylinders. In this structure, the relative rotation of the beams stores deformation energy in the pivots. We refer to Fig. 1 for the geometry of such a structure. The fabrication of pantographic structures is possible by 3D printing or stereolithography technologies [46][47][48][49][50][51].
In a quasi-static loading regime, where inertial terms are negligible, deviation from classical elasticity has been analyzed in ample studies [52][53][54][55]. In the case of a harmonic loading, especially in frequencies where the inertial terms dominate the mechanical response, experimental observations are more recent [56,57]. In [58], a linear second-gradient model is used to model the dynamic behavior of the pantographic structure under harmonic loading.
We hypothetically assume that parameter estimation of metamaterials may be decomposed into two parts. Since the deformation and kinetic energies are independent, we may determine the parameters in each part, separately. For the former, quasi-static loading condition makes kinetic energy negligible. For the latter, loading around eigenfrequencies lets inertial terms dominate. We focus on the former herein and use the formalism in [40], where an automatized optimization problem is utilized to determine the parameters of the mathematical models of metamaterials under static loading. The main goal of the underlying work is to validate this numerical identification procedure and to detect its validity range by analyzing the accuracy of the model at different frequencies. In this paper, a pantographic structure undergoing harmonic loading is considered. The parameters of a macro-scale second-gradient model are identified for the structure. Furthermore, experimental tests are carried out in order to corroborate the numerical results and reveal the limitations of the assumptions used in the model. This paper is structured in the following way. First, in Sect. 2, the geometry and material of the investigated structure are presented. Then, in Sect. 3, the computational method is described, including the formulation and implementation of the strain-gradient model, along with the numerical identification of its parameters. Next, in Sect. 4, the experimental rig and the test procedure are sketched. Finally, the results and error analyses are shown in Sect. 5.

Geometry and material
The mechanical system considered in this paper is a planar pantographic structure with deformable cylindrical pivots. The geometrical dimensions of the structure are shown in Table 1 and Fig. 2.
The test specimens of the pantographic structure are made by selective laser sintering (SLS) 3D printer, Formiga P100, and the material is Polyamide EOS PA 2200. In [59], the mean value of Young's modulus for the laser-sintered polyamide 2200, obtained through a bending test, is reported as 1614 MPa, and in [60], a Poisson's ratio of 0.394 is calculated by performing tensile tests. Here, we need the "surface" mass density since the structure is modeled as a 2D system as it is described in Sect. 3.1. The structure's surface mass density, ρ 0 = 1.266 kg/m 2 , is measured as follows. The mass of a pantographic structure sample is measured by a scale. Then, the mass is divided by the surface area of the structure, i.e., L × l.

Numerical simulation
For the numerical analysis of the pantographic structure under dynamic loading, we make use of a secondgradient model [61,62] suitable for planar pantographic structures. Then, the parameters of this model are determined for the homogenized continuum by virtue of the aforementioned optimization algorithm.

Strain-gradient energy model
In the following, the formulation of a macro-scale nonlinear model [61] for planar pantographic structure is explained briefly. This model is obtained by means of a homogenization procedure from micro-scale.
In the macro-scale model, a homogeneous rectangular plate is modeled instead of the detailed structure, as shown in Fig. 2. Assume that χ maps the initial position X of every material point in the 2D rectangular region to its current position x in R 2 through the displacement u, where (D 1 , D 2 ) are the unit base vectors of a Cartesian coordinate system oriented along the direction of the beams as shown in Fig. 2. The unit vectors e 1 and e 2 are aligned in the direction of the deforming N. Shekarchizadeh et al. ZAMP Fig. 2. The geometry depicted on the pantographic structure, see Table 1 for the geometric dimensions beams, where || · || means the L 2 norm and ∇χ is the deformation gradient, The "stretch of beams," ε, and the "beam geodesic curvature," κ, are defined as where In the vectors c 1 and c 2 , the second-gradient terms of displacement are present. The angle between the beams of the two families is initially 90 degrees, while during the deformation it is given by e 1 and e 2 . The change in the beams angle is the shear distortion, γ = arcsin(e 1 · e 2 ).
The strain energy density, W M , is defined as which accounts for elongation and geodesic bending of the beams and the torsion of the connecting pivots. The positive and constant constitutive parameters K e , K g , and K s are the stretching, geodesic bending, and shear stiffnesses, respectively. The model admits strain-gradient terms in the formulation of the strain energy and therefore is capable of capturing the effects of the complex substructure. The stretching and bending stiffnesses of the beams and the shear stiffness due to the connecting pivots are taken into account.

Weak form
For obtaining the governing equations of the structure, leading to the weak form by the so-called variational formulation, we exploit the least action principle by following [63]. As a matter of fact, for reversible (non-dissipative) systems, the variation of the scalar action functional vanishes for an arbitrary test function δu as follows: Under the assumption that kinetic energy has no additional terms because of the microstructure, we assume the action functional for this problem as where t 0 and t f are the initial and final time of the interval of interest, respectively. B 0 is the continuum body in the reference configuration, and ∂B N 0 is the part of its surface where the traction,t i , is applied. u i is the displacement,u i is the time derivative of the displacement, ρ 0 is the mass density at the reference configuration, and f i is a specific (per unit mass) dead load. Latin indices denote space coordinates. Here and henceforth, we employ Einstein's summation convention over repeated indices.
The first integral in Eq. (10) accounts for three different forms of energy in the system: The first term, 1 2 ρ 0uiui , is the kinetic energy, the second term, W M (ε, κ, γ), is the strain energy, and the third term, ρ 0 f i u i , is the potential energy of the external load, which can be identified with the weight. The second integral is regarding the work done on the boundary of the system. We emphasize that we have used the following assumptions in this reduced-order model: • Kinetic energy is equal to the classical theory; in other words, higher-order inertial terms are neglected. • Strain energy uses the substructure information such that "only" three additional parameters are necessary to accurately describe the response. • Mass density is the mean average analogous to a porous structure.
Since in the application with a given loading, the effect of gravitational force on the deformation is negligible compared to the applied force, we ignore the force term f i . By taking the variation of the action using Eq. (9) by considering isochronous motions [64,65] and by integration by parts on rate of test function, we get the governing equation of the elasticity problem as follows: where andü i is the acceleration. The inertial term, ρ 0üi , has an important role in this problem as we are dealing with dynamic loading. The primitive variable of the partial differential equation (PDE) in Eq. (11) is the displacement u i ; hence, by solving the PDE, the displacement is obtained within the whole domain of the model. In Eq. (11), the derivation of the strain energy with respect to the second gradient of displacement is neglected. This approximation is made in order to have a more efficient code.

Finite element method based simulations
The numerical problem is implemented in FEniCS [66], which is an open-source platform capable of automatic solution of partial differential equations by the finite element method (FEM). We refer to [67] for engineering applications by using FEniCS. A homogeneous 2D rectangular plate is created together with the boundary conditions for simulating the dynamic loading. As shown in Fig. 3, the bottom edge is fixed in the xand y-directions, the top edge is fixed in the x-direction, and a sinusoidal displacement is prescribed on the top edge in the y-direction as where A is the amplitude of loading and f is the frequency of loading over time, t. This sinusoidal function matches the harmonic loading that is applied in the experimental tests. We demonstrate computations for different frequencies of the dynamic loading ranging from 40 to 200 Hz in accordance with the experimental procedure.
For the numerical implementation of the weak form in Eq. (11), the equations need to be discretized in time and space. For the space discretization, we use Hilbertian Sobolev space, H 2 , spanned by quadratic shape functions with compact support; we refer to [68] for details of the finite element method. Simply stated, we use standard polynomial finite elements constructing a vector space where on the so-called Dirichlet boundaries ∂B D the displacement is given. The rectangular domain is meshed with quadratic Lagrange elements for space discretization. For time discretization of the problem, we employ the so-called backward Euler method with a constant time step, Δt, for approximating the second time derivative of displacement as follows: where u 0 i and u 00 i are the displacements calculated, respectively, in one and two steps before the current step.
The derivative of W M is taken symbolically. We refer to [66,69] for an overview of this technique. In FEniCS, for taking the derivative of W M with respect to u i,j , first, u i,j is declared as a variable. Therefore, the derivative of some function can be taken with respect to it.

Material parameter identification
We aim for determining the metamaterial's parameters, K e , K g , and K s , which are manifested by the substructure, herein the pantographic structure as in Fig. 1. To this end, we perform an inverse analysis based on an optimization algorithm following the formalism presented in [70]. The overall methodology is illustrated in Fig. 4.
At first, a micro-model is created and simulated that is a 3D model of the structure by considering all the substructure details. A static bias extension test is simulated by applying a displacement boundary condition on the smaller side of the rectangular sample. The classical elasticity theory is used to solve the problem. After an h-convergence analysis, an FEM mesh of 774,000 degrees of freedom has been selected. The maximum extension is applied in eight progressive steps, and at every step, the strain energy of the whole structure is calculated by the micro-scale model.
Then, an initial guess is made for the parameters of the macro-scale model, and therefore, the strain energy is also calculated by the macro-scale model for the same eight loading steps. The parameters of the homogenized model are determined by enforcing the same energy values in the micro-scale and homogenized model. This process is carried out by setting up an automatized optimization problem that minimizes the difference between the energy values obtained from the micro-and macro-scale models. The employed optimization method utilizes the Trust Region Reflective algorithm [71,72] implemented in Scipy [73] and coded in the Python language.

Experimental tests
The experimental tests have been performed at the Marcus Wallenberg Laboratory for Sound and Vibration Research, KTH Royal Institute of Technology. The experimental rig and the experimental procedure are described in this section. Figure 5 shows the experimental rig. The sample is clamped at one of its short edges-the one that is closest to the ground-to a magnesium table, assumed infinitely stiff regarding the pantographic structure. This configuration realizes in practice the zero-displacement boundary condition imposed in the model. The other short edge is connected to a Brüel & Kjaer type 4809 electrodynamic shaker, suspended in a purpose-built cradle via soft rubber springs. The height of the shaker is adjusted in order to slightly prestretch the specimen (3 mm for the structure of 213 mm long), a necessary measure to avoid the occurrence of buckling instabilities. The long sides of the specimen are perpendicular to the magnesium table and aligned in the direction in which the shaker vibrates. This motion effectively imposes an excitation that is almost along the vertical direction.

Experimental rig
The shaker is driven by a Brüel & Kjaer type 2718 power amplifier, in turn, connected to a NI-9263 output module sitting in a cDAQ-9178 chassis. These devices are connected to a computer via USB and controlled using the Matlab Data Acquisition toolbox. A laser Doppler vibrometer is used to measure the velocity of the moving part of the shaker. The moving part of the shaker is rigidly connected to the short edge of the sample, hence exciting it. The vibrometer is aimed at the moving part via a mirror. The measured velocity over time is compared with the displacement measured using digital image correlation (DIC), and its sole purpose is to assess the accuracy of the camera-based measurements.

Measurement routine
For each frequency within a given range, the measurement routine follows these steps: 1. Four 250 W spotlights are turned on, allowing a short camera exposure (≈ 250 µs) without motion blur. 2. A sinusoidal signal at the given frequency is sent to the power amplifier. Consequently, the sample starts moving. 3. After approximately 2.5 s with at least 100 oscillations, the deformation state is considered at steady state and a trigger is sent to the cameras. Next, the acquisition of 400 stereo frames at a sampling frequency of 2048 Hz is performed. The sampling frequency ensures supersampling for anti-aliasing, while the record length is nine times the minimum imposed excitation frequency within the range. Finally, the video is transferred to persistent storage and the spotlights are turned off.

Data acquisition and digital image correlation (DIC)
The stereo frames acquired are analyzed using digital image correlation (DIC), an optical technique that relates correlation between subsets of pixel intensities to displacements [74]. A variety of techniques and implementations of DIC exist (see [75][76][77][78][79]); herein, we employed a commercial white-light 3D DIC solution (DaVis 8.4) based on the Lucas-Kanade algorithm for optical flow estimation [80][81][82]. The cameras (two Phantom v1612 cameras at resolution of 1280 × 800) are paired to macro lenses with a nominal focal length of 100 mm. The system is calibrated using a two-level calibration plate: A pinhole camera model is employed, accounting for both radial and tangential distortion, yielding an RMS value of the reprojection error smaller than 0.2 px. Subpixel interpolation with a sixth-order spline is used for the subset matching procedure on observation windows of 11 × 11 px with a 2 px step. The displacement measurement has been mapped onto and retrieved from 184 × 560 equidistant locations for the whole pantographic structure, which we call the DIC mesh. The measured amplitudes of the applied sinusoidal loading, corresponding to Eq. (13), for different frequencies are presented in Table 2.

Linear time-invariant data analysis
The data analysis assumes that the system is linear time invariant (LTI). Time invariance results in a periodic output in the case of a periodic input. Linearity denotes that the input and output signals are linearly related by means of the convolution with the impulse response of the system. The main consequence of this assumption is that a sinusoidal imposed displacement on the system will result in a sinusoidal displacement, generally characterized by different amplitudes and phases.
In order to make a comparison between numerical results and experimental observations, we will compare the results from FEM and DIC, respectively. Indeed, both methods provide the displacements in time for each point of the DIC and FEM meshes which, according to the LTI hypothesis, can be modeled in terms of the general sinusoidal function, where A is the amplitude of the signal, f is the frequency of the displacement, ϕ is the phase, C is a bias, and the bold notation refers to the two independent directions given along the short and long sides of the samples in which the oscillation is decomposed. A nonlinear fitting procedure is carried out on the displacements in time. The numerical and experimental displacements for each node of the FEM and DIC mesh are fitted with Eq. (17). Therefore, the time behavior of each point of the system is compressed in the set of parameters (A, f , ϕ, C). The values of these parameters are plotted on the geometry of the system for both experimental and numerical results. Hence, a meaningful comparison is obtained between experiments and simulations of the behavior of the system during the measured time.

Results and discussion
In this section, the identified parameters of the second-gradient model, the numerical simulations, and the experimental results of the dynamic loading on the pantographic structure are presented. Furthermore, the FEM mesh dependency of the numerical model is investigated. At last, the accuracy of the macro-scale model in predicting the mechanical behavior of the structure is quantified by an error measure.
As a result of the inverse analysis, the parameters of the second-gradient model are determined for the studied pantographic structure as presented in Table 3. The displacement plots from numerical computations and the experimental tests for 40, 80, and 140 Hz are presented in Figs. 6, 7, and 8, respectively, where |u x | is the modulus of the (horizontal) displacement in the x-direction and |u y | is the modulus of the (vertical) displacement in the y-direction. The plotted displacements are for the time instant in which the maximum tensile excitation is applied in the fifth period of harmonic loading.
The convergence of the macro-scale model for dealing with dynamic loading is studied here. The vertical displacement (u y ) of the central point of the rectangular sample is computed for different degrees of freedom for the 40 Hz frequency. The h-convergence results with a posteriori analysis for 40 Hz are given in Table 4. The final computations and the results presented above are carried out by using the FEM mesh with 308,482 degrees of freedom; see Table 5 for error analysis for all frequencies with this FEM mesh.
For evaluating the accuracy of the simulations in predicting the experimental results, a mean error percentage is defined as where || · || means the L 2 norm. In Eq. (18), for each nodal displacement in the FEM mesh, u i comp , the corresponding displacement value from the DIC (experimental) mesh, u i exp , has been acquired and compared. The total number of nodes in the sample is n, and A is the amplitude of the harmonic loading. In Eq. (18), the mean value of an error is calculated over all the FEM mesh nodes of the structure and then divided by the amplitude of the excitation to normalize the error value for all the frequencies. The error values are reported in Table 6.
As apparent from the results in Table 6, for frequencies 40, 60, and 80 Hz, the results are adequate under admitting experimental errors introduced in measurements. These frequencies are not close to the eigenfrequencies of the structure. Starting from 100 Hz, the simulation results start deviating significantly from experiments. This deviation is understood that the taken assumptions prevent an appropriate modeling of the response. Two assumptions are of importance: (I) using a planar model, i.e., ignoring the out-of-plane displacements, and (II) not including the higher-order inertial terms in the model.
Since the utilized mathematical model is for planar pantographic structures on the xy-plane, the numerical computations are in two dimensions. The DIC data using two cameras include the out-ofplane displacements as well (in the z-direction) [83]. For an understanding of how large the out-of-plane displacements are, and in order to investigate whether a planar model is a proper choice, we define a measure as where z i is the out-of-plane displacements at the ith node of the DIC mesh. In other words, z-measure is the mean value of modulus of displacement in the z-direction divided by the amplitude of the harmonic excitation. The values of z-measure for all the frequencies are compiled in Table 7, where we realize a noticeably greater value for the 160 Hz test.
The eigenfrequency is expected to be around 160 Hz such that the increased gain in response (ratio of output amplitude to input amplitude) causes an out-of-plane deformation. This deformation may be seen as buckling but not because of an instability rather due to an increase of the mechanical response near the eigenfrequency. Generally speaking, the out-of-plane deformation is simply an artifact of the substructure with so-called non-perfect pivots. Over the thickness, the substructure fails to be symmetric such that a buckling happens as a possible deformation direction in the case of the increased gain. There is no instability or snap-back phenomena observed such that we understand that the deformation is nonplanar near the eigenfrequency. Indeed, this phenomenon justifies the significant discrepancy between the computational and experimental results for the 160 Hz loading. Therefore, as expected, the planar model is not a proper candidate for modeling the problem around the eigenfrequencies.
Even at frequencies of 100 Hz and above, there is a significant deviation between the experimental and computational results. Although the out-of-plane displacement in Table 7 shows that a 2D modeling may  be adequate, this deviation is related to the second assumption in the computational model by ignoring the higher-order inertial terms. From a theoretical point of view, substructure introduces additional dynamical terms [84]. Modeling these terms is known in the literature as well [85][86][87].
Although higher-order inertial terms are existing in theory, their experimental characterization is yet unresolved. The discrepancy in Table 6 may be seen for 100, 120, 140, 160, and 180 Hz. Under the hypothesis that starting from 200 Hz, the error is going to be at the level of 40, 60, and 80 Hz, then we may claim that the higher-order inertial terms are of importance between 100 and 180 Hz. This claim is easily justified with the analogy that the kinetic energy-related inertial term dominates the system response around the eigenfrequency. We examine this possible explanation by changing the weak form in Eq. (16) as follows: where d is the added inertial term. Here, the backward Euler method is employed for time discretization of the second time derivative of the displacement gradient. The term d is a length scale (unit of length), as suggested in [88]. By using fictitious parameters for d in Eq. (20), we repeat the numerical computation of the 140 Hz test. We set different values for d = 0.005, 0.01, and 0.02 m in order to observe the system behavior. The results are plotted in Figs. 9 and 10 . The plots are the modulus of displacement at the instant in which the excitation is at its maximum tensile value in the fifth period of loading.
By comparing Fig. 9a with Fig. 9c and also Fig. 9d with Fig. 9f, we detect that the displacement in the numerical response is smaller in value, and also there is a shift in the location where the maximum displacement is happening. Qualitatively, the same effect has been acquired by including the extra inertial terms in the model. By comparing Fig. 9b with Fig. 9a and also Fig. 9e with Fig. 9d, it is visible that the introduced inertia has caused changes in both the amplitude of displacement and the location of maximum displacement. This effect is in a way explaining the discrepancies shown in Fig. 8. This result shows that the model with higher-order inertia may be of importance for simulating problems with dynamic loading, especially when close to the eigenfrequencies.
Again, we stress that the values chosen for the higher-order inertial terms are chosen arbitrarily in order to comprehend their role in the deformation pattern, and we circumvent ourselves from determining these parameters. In fact, the role of such complicated inertia in a process is not yet fully realized. As shown in Fig. 10, where d values are varied, the system response is highly nonlinear. There is a shift and change in the displacement form at the same time. This result is adequate since the parameter controls second time and second space derivative of the displacement in the energy formulation. Yet at the same time, the characterization of this parameter may be impossible by using one set of data. Hence, it is not clear how to determine or get an estimate of these values. Further studies need to be carried out in this direction.

Conclusion
In this paper, the parameters of a planar strain-gradient model for pantographic structures are determined for dynamic loading, and the accuracy of the numerical method is assessed. The parameters of the model are estimated by virtue of an automatized optimization algorithm. Moreover, experimental tests are performed on 3D-printed specimens of a pantographic structure. The structures undergo forced oscillations under different frequencies. By comparing the experimental and numerical results, it is observed that the computations predicted the displacement field, with acceptable error, in the frequencies 40 to 80 Hz which are far from the eigenfrequencies. In contrast, in higher frequencies, up to 180 Hz, discrepancies were present in the results, which may be traced back to not including higher-order inertial terms. To summarize, the following is achieved in this paper: • Purely computational parameter determination has been utilized for planar pantographic structures with deformable pivots. • Experiments and computations show adequately matching results far from the eigenfrequencies.
• Noticeable out-of-plane displacement is possible around the eigenfrequencies.
• Implementation of higher-order inertia in the computations results in favorable results toward predicting the behavior of system close to the eigenfrequencies, however, the value of the inertial terms has not been determined herein, and we encourage further studies on this matter.