Article Mathematics and Mechanics of Solids 2022, Vol. 27(10) 2201–2217 � The Author(s) 2022 Article reuse guidelines: sagepub.com/journals-permissions DOI: 10.1177/10812865221107623 journals.sagepub.com/home/mms Experimental analysis, discrete modeling and parameter optimization of SLS-printed bi-pantographic structures Jonas Harsch Institute for Nonlinear Mechanics, University of Stuttgart, Stuttgart, Germany Gregor Ganzosch Technische Universität Berlin, Continuum Mechanics and Materials Theory Group, Institute of Mechanics, Berlin, Germany Emilio Barchiesi Department of Architecture, Design and Urban planning (DADU), University of Sassari, Sassari, Italy Alessandro Ciallella International Research Center on Mathematics and Mechanics of Complex Systems (M&MOCS), University of L’Aquila, L’Aquila, Italy Simon R Eugster Institute for Nonlinear Mechanics, University of Stuttgart, Stuttgart, Germany Received 25 March 2022; accepted 31 May 2022 Abstract Making use of experimental data for bias extension, shearing, and point-load tests in large deformation regime for rectan- gular and square bi-pantographic specimens, we perform a numerical identification to fit the a priori parameters of a pla- nar discrete spring model. The main objective of the work is to develop an automatized optimization process based on the Nelder–Mead simplex algorithm for identifying the constitutive parameters of discrete modeling of bi-pantographic structures, as well as assessing its descriptiveness and predictive capacity. The analysis allows to conclude that there exists a single set of parameters for the adopted discrete modeling such that it is descriptive and predictive for several different tests and for a wide range of deformations. Corresponding author: Jonas Harsch, Institute for Nonlinear Mechanics, University of Stuttgart, Pfaffenwaldring 9, 70569 Stuttgart, Germany. Email: harsch@inm.uni-stuttgart.de https://uk.sagepub.com/en-gb/journals-permissions https://doi.dox.org/10.1177/10812865221107623 journals.sagepub.com/home/mms Keywords Bi pantographic structures, experiments, discrete modeling, parameter optimization, large deformations 1. Introduction This paper deals with a generalization of pantographic fabrics [1,2]. Pantographic fabrics are actual fab- rics which are made up of two orthogonal families of parallel beams. Beams of different families are con- nected at their intersection points by means of perfect hinges or, alternatively, by means of elastic hinges, i.e., hinges and rotational springs opposing to the relative rotation of beams concurring to the same hinge. Within each family, beams are at a given distance between each other. Such a distance might be depending on the beam pair or be constant. If each family of beams contains the same number of beams, having the same constant distance one from each other, and the beams have all the same dimen- sions and physical properties, then the pantographic fabric is said to be balanced [3]. Pantographic fab- rics were designed and then realized physically by arranging the two families of beams onto two different parallel planes. The design and realization of connections between two beams have been accomplished by means of cylindrical pivots acting ideally as in-plane elastic hinges, i.e., hinges with tor- sional springs [4,5], and by means of actual hinges which, however, are prone to suffering significant non-elastic phenomena and hence pose complex challenges in manufacturing engineering. While the rea- lization of cylindrical pivots is less challenging from the manufacturing point of view—there are, how- ever, issues related to the plane in which scanning tracks lie and heat treatments for metallic specimens—it has been proved by experimental and computational investigations that the behavior of cylindrical pivots involves complex deformation modes and, particularly, significant non-torsional deformation modes in large deformation. Beams in pantographic fabrics play the role of bending ele- ments. Their utilization in pantographic fabrics has been conceived to confer to the fabric, when regarded as a continuum, second gradient effects in plane due to the dependence of the deformation energy upon the geodesic bending of the material lines. However, the presence of beams does not only entail these desired effects. The fact that in the physical reality beams are arranged on two different planes makes them interact mutually, in a way such that torsional deformation modes [6], as well as in- plane and out-of-plane buckling [7–9], are triggered, even when problems with in-plane boundary condi- tions are considered [10]. In the last years, many research efforts published in the scientific literature have dealt with variations and enhancements of the design of pantographic fabrics. As mentioned before, the design and manufac- turing of perfect hinges in pantographic fabrics have been attempted in the context of three-dimensional (3D) printing using both metallic and polymeric raw material (see in this regard the works [11,12]). It is worth to mention in this introduction that also non-orthogonal families of fibers [13], as well as non- straight fibers [14], have been considered in the scientific literature. The need for performing and under- standing experimental analyses on pantographic fabrics has stimulated many studies in Digital Image Correlation (DIC), manufacturing techniques, and material characterization [15–17]. It is worth to men- tion here that pantographic specimens have been realized also at extremely low scales, including the nano-scale [18–20]. Pantographic materials belong to a class of materials which includes not only pantographic fabrics but, more generally, materials with a microstructure based on the pantographic motif, i.e., a mechanism which is well known from everyday life (pantographic mirrors, expanding fences, scissor lifts, etc.), which is characterized by a zero-energy accordion-like homogeneous extension/compression deforma- tion mode. For example, pantographic beams—namely, strips of pantographic fabrics consisting of the periodic repetition in one direction of a single pantographic cell—have been proposed as an example of complete second gradient beam [21–24]. In this paper, a particular kind of metamaterials is addressed, namely, the so-called Bi-Pantographic Structures (Bi-PSs). Bi-PSs were conceived as a generalization of pantographic fabrics and consist of pantographic fabrics whose beams have in turn a pantographic microstructure and are, actually, pantographic beams (see Figure 1). They were first introduced in the paper [25]. One of the main advantages of additive manufacturing lies in the possibility of easily producing com- plex 3D geometries at, practically, no additional cost. For this reason, the development of 3D-printing techniques witnessed in the last two decades has stimulated research in the field of mechanical 2202 Mathematics and Mechanics of Solids 27(10) metamaterials. These materials assume special properties depending on their microstructure. However, while additive manufacturing is experiencing great progress in the number of raw materials that can be used, see as an instance the recent works on the 3D printing of concrete structures, the main raw materi- als used in additive manufacturing are still polymeric or metallic. When polymeric and metallic materi- als are considered, additive manufacturing and, particularly, powder bed fusion technology, can guarantee a final specimen quality which is comparable to standard manufacturing techniques. Actually, when metallic materials are considered, post-processing of specimens is still needed to remove residual stresses caused by temperature gradients during the printing phase. Generally speaking, printed materials deposited in layers may exhibit anisotropic behavior, due to their lamellar structure, which needs to be evaluated. The interface between the layers, depending on the process, remains a sensitive area, which may represent the mechanically weak points of the structure. It is therefore necessary to establish a study methodology to characterize the complex and anisotropic behavior of this type of material to further improve existing additive manufacturing processes. However, even simple microstructures realized by means of 3D printing have been proven to result in material properties which are extremely interesting from the engineering viewpoint, like extension-bending/exten- sion-torsion coupling [28–31], band gap properties for vibration attenuation [32,33], extreme stiffness [34], fluid-like behavior [35], extreme elastic range [36], damage-tolerant behavior [37], enhanced and/or tailored stability properties [38,39], unidirectional response [40,41], and many others [42]. It is evident that designing microstructures to be 3D printed, possibly by topological optimization techniques [43], requires a thorough quantitative understanding of the mechanical behavior of the struc- ture to be printed. Continuum models are needed to describe the large-scale mechanical behavior of a wide class of materials that, at some spatial scale, possess a microstructure. These models are particu- larly important for bridging spatial scales ranging from that of interactions between elements of the microstructure to the collective behavior of a large number of elements [44]. In the last 10 years, making use of generalized and second gradient continua theories introduced in the early 1960s [45–47], many investigations were carried out dealing with pantographic structures. In particular, the paper [26] pro- vided the first homogenization result for Bi-PS in large deformation regime along with numerical results and experimental measurements for validation purposes. Discrete simulations can be designed that provide trajectories and distribution of strain energies at the element scale. Therefore, discrete models can be used as a more refined description than continuum models. While continuum models are less computationally demanding than discrete models, the formu- lation of continuum models able of capturing the behavior of a material with microstructure, as well as the development of robust numerical procedures [48,49] that work in sometimes-singular cases, is still Figure 1. Quadratic Bi-PS (sample A1) with a CAD detail view of the inner structure on the left-hand side, developed by Barchiesi and colleagues [26,27], printed on a Formagia P 100 (EOS GmbH, Krailling, Germany) using SLS. A rectangle Bi-PS (sample B1) is shown on the right-hand side. The detail view shows a sample with different boundary conditions in order to achieve reinforcement. Harsch et al. 2203 challenging and, in any case, requires the utilization of discrete simulations for validation purposes. Therefore, when exploring the behavior of a material with microstructure, while less appealing in view of getting analytical expressions, the utilization of discrete simulations is of utmost utility. For the under- standing of Bi-PS, the use of discrete simulations has been indeed beneficial. Discrete modeling has been used as microscale description for developing a homogenized continuum description of Bi-PS and, albeit not systematically as in the present work, its experimental fitting has been carried out, together with that of the homogenized continuum model, in the paper [27]. In this paper, motivated by the encouraging results obtained in Barchiesi et al. [27], making use of training and test sets of experimental data, a numerical identification [50–52] is performed to fit the parameters of a planar discrete spring model describing the static behavior of Bi-PS subjected to dif- ferent kind of loads and boundary conditions. The model is fitted using force–displacement curves for bias extension, shearing, and point-load tests in large deformation regime for rectangular and square specimens. The attempt reported in this paper of performing a systematic fitting and valida- tion procedure for Bi-PS is inspired from existing literature on the identification of models describing pantographic fabrics. In De Angelo et al. [53], the in-plane and out-of-plane stiffnesses of a conti- nuum plate model are determined by simulating bias extension and shear tests and employing the Levenberg–Marquardt optimization algorithm. In Shekarchizadeh et al. [52], the numerical identifi- cation of the constitutive parameters of the aforementioned macromodel is carried out through the trust region reflective optimization method. Contrarily to this work, which makes use of experimental data, these two works made use of data generated by a refined 3D model based on Cauchy mechanics and hyperelastic constitutive laws. Developing an automatized optimization process for identifying the constitutive parameters of dis- crete modeling of Bi-PS and testing its descriptiveness and predictive capacity are the main objectives of the work, which is organized as follows. Section 2 introduces the manufacturing technology employed to realize the specimens, the experimental setup, and the discrete modeling of Bi-PS via elastic springs. Section 3 illustrates the minimization problem to be solved by the gradient free Nelder–Mead simplex algorithm for determining an optimal set of model parameters on the basis of experimental measure- ments, and presents the results of the analysis. Finally, conclusions and outlooks are briefly sketched. Building on top of previous work, the analysis allows to conclude that there exists a single set of para- meters for the adopted discrete modeling, which was also employed in previous papers as microscale model for homogenization purposes, such that it is descriptive of several different tests for a wide range of deformations and such that the model is predictive when experimental deformed shapes are compared with numerically computed ones. 2. Materials and methods Two batches of specimens, listed in Table 1, consisting of three quadratic Bi-PS (sample family A, see left of Figure 1) and two rectangle Bi-PS (sample family B, see right of Figure 1), have been printed using a selective laser sintering (SLS) technique. Each sample family was tested in an extension test and shearing test (see Table 1). In addition, a so-called point-force test was applied to the quadratic sample Table 1. Overview of both sample families (A1–A3, B1–B2) including their loading conditions (LCs) (extension, shear, and point-force) and their boundary conditions (BCs) (reinforced or not reinforced). N is the number of unit-cells in vertical direction, M is the number of unit-cells in horizontal direction, and the parameter e describes the diagonal width of a half unit-cell (corresponding to Figure 6). Sample LC BC N M e(mm) A1 (square) Extension Not reinforced 9 9 17 � ffiffiffi 2 p B1 (rectangle) Extension Reinforced 7 11 17 � ffiffiffi 2 p A2 (square) Shear Reinforced 9 9 17 � ffiffiffi 2 p B2 (rectangle) Shear Reinforced 7 11 17 � ffiffiffi 2 p A3 (square) Point-force Not reinforced 9 9 17 � ffiffiffi 2 p 2204 Mathematics and Mechanics of Solids 27(10) (A3). Furthermore, a discrete simulation based on the Hencky-type model is going to be introduced. Its numerical results will be compared with experimental measurement data. All samples have been manufactured by means of SLS. Polyamide (PA) powder (PA 2200, EOS GmbH, Krailling, Germany) with an average grain size of 0.056mm was fused in the commercial 3D printer EOS Formagia P 100 (EOS GmbH) located at the Institute of Mechanics and Printing, University of Technology Warsaw, Poland. In this bottom-up procedure, layers are sintered in the heated, focused laser spot. The specimen is manufactured layer-by-layer by means of a laser beam under inert gas atmosphere. After printing predefined points in the first layer, the table is lowered and shaken before new powder is applied by a wiper. The process repeats until the sample is finished. In order to reduce initial stresses, the sample has to cool down for several hours. All samples have been finished by means of a high pressure cleaner. Printing results are shown in Figure 1. A square Bi-PS (sample A1) with a CAD detail view of the inner structure on the left-hand side and a rectangular Bi-PS (sample B1) on the right-hand side with a detail view of reinforcements can be recognized. In order to overcome the limits of lack of contrast and therefore to enable DIC for the determination of displacements during real-time deformation, both sam- ple families have been speckled with black ink on the outer surface. CAD drawings with their inner and outer geometric parameters are shown in Figure 2. 2.1. Experimental setup All experiments have been performed on an MTS Tytron 250 testing device (MTS Systems Corporation, Eden Prairie, MN, USA) at the Institute of Mechanics, Technical University Berlin, Germany. Three different experiments have been performed: standard extension tests, shearing tests, and a so-called point-force test (see Table 1). Furthermore, in order to investigate the effect of reinforcements at the bor- der of the inner substructure on the outer deformation behavior in more detail, different boundary con- ditions have been taken into account (see Table 1). By adapting a rigid body to the mounting side, a reinforcement was achieved (see detail view on the right-hand side of Figure 1). The schematics of the three different experimental setups are summarized in Figure 3. As shown in Figure 4 for the case of extension load, maximal reaction forces were measured by a device-own loading cell with �f = 6250N and a precision of about 625mN. The displacement Dx was imposed horizontally with a loading rate of _x = 1525mmmin�1. In addition, displacement was mea- sured by means of a non-invasive optical deformation technique—DIC. A commercial Canon EOS 1000D camera with a resolution of about 4272× 2848 pixels recorded one picture every 2 s. Two-dimen- sional (2D)-DIC evaluation was performed by means of GOM Correlate 2017 software (GOM GmbH, Figure 2. CAD drawings of a quadratic Bi-PS (specimen family A) and rectangle Bi-PS (specimen family B). Harsch et al. 2205 Braunschweig, Germany) and are in very good agreement with the machine-code measured by the device-own encoder. 2.2. Discrete model For solving the discrete micromodel, it is convenient to introduce a global, minimal set of generalized coordinates, such that the kinematics of the discrete system is entirely described by the coordinates of the nodal points. The bi-pantographic sheet is modeled as a discrete elastic spring system that is embedded in the 2D Euclidean vector space E 2. The static equilibrium conditions are obtained by the principle of virtual work, which in the case of our elastic system coincides with the principle of stationary Figure 3. Schematics of displacement-controlled (Dx) shear, extension, and point-force tests: (a) extension, (b) shear, and (c) point-force. Figure 4. Extension test of a not reinforced quadratic Bi-PS (A1). Displacement Dx is imposed horizontally on the right side while the reaction force �f is measured by the load cell on the fixed left side. 2206 Mathematics and Mechanics of Solids 27(10) total potential energy. Consequently, only the system’s potential energy and its corresponding variation have to be formulated. The system is composed of extensional and rotational springs only. First, we introduce the potential energies of a standard extensional and rotational spring element. Subsequently, we explain the kinematics of the bi-pantographic sheet, which includes some cumbersome but necessary bookkeeping of the relevant degrees of freedom for each spring contribution. Finally, we state the princi- ple of virtual work for the constrained discrete system subjected to kinematic boundary conditions. The springs are formulated between nodal points, which are depicted as white filled circles (see Figure 5). The position r= xex + yey 2 E 2 of a typical nodal point is commonly addressed by its Cartesian coor- dinates x= (x, y) 2 R 2 with respect to the orthonormal basis vectors ex, ey 2 E 2. If not stated otherwise, R f�tuples are considered in the sense of matrix multiplication as Rf × 1�matrices, i.e., as ‘‘column vec- tors.’’ Let qe = (x1, y1, x2, y2) 2 R 4 be the coordinates of two points interconnected by an extensional spring as depicted in Figure 5(a). Introducing the abbreviations Dx = x2 � x1 and Dy = y2 � y1, the dis- tance between the two points is as follows: l(qe) = ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Dx2 + Dy2 p = ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi (x2 � x1)2 + (y2 � y1)2 q : ð1Þ The derivative with respect to qe is the ‘‘row vector’’: ∂l ∂qe (qe) = 1 l(qe) (�Dx, �Dy,Dx,Dy) 2 R 1×4: ð2Þ Without loss of generality, we restrict ourselves to linear springs (cf. Andreaus et al. [54], equation (4)), which are modeled by means of a quadratic potential with stiffness ke . 0 and undeformed length l0 . 0. That is the function: Ee(ke, l0, qe) = 1 2 ke l(qe)� l0½ �2: ð3Þ The derivative of the potential (3) with respect to qe is the function: ∂Ee ∂qe (ke, l0, q e) = ke l(qe)� l0½ � ∂l ∂qe (qe) =: fe(ke, l0, q e)ð ÞT 2 R 1×4, ð4Þ where the transposed of a matrix is indicated by (�)T. Rotational springs are interactions between three nodal points. Let qr = (x1, y1, x2, y2, x3, y3) 2 R 6 be the Cartesian coordinates of three nodal points as depicted in Figure 5(b). With the abbreviations Dx1 = x2 � x1, Dx2 = x3 � x2, Dy1 = y2 � y1, and Dy2 = y3 � y2, the angles between the ex�axis and the vectors Dr1 = Dx1ex + Dy1ey and Dr2 = Dx2ex + Dy2ey, respectively, are introduced by the relations: Figure 5. (a) Kinematics of an extensional spring and (b) kinematics of a rotational spring. Harsch et al. 2207 f1(qr) = tan�1 Dy1 Dx1 � � , f2(qr) = tan�1 Dy2 Dx2 � � , ð5Þ with the corresponding derivatives: ∂f1 ∂qr (qr) = 1 l1(qr)2 (Dy1, �Dx1, �Dy1,Dx1, 0, 0) 2 R 1×6, ∂f2 ∂qr (qr) = 1 l2(qr)2 (0, 0,Dy2, �Dx2, �Dy2,Dx2) 2 R 1×6: ð6Þ The relative angle Df = f2 � f1 measures the torsion of a rotational spring. Without loss of general- ity, we restrict ourselves to linear springs (cf. dell’Isola et al. [55], equation (2.1)), which are modeled by means of a quadratic potential with stiffness kr . 0 and undeformed relative angle Df0. Their energy is given by the real-valued function: Er(kr,Df0, q r) = 1 2 kr Df(qr)�Df0½ �2: ð7Þ Straightforwardly, the derivative with respect to qr leads to the R1×6 matrix: ∂Er ∂qr (kr,Df0, q r) = kr Df(qr)�Df0½ � ∂f2 ∂qr (qr)� ∂f1 ∂qr (qr) � � =: fr(kr,Df0, q r)ð ÞT: ð8Þ Since tan21 is restricted to the interval ð�p 2 , p 2 Þ, it is not always possible to compute the absolute angles f1 and f2 from equation (5). This can be overcome by numerically computing tan21 via the atan2-function with a range (�p , + p) � R (cf. Barchiesi et al. [27]). Even with this correction, the value of the generalized force law of the torsional spring (8) will be discontinuous at atan2(�1, 0) = p. This singularity can be removed by direct computation of the relative angle Df. For that, a local basis that is aligned with Dr1 is constructed via: eK x = Dr1=kDr1k, eK y = Dr?1 =kDr?1 k, ð9Þ with Dr?1 = (ey � ex � ex � ey)Dr1 being Dr1 rotated counterclockwise by a quarter angle. Introducing the decomposition: Dr2 = Dr k 2e K x + Dr?2 e K y , Dr k 2 = Dr2 � Dr1 kDr1k , Dr?2 = Dr2 � Dr?1 kDr?1 k , ð10Þ and recognizing kDr1k = kDr?1 k , the relative angle can be computed via: Df = tan�1 Dr?2 Dr k 2 ! = tan�1 Dr2 � Dr?1 Dr2 � Dr1 � � = tan�1 Dx2Dy1�Dx1Dy2 Dx2Dx1 + Dy2Dy1 � � : ð11Þ Again, numerically tan21 is evaluated via the atan2-function. By that, no singularities are present for jDfj\p, i.e., every rotational spring’s relative angle is smaller than 180�. Indeed, this is a valid assumption for the numerical simulations of the bi-pantographic sheet shown in the subsequent section. The bi-pantographic sheet is characterized by a periodic substructure, which is captured in the dis- crete model by identical cells composed of extensional and rotational springs. The cells themselves do interact with each other by sharing nodal points with adjacent cells and by additional rotational springs. In the undeformed reference configuration, the cells are quadratic with width ffiffiffi 2 p e (see Figure 6). Each cell has four nodal points as vertices, which are globally addressed by the row-column indices (i, j) within the range i = f0, 1, . . . ,Ng and j = f0, 1, . . . ,Mg. Consequently, the sheet consists of N �M cells with (N + 1)(M + 1) cell vertices. Besides the four vertices, in each cell, there are nine additional nodal points required to capture the complex behavior of the bi-pantographic sheet. In the reference config- uration, these nodal points are arranged symmetrically as depicted in Figure 6. 2208 Mathematics and Mechanics of Solids 27(10) In the deformed configuration, the Cartesian coordinates of the cell vertices are denoted x(i, j) = (x(i, j), y(i, j)) 2 R 2. As depicted in Figure 7, the Cartesian coordinates of the kth nodal point within cell (i, j) are xk (i, j) = (xk (i, j), yk (i, j)) 2 R 2. Note that the coordinates of the cell vertices have different denota- tions depending on their respective application on cell or sheet level. Thus, the f = 2½(N + 1)(M + 1) + 9N �M � generalized coordinates of the discrete system are the Cartesian coordi- nates of all nodal points of the sheet: q= (x(0, 0), . . . , x(M ,N), x 4 (0, 0), . . . , x12 (0, 0), . . . , x4 (M�1,N�1), . . . , x12 (M�1,N�1)) 2 R f : ð12Þ For a compact formulation of the total potential energy of the system, it is convenient to introduce spring coordinates, i.e., sets of coordinates that involve only the coordinates of the relevant nodal points for each spring. We begin with the interactions in each cell (i, j) by considering Figure 7 together with Table 2, where the explicit assignments of the nodal coordinates to the spring coordinates are specified. The axial stiffnesses of the fibers in the sheet are modeled by 16 extensional springs each with stiffness ke and undeformed length l0 = e=(2 cosb) = e= ffiffiffi 3 p , where we have used b = p 4 � a = p 6 . The coordinates required for the kth spring are given by qe(i, j)k 2 R 4. The resistance of the pins with respect to torsion is captured by eight rotational springs, depicted in Figure 7 by green arcs, with stiffness kr1. The energy of Figure 6. Reference configuration of the discrete model. The (i, j)th vertex and the (i, j)th cell are highlighted by a red dot and a red square, respectively. A close up shows the arrangement of the nodal points in the undeformed reference configuration. Figure 7. Deformed configuration of the discrete model. The (i, j)th vertex and the (i, j)th cell are highlighted by a red dot and a red quadrilateral, respectively. A close up shows the arrangement of the nodal points in the deformed configuration and their interaction by extensional and rotational springs. Harsch et al. 2209 the lth rotational spring is formulated with qr1(i, j)l 2 R 6 and the undeformed angle f0 = (�1)lp=3. The bending stiffness of the fibers within each cell is included by four rotational springs with stiffness kr2 and f0 = 0 influencing the coordinates qr2(i, j)m 2 R 6. The cells interact with each other by sharing nodal points with adjacent cells. However, we also have to account for the bending stiffness of the fibers at each cell vertex for i = 1, . . . ,N � 1 and j = 1, . . . ,M � 1. These bending stiffnesses are realized by rota- tional springs with stiffness kr2 and corresponding spring coordinates qr3(i, j)m 2 R 6, which are specified in Figure 8. To extract the spring coordinates from the generalized coordinates (12), the Boolean connec- tivity matrices Ce (i, j)k 2 R 4× f and Cr1 (i, j)l,C r2 (i, j)m,C r3 (i, j)m 2 R 6× f are defined by the relations: qe(i, j)k =Ce (i, j)kq , qr1(i, j)l =Cr1 (i, j)lq , qr2(i, j)m =Cr2 (i, j)mq , qr3(i, j)m =Cr3 (i, j)mq : ð13Þ Using the spring coordinates of Table 2 and Figure 8 together with the potential energies (3) and (7), the total potential energy of the discrete bi-pantographic sheet is as follows: E(q) = XN�1 i = 0 XM�1 j = 0 X15 k = 0 Ee(ke, l0,C e (i, j)kq) + X7 l = 0 Er kr1, (�1)lp 3 ,Cr1 (i, j)lq � � + X3 m = 0 Er kr2, 0,C r2 (i, j)mq � �! + XN�1 i = 1 XM�1 j = 1 X3 m = 0 Er kr2, 0,Cr3 (i, j)mq � � : ð14Þ Let q̂(e) be a function of e that includes the actual coordinates q for static equilibrium in the case of e= 0, i.e., q̂(0) = q. Then, the variation of the total potential energy E induced by q̂ is as follows: dE(q) = ∂E ∂q (q) dq̂ de (0) = fT(q)dq , ð15Þ where dq= dq̂ de (0) are the virtual displacements, and f(q) = (∂e=∂q)T(q) are the internal generalized forces. Using equations (4) and (8) together with the total potential energy (14), the internal generalized forces of the bi-pantographic sheet are obtained by: fT(q) = ∂E ∂q (q) = = XN�1 i = 0 XM�1 j = 0 X15 k = 0 fe(ke, l0,C e (i, j)kq)Ce (i, j), k + X7 l = 0 fr(kr1, (�1)lp 3 ,Cr1 (i, j)lq)Cr1 (i, j)l + X3 m = 0 fr(kr2, 0,Cr2 (i, j)mq)Cr2 (i, j)m ! + XN�1 i = 0 XM�1 j = 0 X3 m = 0 fr(kr2, 0,C r3 (i, j)mq)Cr3 (i, j)m: ð16Þ Table 2. Assignments of coordinates of nodal points to spring coordinates within cell (i, j). Spring coordinates Assignment of nodal points qe (i, j)k = (x p (i, j), x q (i, j)) fp, qg 2 f0, 5g k = 0 , f0, 6g k = 1 , f5, 4g k = 2 , f6, 4g k = 3 , f4, 7g k = 4 , f4, 8g k = 5 , f7, 1g k = 6 , f8, 1g k = 7 , f4, 9g k = 8 , f4, 10g k = 9 , f9, 2g k = 10 , f10, 2g k = 11 , f4, 11g k = 12 , f4, 12g k = 13 , f11, 3g k = 14 , f12, 3g k = 15 qr1 (i, j)l = (x p (i, j), x q (i, j), xs (i, j)) fp, q, sg 2 f0, 6, 4g l = 0 , f0, 5, 4g l = 1 , f4, 7, 1g l = 2 , f4, 8, 1g l = 3 , f4, 9, 2g l = 4 , f4, 10, 2g l = 5 , f3, 12, 4g l = 6 , f3, 11, 4g l = 7 qr2 (i, j)m = (x p (i, j), x q (i, j), xs (i, j)) fp, q, sg 2 f5, 4, 9g m = 0 , f6, 4, 10g m = 1 , f12, 4, 8g m = 2 , f11, 4, 7g m = 3 2210 Mathematics and Mechanics of Solids 27(10) Kinematic boundary conditions are imposed by perfect bilateral constraints 0= g(q) 2 R m with a vir- tual work contribution dW c = dgT l = dqTW(q) l , whereW(q)T = ∂g ∂q (q) 2 R m× f is the matrix of gener- alized force directions and l 2 R m is the vector of constraint forces. In such a manner, the position of the most outer nodes is constrained to zero or a prescribed displacement, depending on the respective left and right boundaries. In addition, the central node x0 i, j is constrained on these boundaries in order to model the reinforcement using an additional pivot as shown in the first, second, and fourth experiments in Figure 10. This assumption is fostered by observing that the stiffness of reinforcement is of magni- tudes larger compared to the extensional stiffness of the axial springs. For the point-load experiment, all nodes on the left-hand side are connected to a 2D rigid body. This body is constrained such that it is only allowed to rotate around a pivoted point with relative horizontal (25.1mm) and vertical (24.35mm) displacement with respect to the lower left node of the specimen. The discrete system is now in static equilibrium iff the total virtual work of internal generalized forces and constraint forces vanishes for all virtual displacements, i.e., dqT(f(q) +W(q)l) = 0 8dq 2 R f , ð17Þ and the constraint equations are g(q) = 0 are satisfied. Thus, the equilibrium configuration is determined by the set of nonlinear equations: h(q,l) = f(q) +W(q)l g(q) � = 0, ð18Þ which can be solved, at least locally, by a Newton–Raphson iteration scheme. 3. Experiments This section presents the results of both, the simulation and the real world experiments. For identifica- tion of the unknown stiffness parameters, a parameter optimization procedure will be presented in the beginning. Subsequently, it is shown that the numerical model is capable for description of the panto- graphic sheet using a single set of experimentally optimized parameters. 3.1. Parameter optimization The yet unknown stiffness parameters b = (ke, kr1, kr2) are identified experimentally in the following way and reported in Table 3. Let the Euclidean error between the simulated force fij and the measured one f ij (see Figure 11) normalized by the maximum force f ijmax within each experiment be denoted as: dij = fij � f ij f ijmax : ð19Þ Figure 8. Assignments of coordinates of nodal points to spring coordinates for interaction between cells. Harsch et al. 2211 This error was minimized during an optimization procedure using all N = 5 experiments. Within each experiment M = 20, different static equilibria were solved. The optimal values b� were obtained by introducing the cost function: K(b) = 1 2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiXN i = 1 XM j = 1 d2 ij N vuut , ð20Þ and solving the optimization problem: b�= argmin b2R3 K(b) , s:t: h(q,l) = 0, ð21Þ using a gradient free Nelder–Mead simplex algorithm (see Gao and Han [56]). By heuristic, the initial parameters were set to ke = 20Nm�1, kr1 = 10Nm, and kr2 = 50Nm. In order to study the sensitivity of the cost function a grid search was performed. For that the value of the cost function is computed for the optimal parameters b� varied by 630%. On the left-hand side of Figure 9, the loss surface of the cost function is depicted for the optimal parameter ke = 15:03Nm�1 but varied kr1 and kr2. It can clearly be seen that for the optimal value of ke, a variation of kr1 and kr2 increases the global error. Led by that observation the cost function is evaluated for the optimal values kr1 = 14:53Nm and kr2 = 22:10Nm but varied values of ke (see right-hand side of Figure 9). Again, for the optimal values, every variation in ke leads to an increasing cost function, which shows that the found set of optimal parameters is a local minimizer for the given cost function (20). 3.2. Experiment versus simulation In Figure 10, the initial and deformed configurations are shown for all five experiments. It can be seen that the numerically obtained deformed configurations are in excellent agreement with the experiments, neglecting optical distortion effects introduced by the camera system. There is only a single outlier, the Table 3. Optimized stiffness parameters and respective relative error. ke kr1 kr2 K(b�) 15.03 Nm21 14.53 Nm 22.10 Nm 1:497× 10�1 N Figure 9. Loss surface for optimal ke and varied kr1, kr2 (left). Loss function for fixed kr1 and kr2 but varied ke (right). 2212 Mathematics and Mechanics of Solids 27(10) Figure 10. Initial and deformed configurations—simulated versus experimental results. First and second rows show, respectively, the extension and shear deformation of the rectangular specimen. Third, fourth, and sixth rows show, respectively, the extension, shear, and point deformation of the quadratic specimen. The simulated configurations are visualized by overlayed blue lines. All simulations are performed with the same parameter set given in Table 3. Harsch et al. 2213 Figure 11. Force–displacement curves of the performed experiments. (Top left) Extension rectangular specimen (B1), (top right) shear rectangular specimen (B2), (center left) extension quadratic specimen (A1), (center right) shear quadratic specimen (B2), and (bottom left) point test quadratic specimen (B3). All simulations are performed with the same parameter set given in Table 3. 2214 Mathematics and Mechanics of Solids 27(10) rectangle specimen during the extension experiment. It can be observed that during the numerical experi- ment, some elements of the pantographic sheet intersect each other. Obviously, this is no physical beha- vior and stems from the fact that self-contact and friction within the pantographic sheet is not included to the numerical model. By that, the model underestimates the required force and yields wrong displace- ment predictions since self-contact would prevent the torsional springs from deforming further. The force–displacement curves of the very same experiments are reported in Figure 11. Once again, for all but the first experiment, the simulated results are in excellent agreement with the experiments. This can again be explained by the missing self-contact, since the specimen is allowed to penetrating itself. If this would be prohibited by an appropriate contact formulation, the load on the extensional springs would increase and lead to higher internal forces as observed in the experiments. Furthermore, in the last two experiments, the forces are overestimated by the numerical model for very large displacements. 4. Discussion and outlook Within this publication, we presented a set of five displacement-controlled experiments of pantographic sheets. They demonstrated three different kinds of deformation states, namely, extension, shear, and point loads. All applied displacements are in the regime of medium to large elastic deformations. Furthermore, a new formulation for the computation of the intermediate angle required for the tor- sional springs of the discrete model removed the presence of geometric singularities of the existing model [27]. The results presented in the previous section have shown that the presented discrete model is a both powerful and easy tool for description of large deformations of pantographic structures. Good predictions of the deformations and the required forces are obtained. Future investigations should clar- ify if the usage of nonlinear material models for the extensional and torsional springs lead to even more accurate predictions. In order to include all relevant physical effects, frictional self-contact [57,58] has to be included to the discrete model. By such an improved numerical model, the presented investigations can be extended the area of extremely large elastic deformations. Furthermore, dynamic analysis and nonlinear wave propagation of pantographic fabrics [59] can be investigated using the present model. From the experimental point of view, optical distortion effects introduced by the camera system have to be reduced or avoided by an appropriate post-processing step. Finally, the individual nodal displace- ments and further derived quantities, e.g., intermediate angles, can be used for a more sophisticated opti- mization procedure. This requires the application of DIC (cf. [27,60]) or a full-field magnetic resonance imaging (MRI)/computed tomography (CT) scan. Funding The author(s) received no financial support for the research, authorship, and/or publication of this article. ORCID iDs Jonas Harsch https://orcid.org/0000-0003-1065-2216 Simon R Eugster https://orcid.org/0000-0002-4562-1287 References 1. Placidi, L, Barchiesi, E, Turco, E, et al. A review on 2D models for the description of pantographic fabrics. Z Angew Math Phys 2016; 67(5): 1–20. 2. Seppecher, P, Alibert, JJ, Lekszycki, T, et al. Pantographic metamaterials: an example of mathematically driven design and of its technological challenges. Contin Mech Thermodyn 2019; 31(4): 851–884. 3. Cuomo, M, Boutin, C, Contrafatto, L, et al. Effective anisotropic properties of fibre network sheets. Eur J Mech-A/Solid 2022; 93: 104492. 4. Hild, F, and Misra, A. Multiscale DIC applied to pantographic structures. Exp Mech 2021; 61(2): 431–443. 5. Spagnuolo, M, Andreaus, U, Misra, A, et al. Mesoscale modeling and experimental analyses for pantographic cells: effect of hinge deformation. Mech Mater 2021; 160: 103924. 6. Giorgio, I. Numerical identification procedure between a micro-Cauchy model and a macro-second gradient model for planar pantographic structures. Z Angew Math Phys 2016; 67(4): 1–17. Harsch et al. 2215 https://orcid.org/0000-0003-1065-2216 https://orcid.org/0000-0002-4562-1287 7. Giorgio, I, Rizzi, N, and Turco, E. Continuum modelling of pantographic sheets for out-of-plane bifurcation and vibrational analysis. Proc R Soc A: Math Phys Eng Sci 2017; 473(2207): 20170636. 8. Giorgio, I, Rizzi, NL, Andreaus, U, et al. A two-dimensional continuum model of pantographic sheets moving in a 3D space and accounting for the offset and relative rotations of the fibers. Math Mech Compl Syst 2019; 7(4): 311–325. 9. Spagnuolo, M, Yildizdag, ME, Pinelli, X, et al. Out-of-plane deformation reduction via inelastic hinges in fibrous metamaterials and simplified damage approach. Math Mech Solid 2022; 27(6): 1011–1031. 10. Yang, H, Ganzosch, G, Giorgio, I, et al. Material characterization and computations of a polymeric metamaterial with a pantographic substructure. Z Angew Math Phys 2018; 69(4): 1–16. 11. Spagnuolo, M, Peyre, P, and Dupuy, C. Phenomenological aspects of quasi-perfect pivots in metallic pantographic structures. Mech Res Commun 2019; 101: 103415. 12. Golaszewski, M, Grygoruk, R, Giorgio, I, et al. Metamaterials with relative displacements in their microstructure: technological challenges in 3D printing, experiments and numerical predictions. Contin Mech Thermodyn 2019; 31(4): 1015–1034. 13. Turco, E, Golaszewski, M, Giorgio, I, et al. Pantographic lattices with non-orthogonal fibres: experiments and their numerical simulations. Compos Part B: Eng 2017; 118: 1–14. 14. Scerrato, D, Giorgio, I, and Rizzi, NL. Three-dimensional instabilities of pantographic sheets with parabolic lattices: numerical investigations. Z Angew Math Phys 2016; 67(3): 1–19. 15. Seppecher, P, Spagnuolo, M, Barchiesi, E, et al. Advances in pantographic structures: design, manufacturing, models, experiments and image analyses. Contin Mech Thermodyn 2019; 31(4): 1231–1282. 16. Juritza, A, Yang, H, and Ganzosch, G. Qualitative investigations of experiments performed on 3D-FDM-printed pantographic structures made out of PLA. In: Abali, B, Altenbach, H, dell’Isola, F, et al. (eds) New achievements in continuum mechanics and thermodynamics. Cham: Springer, 2019, pp. 197–209. 17. Trippel, A, Stilz, M, Gutmann, F, et al. Device for characterizing rotational joints in metamaterials. Mech Res Commun 2020; 104: 103501. 18. Nejadsadeghi, N, Laudato, M, De Angelo, M, et al. Mechanical behavior investigation of 3D printed pantographic unit cells via tension and compression tests. In: Abali, B, and Giorgio, I (eds) Developments and novel approaches in biomechanics and metamaterials. Cham: Springer, 2020, pp. 409–422. 19. dell’Isola, F, Turco, E, Misra, A, et al. Force–displacement relationship in micro-metric pantographs: experiments and numerical simulations. Compt Rend Mécan 2019; 347(5): 397–405. 20. Vangelatos, Z, Yildizdag, ME, Giorgio, I, et al. Investigating the mechanical response of microscale pantographic structures fabricated by multiphoton lithography. Ext Mech Lett 2021; 43: 101202. 21. Barchiesi, E, Eugster, SR, Placidi, L, et al. Pantographic beam: a complete second gradient 1D-continuum in plane. Z Angew Math Phys 2019; 70(5): 1–24. 22. Turco, E, and Barchiesi, E. Equilibrium paths of Hencky pantographic beams in a three-point bending problem. Math Mech Compl Syst 2019; 7(4): 287–310. 23. Barchiesi, E, Laudato, M, and Di Cosmo, F. Wave dispersion in non-linear pantographic beams. Mech Res Commun 2018; 94: 128–132. 24. Turco, E. Stepwise analysis of pantographic beams subjected to impulsive loads. Math Mech Solid 2021; 26(1): 62–79. 25. Seppecher, P, Alibert, JJ, and Isola, FD. Linear elastic trusses leading to continua with exotic mechanical interactions. J Phys: Conf Ser 319: 012018. 26. Barchiesi, E, Eugster, SR, dell’Isola, F, et al. Large in-plane elastic deformations of bi-pantographic fabrics: asymptotic homogenization and experimental validation. Math Mech Solid 2020; 25(3): 739–767. 27. Barchiesi, E, Harsch, J, Ganzosch, G, et al. Discrete versus homogenized continuum modeling in finite deformation bias extension test of bi-pantographic fabrics. Contin Mech Thermodyn. Epub ahead of print 12 September 2020. DOI: 10.1007/ s00161-020-00917-w. 28. Giorgio, I, dell’Isola, F, and Misra, A. Chirality in 2D Cosserat media related to stretch-micro-rotation coupling with links to granular micromechanics. Int J Solid Struct 2020; 202: 28–38. 29. Misra, A, Nejadsadeghi, N, De Angelo, M, et al. Chiral metamaterial predicted by granular micromechanics: verified with 1D example synthesized using additive manufacturing. Contin Mech Thermodyn 2020; 32(5): 1497–1513. 30. De Angelo, M, Placidi, L, Nejadsadeghi, N, et al. Non-standard Timoshenko beam model for chiral metamaterial: identification of stiffness parameters. Mech Res Commun 2020; 103: 103462. 31. Fernandez-Corbaton, I, Rockstuhl, C, Ziemke, P, et al. New twists of 3D chiral metamaterials. Adv Mater 2019; 31(26): 1807742. 32. Chen, Y, Huang, G, and Sun, C. Band gap control in an active elastic metamaterial with negative capacitance piezoelectric shunting. J Vib Acoust 2014; 136(6): 061008. 33. El Sherbiny, MG, and Placidi, L. Discrete and continuous aspects of some metamaterial elastic structures with band gaps. Arch Appl Mech 2018; 88(10): 1725–1742. 34. Wang, Y, and Sigmund, O. Quasiperiodic mechanical metamaterials with extreme isotropic stiffness. Ext Mech Lett 2020; 34: 100596. 2216 Mathematics and Mechanics of Solids 27(10) 35. Kadic, M, Bückmann, T, Schittny, R, et al. On anisotropic versions of three-dimensional pentamode metamaterials. N J Phys 2013; 15(2): 023029. 36. Mirzaali, M, Hedayati, R, Vena, P, et al. Rational design of soft mechanical metamaterials: independent tailoring of elastic properties with randomness. Appl Phys Lett 2017; 111(5): 051903. 37. Turco, E, dell’Isola, F, Rizzi, NL, et al. Fiber rupture in sheared planar pantographic sheets: numerical and experimental evidence. Mech Res Commun 2016; 76: 86–90. 38. Vangelatos, Z, Melissinaki, V, Farsari, M, et al. Intertwined microlattices greatly enhance the performance of mechanical metamaterials. Math Mech Solid 2019; 24(8): 2636–2648. 39. Vangelatos, Z, Gu, GX, and Grigoropoulos, CP. Architected metamaterials with tailored 3D buckling mechanisms at the microscale. Ext Mech Lett 2019; 33: 100580. 40. Barchiesi, E, dell’Isola, F, Bersani, AM, et al. Equilibria determination of elastic articulated duoskelion beams in 2D via a Riks-type algorithm. Int J Non-Lin Mech 2021; 128: 103628. 41. Turco, E, Barchiesi, E, and dell’Isola, F. In-plane dynamic buckling of duoskelion beam-like structures: discrete modeling and numerical results. Math Mech Solid 2021; 27: 1164–1184. 42. Barchiesi, E, Spagnuolo, M, and Placidi, L. Mechanical metamaterials: a state of the art. Math Mech Solid 2019; 24(1): 212–234. 43. Djourachkovitch, T, Blal, N, Hamila, N, et al. Multiscale topology optimization of 3D structures: a micro-architectured materials database assisted strategy. Comput Struct 2021; 255: 106574. 44. Abdoul-Anziz, H, and Seppecher, P. Strain gradient and generalized continua obtained by homogenizing frame lattices. Math Mech Compl Syst 2018; 6(3): 213–250. 45. Mindlin, RD. Microstructure in linear elasticity. Technical report, Department of Civil Engineering and Engineering Mechanics, Columbia University, New York, 1963. 46. dell’Isola, F, Corte, AD, and Giorgio, I. Higher-gradient continua: the legacy of Piola, Mindlin, Sedov and Toupin and some future research perspectives. Math Mech Solid 2017; 22(4): 852–872. 47. dell’Isola, F, Andreaus, U, and Placidi, L. At the origins and in the vanguard of peridynamics, non-local and higher- gradient continuum mechanics: an underestimated and still topical contribution of Gabrio Piola. Math Mech Solid 2015; 20(8): 887–928. 48. Abali, BE, Müller, WH, and Eremeyev, VA. Strain gradient elasticity with geometric nonlinearities and its computational evaluation. Mech Adv Mat Mod Process 2015; 1(1): 1–11. 49. Abali, BE, Müller, WH, et al. Theory and computation of higher gradient elasticity theories based on action principles. Arch Appl Mech 2017; 87(9): 1495–1510. 50. Barchiesi, E, Misra, A, Placidi, L, et al. Granular micromechanics-based identification of isotropic strain gradient parameters for elastic geometrically nonlinear deformations. ZAMM-J Appl Math Mech/Z Angew Math Mech 2021; 101(11): e202100059. 51. Placidi, L, Andreaus, U, and Giorgio, I. Identification of two-dimensional pantographic structure via a linear D4 orthotropic second gradient elastic model. J Eng Math 2017; 103(1): 1–21. 52. Shekarchizadeh, N, Abali, BE, Barchiesi, E, et al. Inverse analysis of metamaterials and parameter determination by means of an automatized optimization problem. ZAMM-J Appl Math Mech/Z Angew Math Mech 2021; 101(8): e202000277. 53. De Angelo, M, Barchiesi, E, Giorgio, I, et al. Numerical identification of constitutive parameters in reduced-order bi- dimensional models for pantographic structures: application to out-of-plane buckling. Arch Appl Mech 2019; 89(7): 1333–1358. 54. Andreaus, U, Spagnuolo, M, Lekszycki, T, et al. A Ritz approach for the static analysis of planar pantographic structures modeled with nonlinear Euler–Bernoulli beams. Contin Mech Thermodyn 2018; 30(5): 1103–1123. 55. dell’Isola, F, Giorgio, I, Pawlikowski, M, et al. Large deformations of planar extensible beams and pantographic lattices: heuristic homogenization, experimental and numerical examples of equilibrium. Proc R Soc A: Math Phys Eng Sci 2016; 472(2185): 20150790. 56. Gao, F, and Han, L. Implementing the Nelder–Mead simplex algorithm with adaptive parameters. Comput Optim Appl 2012; 51(1): 259–277. 57. Glocker, C. Set-valued force laws: Dynamics of non-smooth systems, lecture notes in applied mechanics, vol. 1. Cham: Springer, 2001. 58. Leine, RI, and Nijmeijer, H. Dynamics and bifurcations of non-smooth mechanical systems, lecture notes in applied and computational mechanics, vol. 18. Berlin; Heidelberg: Springer, 2004. 59. Eugster, SR. Numerical analysis of nonlinear wave propagation in a pantographic sheet. Math Mech Compl Syst 2021; 9(3): 293–310. 60. Barchiesi, E. On the validation of homogenized modeling for bi-pantographic metamaterials via digital image correlation. Int J Solid Struct 2021; 208: 49–62. Harsch et al. 2217