Discrete versus homogenized continuum modeling in finite deformation bias extension test of bi-pantographic fabrics

A 2D-continuum model describing finite deformations in plane of discrete bi-pantographic fabrics has been recently obtained by applying an asymptotic procedure based on a set of local generalized coordinates. Rectangular bi-pantographic prototypes were additively manufactured by selective laser sintering using polyamide as raw material. Displacement-controlled bias extension tests were performed on such specimens for total elastic deformations up to ca. 25%. Experimental force measurements, complemented by discrete displacement measurements obtained by local digital image correlation, were used to fit the continuum model. In the present paper, a global and minimal set of generalized coordinates, alternative to the one used for the homogenization, is introduced for the discrete model. The mechanical constitutive parameters appearing in the discrete model are then found by means of collected experimental data. Finally, a comparison between experiments, the discrete and the continuum model is presented. It is concluded that (a) the discrete model and the experimental data are in excellent agreement, and that (b) the continuum retains the relevant phenomenology of the discrete system even for a rather low number of cells.

Additively manufactured bi-pantographic rectangular specimen of discrete pantographic beams (see Fig. 1) leading at macroscopic scale to second gradient materials in plane. (For a representative account of second gradient and generalized continua, the reader is referred to [30][31][32][33][34][35][36][37][38][39].) For such materials, the deformation energy density depends upon the second gradient of the deformation. In bi-pantographic structures, such a dependence is given through the rate of change in orientation and stretch of material lines directed along the constituting pantographic beams.
In the spirit of meta-materials [40][41][42], rectangular bi-pantographic specimens shown in Fig. 1 were designed and additively manufactured [43][44][45] by selective laser sintering (SLS) aimed at obeying the theoretical discrete model. Displacement-controlled bias extension tests were performed on these specimens for total elastic deformations up to ca. 25%. Experimental deformation measurements obtained by local digital image correlation (DIC) and force-displacement measurements were used to fit the continuum model.
In Sect. 2, we briefly recall the main results regarding the homogenization of the discrete model, including micro-macro correspondences needed to compare numerical results obtained for descriptions at different scales, i.e., discrete at micro-scale and continuum at macro-scale. The bias extension test considered in the present study is further introduced. Motivated by computational and implementation convenience, we successively introduce in Sect. 3 a minimal set of generalized coordinates for the discrete model, which is alternative to the one used for the homogenization. In Sect. 4, we fit the constitutive parameters of the discrete model by using collected experimental data. We then compare experiments, micro-and macro-modeling. Finally, in Sect. 5, some conclusions are drawn.

Homogenized continuum
In the undeformed configuration (Fig. 2b), the discrete structure is formed by cells arranged upon straight lines within the reference domain , in direction of the unit basis vectors e x , e y ∈ E 2 , being E 2 the two-dimensional Euclidean vector space. The set ⊆ R 2 has boundary ∂ being the disjoint union of N ∈ N smooth line sets ∂ k , k ∈ [1; N ], pairwise intersecting in distinct vertices (see Fig. 2a). The quantity p i, j ∈ R 2 (see Fig. 2b, d) is the current position of the point at position P i, j in the reference configuration. In Fig. 2c, elastic elements are colored in black (extensional Hooke elastic springs, stiffness k E ), red (rotational Hooke elastic springs, stiffness k F ), blue (rotational Hooke elastic springs, stiffness k F ) and green (rotational Hooke elastic springs, stiffness k S ). The homogenization target is a 2D-continuum whose kinematics is characterized by its placement function χ : → R 2 , or, equivalently, by its displacement function u : → E 2 such that u(x, y) = χ(x, y) − x e x − y e y , where x and y denote the coordinates of the undeformed configuration given in the basis {e x , e y }. For the sake of conciseness, we introduce the index α ranging within the set {x, y}. We also introduce the maps ρ α : → R + and ϑ α : → [0, 2π) aimed at rewriting the tangent vector field ∂χ ∂α to deformed material lines directed along e α in the reference configuration as Hence, the quantity ρ α , henceforth referred to as α-stretch, is the length of the vector ∂χ /∂α tangent to the deformed material lines directed along e α in the reference configuration. Letting ε → 0 in Fig. 2b while keeping fixed the overall dimension of the system, and suitably scaling the force elements' stiffnesses, Barchiesi et al. [29] have shown that the deformation energy for the homogenized macro-model reads as being K S , K E and K F the scaled macro-stiffnesses corresponding to the micro-stiffnesses k S , k E and k F , respectively. It is worth recalling [46] that for the above continuum deformation energy we can prescribe conditions at the smooth boundaries on (a) the normal placement gradient, (b) the placement function, and conditions on the placement function at vertices, i.e., singular points of the boundary. To compare the numerical results of the micro-and macro-model, beyond the micro-macro identification of independent kinematic descriptors the following micro-macro correspondences, obtained by neglecting non-leading ε-terms in Taylor expansions of corresponding continuum quantities evaluated at discrete points, are considered For the shearing angle (see Fig. 3), we can identify the micro-macro correspondence π/2 − arccos ∇χ e x · ∇χ e y ∇χ e x ∇χ e y (x,y)=(x i ,y j ) for the undeformed (left) and deformed (right) configurations Introducing for the discrete system the nodal displacements u i, j ∈ E 2 such that u i, j = p i, j − P i, j , by the micro-macro identification (3), we have u(x i , y j ) = u i, j . Therefore, we can define the following micro-macro correspondences at the boundaries (cf. Fig. 4) The bias extension test 1 is commonly used in textile mechanics as a standard experiment to measure the in-plane combined shearing and tensile response of materials made up of two families of fibers. Such a monoaxial extension test is performed on rectangle-shaped specimens having long sides aligned with the loading direction. Constituting fibers-pantographic ones in the present study-are initially oriented at ± 45-degrees to the loading direction. In the present study, the test is modeled by a domain having the shape of a rectangle (with sides 187 mm × 119 mm and γ in Fig. 2b equal to π/6, see Fig. 4) and being subjected to the following essential boundary conditions: Shearing angles at points B and A in Fig. 4 are connected to the transverse contraction of the specimen (hence shearing stiffness) and fiber bending (hence effective fiber bending stiffness), respectively. Indeed, shearing deformation is maximum in B and, since the gray dashed triangle adjacent to the (left) shortest side behaves quasi-rigidly, fiber bending is maximum in A. Fiber extension is maximum at specimen corners-close to corners in the present study-and, therefore, effective fiber extensional stiffness is strongly related to fiber extension at specimen corners. For the discrete system, we have ε=12.02 mm, see Fig. 6. The correspondence Table 1 Boundary conditions for micro-and macro-model between boundary conditions (7) for the continuum model and those for the discrete model are schematized in Table 1.

Micro-model revisited
For solving the discrete micro-model directly and without making any of the hypotheses assumed for the derivation of the continuum model, it is much more convenient to introduce an alternative global, minimal set of generalized coordinates than the one used for the homogenization, i.e., the kinematics of the discrete system is entirely described by the coordinates of the nodal points. The bi-pantographic fabric is modeled as a discrete elastic spring system that is embedded in the two-dimensional 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 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 fabric, which includes some cumbersome but necessary bookkeeping of the relevant degrees of freedom for each spring contribution. Lastly, 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 Fig. 5. The position r = ζ e ζ + ς e ς ∈ E 2 of a typical nodal point is commonly addressed by its Cartesian coordinates x = (ζ, ς ) ∈ R 2 with respect to the orthonormal basis vectors e ζ , e ς ∈ E 2 . If not stated otherwise, R ftuples are considered in the sense of matrix multiplication as R f ×1 -matrices, i.e., as "column vectors". Let q e = (ζ 1 , ς 1 , ζ 2 , ς 2 ) ∈ R 4 be the coordinates of two points interconnected by an extensional spring as depicted in Fig. 5a. Introducing the abbreviations ζ = ζ 2 − ζ 1 and ς = ς 2 − ς 1 , the distance between the two points is The derivative with respect to q e is the "row vector" The potential of an extensional spring with stiffness k e > 0 and undeformed length l 0 > 0 is the function The derivative of the potential (10) with respect to q e is the function Rotational springs are interactions between three nodal points. Let q r = (ζ 1 , ς 1 , ζ 2 , ς 2 , ζ 3 , ς 3 ) ∈ R 6 be the Cartesian coordinates of three nodal points as depicted in Fig. 5b. With the abbreviations ζ 1 = ζ 2 − ζ 1 , ζ 2 = ζ 3 − ζ 2 , ς 1 = ς 2 − ς 1 and ς 2 = ς 3 − ς 2 , the distances between the respective points are with the corresponding derivatives The angles between the e ζ -axis and the vectors ζ 1 e ζ + ς 1 e ς and ζ 2 e ζ + ς 2 e ς , respectively, are introduced by the relations with the corresponding derivatives Note, numerically tan −1 is implemented by an atan2-function with a range (−π, +π] ⊂ R. The potential energy of a rotational spring with stiffness k r > 0 and undeformed angle −π < φ 0 ≤ π is Straightforwardly, the derivative with respect to q r leads to the R 1×6 matrix The bi-pantographic fabric 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  Table 2 Assignments of coordinates of nodal points to spring coordinates within cell (i, j)

Spring coordinates
Assignment of nodal points Besides the 4 vertices, in each cell, there are 9 additional nodal points required to capture the complex behavior of the bi-pantographic fabric. In the reference configuration, these nodal points are arranged symmetrically as depicted in Fig. 6.
In the deformed configuration, the Cartesian coordinates of the cell vertices are denoted x (i, j) = (ζ (i, j) , ς (i, j) ) ∈ R 2 . As depicted in Fig. 7, the Cartesian coordinates of the kth nodal point within cell (i, j) are  1,N −1) , . . . , 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 Fig. 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 fabric are modeled by 16 extensional springs each with stiffness k e and undeformed length l 0 = √ 3 3 ε. The coordinates required for the kth spring are given by q e (i, j)k ∈ R 4 . The resistance of the pins with respect to torsion is captured by 8 rotational springs, depicted in Fig. 7 by green arcs, with stiffness k r1 . The energy of the lth rotational spring is formulated with q r1 (i, j)l ∈ R 6 and the undeformed angle φ 0 = (−1) l π/3. The bending stiffness of the fibers within each cell is included by 4 rotational springs with stiffness k r2 influencing the coordinates q r2 (i, j)m ∈ 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  = 1, . . . , N − 1 and j = 1, . . . , M − 1. These bending stiffnesses are realized by rotational springs with stiffness k r3 = k r2 and corresponding spring coordinates q r3 (i, j)m ∈ R 6 , which are specified in Fig. 8. To extract the spring coordinates from the generalized coordinates (18), the Boolean connectivity matrices Using the spring coordinates of Table 2 and Fig. 8 together with the potential energies (10) and (16), the total potential energy of the discrete bi-pantographic fabric is Letq(ε) be a function of ε that includes the actual coordinates q for static equilibrium in the case of ε = 0, i.e.,q(0) = q. Then, the variation of the total potential energy E induced byq is where δq = dq dε (0) are the virtual displacements, f(q) = (∂E/∂q) T (q) are the internal generalized forces and the transposed of a matrix is indicated by (·) T . Using (11) and (17) together with the total potential energy (20), the internal generalized forces of the bi-pantographic fabric are obtained by Kinematic boundary conditions are imposed by perfect bilateral constraints 0 = g(q) ∈ R m with a virtual work contribution δW c = δg T λ = δq T W(q)λ, where W(q) T = ∂g ∂q (q) ∈ R m× f is the matrix of generalized force directions and λ ∈ R m the vector of constraint forces. The discrete system is now in static equilibrium if and only if the total virtual work of internal generalized forces and constraint forces vanishes for all virtual displacements, i.e., δq T (f(q) + W(q)λ) = 0 ∀δq ∈ R f (23) and the constraint equations g(q) = 0 are satisfied. Thus, the equilibrium configuration is determined by the set of nonlinear equations f(q) + W(q)λ g(q) = 0 ( 2 4 ) which can be solved, at least locally, by a Newton-Raphson iteration scheme.

Results
Specimens were 3D-printed using selective laser sintering (SLS). Polyamide powder PA2200 was used as raw material. A top view of the manufactured specimen, with a zoom on relevant details, is shown in Fig. 1. The micro-model's kinematics prescriptions in the third and sixth rows of Table 1 are realized by connecting adjacent hinge axes in the vicinity of gripping areas, see Fig. 1. Stocky rhomboidal elements are used to this end, which can be assumed to be rigid for the studied load cases. An increasing displacement was prescribed horizontally on the right side of the specimen with a loading rate of 15 mm/min meant to be quasi-static in relation to the weight of the specimen. Kinematic experimental acquisition was achieved by non-stereo digital image correlation (DIC). Pictures of the surface during deformation were acquired (0.5 fps, i.e., nominally for 1 mm displacement increments) by means of a Canon EOS 600D camera with a definition of 4272 × 2848 pixels and an 8-bit dynamic range. The sought kinematics is given by the in-plane displacements of the hinges at positions p i, j 's of the bi-pantographic structure. The displacements of these (discrete) points are retrieved via local DIC (i.e., using zones of interest or ZOIs [47] centered on each hinge). The simplest approach seeks the rigid body translation of each considered ZOI, as originally performed in particle image velocimetry [48,49]. An in-depth illustration of the application of the above-mentioned DIC techniques to the present study is beyond the scopes of this article. Therefore, it is omitted. For more details on the design, manufacturing and experimental testing of bi-pantographic specimens, we refer to [29].
Parameters for the discrete (k F , k E and k S ) model and the continuum (K F , K E and K S ) one were independently found by fitting the three curves (see Fig. 9). The first one ( Fig. 9 (left)) is the total reaction force versus u. The second one ( Fig. 9 (center)) is the shearing angle at point A (cf. Fig. 4) versusū. Finally, the third one ( Fig. 9 (right)) is the shearing angle at point B, which is strongly related to the transverse contraction of the specimen, againstū.
Parameters' value for discrete and continuum models as found by way of the fitting procedure described above is reported in Table 3. Deformed configurations as predicted by modeling, i.e., p i, j 's for the discrete Fig. 10 Deformed configurations as predicted by the discrete and continuum models, i.e., p i, j 's for the discrete model and χ (x i , y j ) for the continuum one, for differentū's are compared with experimentally measured ones. Abscissas and ordinates are expressed in mm model and χ(x i , y j ) for the continuum model, are compared for differentū's with experimental data, and are plotted in Fig. 10. Discrete and continuum modeling substantially overlap, and the agreement with experiments is remarkably good. It is worth noticing that experimental data are noticeably not-completely symmetric, which might be a point to address in future. A contour plot of the y-stretch ρ y is shown in Fig. 11 for discrete and continuum models. Again, micro-macro agreement is very good even for such a rather small number of cells.

Conclusion
In the present article, we performed a preliminary comparison between discrete and homogenized continuum modeling in the finite deformation regime analysis of a recently designed specimen with bi-pantographic micro-structure. From obtained evidences, we conclude that the metamaterials synthesis problem has been successfully accomplished and that the agreement between homogenized continuum and discrete system is very good already for low numbers of cells. Possible outlooks include the investigation of lower-scale approaches based on granular micro-mechanics [50][51][52], the use of semi-discrete descriptions in which continuous beam elements replacing assemblies like those in Fig. 2d are solved by iso-geometric analyses [53][54][55][56][57][58][59] and the study of different boundary prescriptions [60]. Fig. 11 y-stretch ρ y for discrete and continuum models. Abscissas and ordinates are expressed in mm