Experimental analysis,discrete modeling and parameter optimization of SLS-printed bi-pantographic structures

Abstract

Making use of experimental data for bias extension, shearing, and point-load tests in large deformation regime for rectangular and square bi-pantographic specimens, we perform a numerical identification to fit the a priori parameters of a planar 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.

1. Introduction

This paper deals with a generalization of pantographic fabrics [1,2]. Pantographic fabrics are actual fabrics which are made up of two orthogonal families of parallel beams. Beams of different families are connected 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 dimensions and physical properties, then the pantographic fabric is said to be balanced [3]. Pantographic fabrics 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 torsional 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 realization of cylindrical pivots is less challenging from the manufacturing point of view—there are, however, 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 elements. 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 conditions 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 manufacturing 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 understanding 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 mention 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 deformation 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].

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.

One of the main advantages of additive manufacturing lies in the possibility of easily producing complex 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 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 materials used in additive manufacturing are still polymeric or metallic. When polymeric and metallic materials 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/extension-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 structure 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 particularly 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] provided 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 formulation 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 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 understanding 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 different 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 validation 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 continuum 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 identification 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 discrete 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 measurements, 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 parameters 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 (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.

Table 1.

Overview of both sample families (A1–A3, B1–B2) including their loading conditions (LCs) (extension, shear, andpoint-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 $\epsilon$ describes the diagonal width of a half unit-cell (corresponding to Figure 6).

Sample	LC	BC	N	M	$\epsilon (\mathrm{mm})$
A1 (square)	Extension	Not reinforced	9	9	$17/\sqrt{2}$
B1 (rectangle)	Extension	Reinforced	7	11	$17/\sqrt{2}$
A2 (square)	Shear	Reinforced	9	9	$17/\sqrt{2}$
B2 (rectangle)	Shear	Reinforced	7	11	$17/\sqrt{2}$
A3 (square)	Point-force	Not reinforced	9	9	$17/\sqrt{2}$

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.056 mm 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 sample 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.

Figure 2.

CAD drawings of a quadratic Bi-PS (specimen family A) and rectangle Bi-PS (specimen family B).

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 border of the inner substructure on the outer deformation behavior in more detail, different boundary conditions 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.

Figure 3.

Schematics of displacement-controlled $(\mathrm{\Delta }x)$ shear, extension, and point-force tests: (a) extension, (b) shear, and(c) point-force.

As shown in Figure 4 for the case of extension load, maximal reaction forces were measured by a device-own loading cell with $\overline{f}=\pm 250\mathrm{N}$ and a precision of about ±25 mN. The displacement $\mathrm{\Delta }x$ was imposed horizontally with a loading rate of $\stackrel{·}{x}=1525\mathrm{mm}\mathrm{mi}{\mathrm{n}}^{-1}$ . In addition, displacement was measured by means of a non-invasive optical deformation technique—DIC. A commercial Canon EOS 1000D camera with a resolution of about $4272\times 2848$ pixels recorded one picture every 2 s. Two-dimensional (2D)-DIC evaluation was performed by means of GOM Correlate 2017 software (GOM GmbH, Braunschweig, Germany) and are in very good agreement with the machine-code measured by the device-own encoder.

Figure 4.

Extension test of a not reinforced quadratic Bi-PS (A1). Displacement $\mathrm{\Delta }x$ is imposed horizontally on the right side while the reaction force $\overline{f}$ is measured by the load cell on the fixed left side.

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 ${\mathbb{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 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 principle 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 $\mathbf{r}=x{\mathbf{e}}_x+y{\mathbf{e}}_y\in {\mathbb{E}}^2$ of a typical nodal point is commonly addressed by its Cartesian coordinates $\mathbf{x}=(x,y)\in {\mathbb{R}}^2$ with respect to the orthonormal basis vectors ${\mathbf{e}}_x,{\mathbf{e}}_y\in {\mathbb{E}}^2$ . If not stated otherwise, ${\mathbb{R}}^f-\mathrm{tuples}$ are considered in the sense of matrix multiplication as ${\mathbb{R}}^{f\times 1}-\mathrm{matrices}$ , i.e., as “column vectors.” Let ${\mathbf{q}}^{\mathrm{e}}=(x_1,y_1,x_2,y_2)\in {\mathbb{R}}^4$ be the coordinates of two points interconnected by an extensional spring as depicted in Figure 5(a). Introducing the abbreviations $\mathrm{\Delta }x=x_2-x_1$ and $\mathrm{\Delta }y=y_2-y_1$ , the distance between the two points is as follows:

$$l({\mathbf{q}}^{\mathrm{e}})=\sqrt{\mathrm{\Delta }{x}^2+\mathrm{\Delta }{y}^2}=\sqrt{{(x_2-x_1)}^2+{(y_2-y_1)}^2}.$$

(1)

The derivative with respect to ${\mathbf{q}}^{\mathrm{e}}$ is the “row vector”:

$$\frac{\partial l}{\partial {\mathbf{q}}^{\mathrm{e}}}({\mathbf{q}}^{\mathrm{e}})=\frac{1}{l({\mathbf{q}}^{\mathrm{e}})}(-\mathrm{\Delta }x,-\mathrm{\Delta }y,\mathrm{\Delta }x,\mathrm{\Delta }y)\in {\mathbb{R}}^{1\times 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 $k_{\mathrm{e}}>0$ and undeformed length $l_0>0$ . That is the function:

$${\mathcal{E}}^{\mathrm{e}}(k_{\mathrm{e}},l_0,{\mathbf{q}}^{\mathrm{e}})=\frac{1}{2}{k}_{\mathrm{e}}{[l({\mathbf{q}}^{\mathrm{e}})-l_0]}^2.$$

(3)

The derivative of the potential (3) with respect to ${\mathbf{q}}^{\mathrm{e}}$ is the function:

$$\frac{\partial {\mathcal{E}}^{\mathrm{e}}}{\partial {\mathbf{q}}^{\mathrm{e}}}(k_{\mathrm{e}},l_0,{\mathbf{q}}^{\mathrm{e}})=k_{\mathrm{e}}[l({\mathbf{q}}^{\mathrm{e}})-l_0]\frac{\partial l}{\partial {\mathbf{q}}^{\mathrm{e}}}({\mathbf{q}}^{\mathrm{e}})=:{({\mathbf{f}}^{\mathrm{e}}(k_{\mathrm{e}},l_0,{\mathbf{q}}^{\mathrm{e}}))}^{\mathrm{T}}\in {\mathbb{R}}^{1\times 4},$$

(4)

where the transposed of a matrix is indicated by $(·{)}^{\mathrm{T}}$ .

Figure 5.

(a) Kinematics of an extensional spring and (b) kinematics of a rotational spring.

Rotational springs are interactions between three nodal points. Let ${\mathbf{q}}^{\mathrm{r}}=(x_1,y_1,x_2,y_2,x_3,y_3)\in {\mathbb{R}}^6$ be the Cartesian coordinates of three nodal points as depicted in Figure 5(b). With the abbreviations $\mathrm{\Delta }{x}_1=x_2-x_1$ , $\mathrm{\Delta }{x}_2=x_3-x_2$ , $\mathrm{\Delta }{y}_1=y_2-y_1$ , and $\mathrm{\Delta }{y}_2=y_3-y_2$ , the angles between the ${\mathbf{e}}_x-\mathrm{axis}$ and the vectors $\mathrm{\Delta }{\mathbf{r}}_1=\mathrm{\Delta }{x}_1{\mathbf{e}}_x+\mathrm{\Delta }{y}_1{\mathbf{e}}_y$ and $\mathrm{\Delta }{\mathbf{r}}_2=\mathrm{\Delta }{x}_2{\mathbf{e}}_x+\mathrm{\Delta }{y}_2{\mathbf{e}}_y$ , respectively, are introduced by the relations:

$${\varphi }_1({\mathbf{q}}^{\mathrm{r}})={\tan }^{-1}\left(\frac{\mathrm{\Delta }{y}_1}{\mathrm{\Delta }{x}_1}\right),{\varphi }_2({\mathbf{q}}^{\mathrm{r}})={\tan }^{-1}\left(\frac{\mathrm{\Delta }{y}_2}{\mathrm{\Delta }{x}_2}\right),$$

(5)

with the corresponding derivatives:

$$\begin{array}{cc}\multicolumn{1}{c}{\frac{\partial {\varphi }_1}{\partial {\mathbf{q}}^{\mathrm{r}}}({\mathbf{q}}^{\mathrm{r}})}& =\frac{1}{l_1{({\mathbf{q}}^{\mathrm{r}})}^2}(\mathrm{\Delta }{y}_1,-\mathrm{\Delta }{x}_1,-\mathrm{\Delta }{y}_1,\mathrm{\Delta }{x}_1,0,0)\in {\mathbb{R}}^{1\times 6},\hfill \\ \multicolumn{1}{c}{\frac{\partial {\varphi }_2}{\partial {\mathbf{q}}^{\mathrm{r}}}(q^{\mathrm{r}})}& =\frac{1}{l_2{({\mathbf{q}}^{\mathrm{r}})}^2}(0,0,\mathrm{\Delta }{y}_2,-\mathrm{\Delta }{x}_2,-\mathrm{\Delta }{y}_2,\mathrm{\Delta }{x}_2)\in {\mathbb{R}}^{1\times 6}.\hfill \end{array}$$

(6)

The relative angle $\mathrm{\Delta }\varphi ={\varphi }_2-{\varphi }_1$ measures the torsion of a rotational spring. Without loss of generality, 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 $k_{\mathrm{r}}>0$ and undeformed relative angle $\mathrm{\Delta }{\varphi }_0$ . Their energy is given by the real-valued function:

$${\mathcal{E}}^{\mathrm{r}}(k_{\mathrm{r}},\mathrm{\Delta }{\varphi }_0,{\mathbf{q}}^{\mathrm{r}})=\frac{1}{2}{k}_{\mathrm{r}}{[\mathrm{\Delta }\varphi ({\mathbf{q}}^{\mathrm{r}})-\mathrm{\Delta }{\varphi }_0]}^2.$$

(7)

Straightforwardly, the derivative with respect to ${\mathbf{q}}^{\mathrm{r}}$ leads to the ${\mathbb{R}}^{1\times 6}$ matrix:

$$\frac{\partial {\mathcal{E}}^{\mathrm{r}}}{\partial {\mathbf{q}}^{\mathrm{r}}}(k_{\mathrm{r}},\mathrm{\Delta }{\varphi }_0,{\mathbf{q}}^{\mathrm{r}})=k_{\mathrm{r}}[\mathrm{\Delta }\varphi ({\mathbf{q}}^{\mathrm{r}})-\mathrm{\Delta }{\varphi }_0](\frac{\partial {\varphi }_2}{\partial {\mathbf{q}}^{\mathrm{r}}}({\mathbf{q}}^{\mathrm{r}})-\frac{\partial {\varphi }_1}{\partial {\mathbf{q}}^{\mathrm{r}}}({\mathbf{q}}^{\mathrm{r}}))=:{({\mathbf{f}}^{\mathrm{r}}(k_{\mathrm{r}},\mathrm{\Delta }{\varphi }_0,{\mathbf{q}}^{\mathrm{r}}))}^{\mathrm{T}}.$$

(8)

Since tan⁻¹ is restricted to the interval $(-\frac{\pi }{2},\frac{\pi }{2})$ , it is not always possible to compute the absolute angles ${\varphi }_1$ and ${\varphi }_2$ from equation (5). This can be overcome by numerically computing tan⁻¹ via the atan2-function with a range $(-\pi ,+\pi )\subset \mathbb{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)=\pi$ . This singularity can be removed by direct computation of the relative angle $\mathrm{\Delta }\varphi$ . For that, a local basis that is aligned with $\mathrm{\Delta }{\mathbf{r}}_1$ is constructed via:

$${\mathbf{e}}_x^K=\mathrm{\Delta }{\mathbf{r}}_1/\parallel \mathrm{\Delta }{\mathbf{r}}_1\parallel ,{\mathbf{e}}_y^K=\mathrm{\Delta }{\mathbf{r}}_1^{\perp }/\parallel \mathrm{\Delta }{\mathbf{r}}_1^{\perp }\parallel ,$$

(9)

with $\mathrm{\Delta }{\mathbf{r}}_1^{\perp }=({\mathbf{e}}_y\otimes {\mathbf{e}}_x-{\mathbf{e}}_x\otimes {\mathbf{e}}_y)\mathrm{\Delta }{\mathbf{r}}_1$ being $\mathrm{\Delta }{\mathbf{r}}_1$ rotated counterclockwise by a quarter angle. Introducing the decomposition:

$$\mathrm{\Delta }{\mathbf{r}}_2=\mathrm{\Delta }{r}_2^{\parallel }{\mathbf{e}}_x^K+\mathrm{\Delta }{r}_2^{\perp }{\mathbf{e}}_y^K,\mathrm{\Delta }{r}_2^{\parallel }=\frac{\mathrm{\Delta }{\mathbf{r}}_2·\mathrm{\Delta }{\mathbf{r}}_1}{\parallel \mathrm{\Delta }{\mathbf{r}}_1\parallel },\mathrm{\Delta }{r}_2^{\perp }=\frac{\mathrm{\Delta }{\mathbf{r}}_2·\mathrm{\Delta }{\mathbf{r}}_1^{\perp }}{\parallel \mathrm{\Delta }{\mathbf{r}}_1^{\perp }\parallel },$$

(10)

and recognizing $\parallel \mathrm{\Delta }{\mathbf{r}}_1\parallel =\parallel \mathrm{\Delta }{\mathbf{r}}_1^{\perp }\parallel$ , the relative angle can be computed via:

$$\mathrm{\Delta }\varphi ={\tan }^{-1}\left(\frac{\mathrm{\Delta }{r}_2^{\perp }}{\mathrm{\Delta }{r}_2^{\parallel }}\right)={\tan }^{-1}\left(\frac{\mathrm{\Delta }{\mathbf{r}}_2·\mathrm{\Delta }{\mathbf{r}}_1^{\perp }}{\mathrm{\Delta }{\mathbf{r}}_2·\mathrm{\Delta }{\mathbf{r}}_1}\right)={\tan }^{-1}\left(\frac{\mathrm{\Delta }{x}_2\mathrm{\Delta }{y}_1-\mathrm{\Delta }{x}_1\mathrm{\Delta }{y}_2}{\mathrm{\Delta }{x}_2\mathrm{\Delta }{x}_1+\mathrm{\Delta }{y}_2\mathrm{\Delta }{y}_1}\right).$$

(11)

Again, numerically tan⁻¹ is evaluated via the atan2-function. By that, no singularities are present for $|\mathrm{\Delta }\varphi |<\pi$ , 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 discrete 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 $\sqrt{2}\epsilon$ (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=\{0,1,\dots ,N\}$ and $j=\{0,1,\dots ,M\}$ . 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 configuration, these nodal points are arranged symmetrically as depicted in Figure 6.

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.

In the deformed configuration, the Cartesian coordinates of the cell vertices are denoted ${\mathbf{x}}_{(i,j)}=(x_{(i,j)},y_{(i,j)})\in {\mathbb{R}}^2$ . As depicted in Figure 7, the Cartesian coordinates of the kth nodal point within cell (i, j) are ${\mathbf{x}}_{(i,j)}^k=(x_{(i,j)}^k,y_{(i,j)}^k)\in {\mathbb{R}}^2$ . Note that the coordinates of the cell vertices have different denotations 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 coordinates of all nodal points of the sheet:

$$\mathbf{q}=({\mathbf{x}}_{(0,0)},\dots ,{\mathbf{x}}_{(M,N)},{\mathbf{x}}_{(0,0)}^4,\dots ,{\mathbf{x}}_{(0,0)}^{12},\dots ,{\mathbf{x}}_{(M-1,N-1)}^4,\dots ,{\mathbf{x}}_{(M-1,N-1)}^{12})\in {\mathbb{R}}^f.$$

(12)

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.

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 $k_{\mathrm{e}}$ and undeformed length $l_0=\epsilon /(2\cos \beta )=\epsilon /\sqrt{3}$ , where we have used $\beta =\frac{\pi }{4}-\alpha =\frac{\pi }{6}$ . The coordinates required for the kth spring are given by ${\mathbf{q}}_{(i,j)k}^{\mathrm{e}}\in {\mathbb{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 $k_{\mathrm{r}1}$ . The energy of the lth rotational spring is formulated with ${\mathbf{q}}_{(i,j)l}^{\mathrm{r}1}\in {\mathbb{R}}^6$ and the undeformed angle ${\varphi }_0=(-1{)}^l\pi /3$ . The bending stiffness of the fibers within each cell is included by four rotational springs with stiffness $k_{\mathrm{r}2}$ and ${\varphi }_0=0$ influencing the coordinates ${\mathbf{q}}_{(i,j)m}^{\mathrm{r}2}\in {\mathbb{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,\dots ,N-1$ and $j=1,\dots ,M-1$ . These bending stiffnesses are realized by rotational springs with stiffness $k_{\mathrm{r}2}$ and corresponding spring coordinates ${\mathbf{q}}_{(i,j)m}^{\mathrm{r}3}\in {\mathbb{R}}^6$ , which are specified in Figure 8. To extract the spring coordinates from the generalized coordinates (12), the Boolean connectivity matrices ${\mathbf{C}}_{(i,j)k}^{\mathrm{e}}\in {\mathbb{R}}^{4\times f}$ and ${\mathbf{C}}_{(i,j)l}^{\mathrm{r}1},{\mathbf{C}}_{(i,j)m}^{\mathrm{r}2},{\mathbf{C}}_{(i,j)m}^{\mathrm{r}3}\in {\mathbb{R}}^{6\times f}$ are defined by the relations:

$${\mathbf{q}}_{(i,j)k}^{\mathrm{e}}={\mathbf{C}}_{(i,j)k}^{\mathrm{e}}\mathbf{q},{\mathbf{q}}_{(i,j)l}^{\mathrm{r}1}={\mathbf{C}}_{(i,j)l}^{\mathrm{r}1}\mathbf{q},{\mathbf{q}}_{(i,j)m}^{\mathrm{r}2}={\mathbf{C}}_{(i,j)m}^{\mathrm{r}2}\mathbf{q},{\mathbf{q}}_{(i,j)m}^{\mathrm{r}3}={\mathbf{C}}_{(i,j)m}^{\mathrm{r}3}\mathbf{q}.$$

(13)

Table 2.

Assignments of coordinates of nodal points to spring coordinates within cell (i, j).

Spring coordinates	Assignment of nodal points
${\mathbf{q}}_{(i,j)k}^{\mathrm{e}}=({\mathbf{x}}_{(i,j)}^p,{\mathbf{x}}_{(i,j)}^q)$	$\{p,q\}\in \stackrel{k=0}{\{0,5\}},\stackrel{k=1}{\{0,6\}},\stackrel{k=2}{\{5,4\}},\stackrel{k=3}{\{6,4\}},\stackrel{k=4}{\{4,7\}},\stackrel{k=5}{\{4,8\}},\stackrel{k=6}{\{7,1\}},$ $\stackrel{k=7}{\{8,1\}},\stackrel{k=8}{\{4,9\}},\stackrel{k=9}{\{4,10\}},\stackrel{k=10}{\{9,2\}},\stackrel{k=11}{\{10,2\}},\stackrel{k=12}{\{4,11\}},$ $\stackrel{k=13}{\{4,12\}},\stackrel{k=14}{\{11,3\}},\stackrel{k=15}{\{12,3\}}$
${\mathbf{q}}_{(i,j)l}^{\mathrm{r}1}=({\mathbf{x}}_{(i,j)}^p,{\mathbf{x}}_{(i,j)}^q,{\mathbf{x}}_{(i,j)}^s)$	$\{p,q,s\}\in \stackrel{l=0}{\{0,6,4\}},\stackrel{l=1}{\{0,5,4\}},\stackrel{l=2}{\{4,7,1\}},\stackrel{l=3}{\{4,8,1\}},\stackrel{l=4}{\{4,9,2\}},$ $\stackrel{l=5}{\{4,10,2\}},\stackrel{l=6}{\{3,12,4\}},\stackrel{l=7}{\{3,11,4\}}$
${\mathbf{q}}_{(i,j)m}^{\mathrm{r}2}=({\mathbf{x}}_{(i,j)}^p,{\mathbf{x}}_{(i,j)}^q,{\mathbf{x}}_{(i,j)}^s)$	$\{p,q,s\}\in \stackrel{m=0}{\{5,4,9\}},\stackrel{m=1}{\{6,4,10\}},\stackrel{m=2}{\{12,4,8\}},\stackrel{m=3}{\{11,4,7\}}$

Figure 8.

Assignments of coordinates of nodal points to spring coordinates for interaction between cells.

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:

$$\begin{array}{c}\mathcal{E}(\mathbf{q})=\sum _{i=0}^{N-1}\sum _{j=0}^{M-1}(\sum _{k=0}^{15}{\mathcal{E}}^{\mathrm{e}}(k_{\mathrm{e}},l_0,{\mathbf{C}}_{(i,j)k}^{\mathrm{e}}\mathbf{q})+\sum _{l=0}^7{\mathcal{E}}^{\mathrm{r}}(k_{\mathrm{r}1},(-{1)}^l\frac{\pi }{3},{\mathbf{C}}_{(i,j)l}^{\mathrm{r}1}\mathbf{q})\hfill \\ +\sum _{m=0}^3{\mathcal{E}}^{\mathrm{r}}(k_{\mathrm{r}2},0,{\mathbf{C}}_{(i,j)m}^{\mathrm{r}2}\mathbf{q}))+\sum _{i=1}^{N-1}\sum _{j=1}^{M-1}\sum _{m=0}^3{\mathcal{E}}^{\mathrm{r}}(k_{\mathrm{r}2},0,{\mathbf{C}}_{(i,j)m}^{\mathrm{r}3}\mathbf{q}).\hfill \end{array}$$

(14)

Let $\hat{\mathbf{q}}(\epsilon )$ be a function of $\epsilon$ that includes the actual coordinates $\mathbf{q}$ for static equilibrium in the case of $\epsilon =0$ , i.e., $\hat{\mathbf{q}}(0)=\mathbf{q}$ . Then, the variation of the total potential energy $\mathcal{E}$ induced by $\hat{\mathbf{q}}$ is as follows:

$$\delta \mathcal{E}(\mathbf{q})=\frac{\partial \mathcal{E}}{\partial \mathbf{q}}(\mathbf{q})\frac{\mathrm{d}\hat{\mathbf{q}}}{\mathrm{d}\epsilon }(0)={\mathbf{f}}^{\mathrm{T}}(\mathbf{q})\delta \mathbf{q},$$

(15)

where $\delta \mathbf{q}=\frac{\mathrm{d}\hat{\mathbf{q}}}{\mathrm{d}\epsilon }(0)$ are the virtual displacements, and $\mathbf{f}(\mathbf{q})=(\partial \epsilon /\partial \mathbf{q}{)}^{\mathrm{T}}(\mathbf{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:

$$\begin{array}{c}\multicolumn{1}{c}{{\mathbf{f}}^{\mathrm{T}}(\mathbf{q})=\frac{\partial \mathcal{E}}{\partial \mathbf{q}}(\mathbf{q})=}\\ \multicolumn{1}{c}{=\sum _{i=0}^{N-1}\sum _{j=0}^{M-1}(\sum _{k=0}^{15}{\mathbf{f}}^{\mathrm{e}}(k_{\mathrm{e}},l_0,{\mathbf{C}}_{(i,j)k}^{\mathrm{e}}\mathbf{q}){\mathbf{C}}_{(i,j),k}^{\mathrm{e}}+\sum _{l=0}^7{\mathbf{f}}^{\mathrm{r}}(k_{\mathrm{r}1},(-{1)}^l\frac{\pi }{3},{\mathbf{C}}_{(i,j)l}^{\mathrm{r}1}\mathbf{q}){\mathbf{C}}_{(i,j)l}^{\mathrm{r}1}}\\ \multicolumn{1}{c}{+\sum _{m=0}^3{\mathbf{f}}^{\mathrm{r}}(k_{\mathrm{r}2},0,{\mathbf{C}}_{(i,j)m}^{\mathrm{r}2}\mathbf{q}){\mathbf{C}}_{(i,j)m}^{\mathrm{r}2})+\sum _{i=0}^{N-1}\sum _{j=0}^{M-1}\sum _{m=0}^3{\mathbf{f}}^{\mathrm{r}}(k_{\mathrm{r}2},0,{\mathbf{C}}_{(i,j)m}^{\mathrm{r}3}\mathbf{q}){\mathbf{C}}_{(i,j)m}^{\mathrm{r}3}.}\end{array}$$

(16)

Kinematic boundary conditions are imposed by perfect bilateral constraints $0=\mathbf{g}(\mathbf{q})\in {\mathbb{R}}^m$ with a virtual work contribution $\delta W^{\mathrm{c}}=\delta {\mathbf{g}}^{\mathrm{T}}\mathbf{\lambda }=\delta {\mathbf{q}}^{\mathrm{T}}\mathbf{W}(\mathbf{q})\mathbf{\lambda }$ , where $\mathbf{W}(\mathbf{q}{)}^{\mathrm{T}}=\frac{\partial \mathbf{g}}{\partial \mathbf{q}}(\mathbf{q})\in {\mathbb{R}}^{m\times f}$ is the matrix of generalized force directions and $\mathbf{\lambda }\in {\mathbb{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 ${\mathbf{x}}_{i,j}^0$ 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 magnitudes 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 (−5.1 mm) and vertical (−4.35 mm) 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.,

$$\delta {\mathbf{q}}^{\mathrm{T}}(\mathbf{f}(\mathbf{q})+\mathbf{W}(\mathbf{q})\mathbf{\lambda })=0\forall \delta \mathbf{q}\in {\mathbb{R}}^f,$$

(17)

and the constraint equations are $\mathbf{g}(\mathbf{q})=0$ are satisfied. Thus, the equilibrium configuration is determined by the set of nonlinear equations:

$$\mathbf{h}(\mathbf{q},\mathbf{\lambda })=\left[\begin{array}{c}\mathbf{f}(\mathbf{q})+\mathbf{W}(\mathbf{q})\mathbf{\lambda }\hfill \\ \multicolumn{1}{c}{\mathbf{g}(\mathbf{q})}\end{array}\right]=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 identification 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 pantographic sheet using a single set of experimentally optimized parameters.

3.1. Parameter optimization

The yet unknown stiffness parameters $\beta =(k_{\mathrm{e}},k_{\mathrm{r}1},k_{\mathrm{r}2})$ are identified experimentally in the following way and reported in Table 3. Let the Euclidean error between the simulated force $f_{\mathrm{ij}}$ and the measured one ${\overline{f}}_{\mathrm{ij}}$ (see Figure 11) normalized by the maximum force ${\overline{f}}_{i{j}_{\max }}$ within each experiment be denoted as:

$${\delta }_{\mathrm{ij}}=\frac{f_{\mathrm{ij}}-{\overline{f}}_{\mathrm{ij}}}{{\overline{f}}_{i{j}_{\max }}}.$$

(19)

Table 3.

Optimized stiffness parameters and respective relative error.

$k_{\mathrm{e}}$	$k_{\mathrm{r}1}$	$k_{\mathrm{r}2}$	$K({\beta }^{*})$
15.03 Nm−1	14.53 Nm	22.10 Nm	$1.497\times {10}^{-1}\mathrm{N}$

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 ${\beta }^{*}$ were obtained by introducing the cost function:

$$K(\beta )=\frac{1}{2}\sqrt{\sum _{i=1}^N\sum _{j=1}^M\frac{{\delta }_{\mathrm{ij}}^2}{N}},$$

(20)

and solving the optimization problem:

$${\beta }^{*}=\underset{\beta \in {\mathbb{R}}^3}{\mathrm{argmin}}K(\beta ),\mathrm{s}.\mathrm{t}.\mathbf{h}(\mathbf{q},\mathbf{\lambda }\}!)=0,$$

(21)

using a gradient free Nelder–Mead simplex algorithm (see Gao and Han [56]). By heuristic, the initial parameters were set to $k_{\mathrm{e}}=20\mathrm{N}{\mathrm{m}}^{-1}$ , $k_{\mathrm{r}1}=10\mathrm{N}\mathrm{m}$ , and $k_{\mathrm{r}2}=50\mathrm{N}\mathrm{m}$ .

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 ${\beta }^{*}$ varied by $\pm 30\mathrm{\%}$ . On the left-hand side of Figure 9, the loss surface of the cost function is depicted for the optimal parameter $k_{\mathrm{e}}=15.03\mathrm{N}{\mathrm{m}}^{-1}$ but varied $k_{\mathrm{r}1}$ and $k_{\mathrm{r}2}$ . It can clearly be seen that for the optimal value of $k_{\mathrm{e}}$ , a variation of $k_{\mathrm{r}1}$ and $k_{\mathrm{r}2}$ increases the global error. Led by that observation the cost function is evaluated for the optimal values $k_{\mathrm{r}1}=14.53\mathrm{N}\mathrm{m}$ and $k_{\mathrm{r}2}=22.10\mathrm{N}\mathrm{m}$ but varied values of $k_{\mathrm{e}}$ (see right-hand side of Figure 9). Again, for the optimal values, every variation in $k_{\mathrm{e}}$ 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).

Figure 9.

Loss surface for optimal $k_{\mathrm{e}}$ and varied $k_{\mathrm{r}1}$ , $k_{\mathrm{r}2}$ (left). Loss function for fixed $k_{\mathrm{r}1}$ and $k_{\mathrm{r}2}$ but varied $k_{\mathrm{e}}$ (right).

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 rectangle specimen during the extension experiment. It can be observed that during the numerical experiment, some elements of the pantographic sheet intersect each other. Obviously, this is no physical behavior 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 displacement 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.

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.

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.

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 torsional 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 clarify 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 displacements and further derived quantities, e.g., intermediate angles, can be used for a more sophisticated optimization procedure. This requires the application of DIC (cf. [27,60]) or a full-field magnetic resonance imaging (MRI)/computed tomography (CT) scan.