Symmetric-in-Plane Compression of Polyamide Pantographic Fabrics—Modelling, Experiments and Numerical Exploration

: Symmetric in-plane compression of a pantographic lattice structure is modelled and simulated, and the results are compared to previously available experimental data. Said experimental results had shown a peculiar behaviour: depending on the ﬁber density, the deformed shape could present either one or two swellings under compression. The present article is a preliminary modelling attempt aiming at capturing that behaviour numerically.


Model-Driven Design and Synthesis of Pantographic Structures
Pantographic structures are an example of so-called metamaterials [44,45], i.e., they are the result of a model-driven design: given a set of equations, pantographic structures were synthesised in order to follow the very behaviour that is defined by those starting equations. Here, the considered theoretical behaviour is that of an elastic material whose deformation energy may depend only on strain gradient under certain conditions. In other words, pantographic structures are an answer to the synthesis problem of a material that may be modeled as a second gradient continuum.
As a consequence of the peculiar origin of pantographic structures, their study needs a synergy between: first, the starting theoretical model; second, the numerical modelling and simulation of that theoretical model; and last, the manufacturing of, and experimental tests on, real specimen. Regarding that last aspect, studies have been made to understand the influence of certain 3D-printing processes and constitutive materials [46][47][48][49] on the mechanical behaviour of the resulting specimena [50]. Other studies Symmetry 2020, 12 focused on how the design of different mechanical elements inside the structure may affect their global behaviour [51]. Those studies considered different scales, from micrometers [52] to millimetres [53], as well as how the different behaviours at different scales may be linked [54,55]. Studying those various scales involved new methodologies which have been developed for the (micro)mechanics of granular media [56][57][58]. The experimental studies yielded interesting results, such as the ability of pantographic structures to deform elastically for high deformations and how their qualitative behaviour does not depend on their constitutive material [59]. New phenomena have also been observed, like a reversed Poynting effect [60] which led to less known features of the model. There were also advances in the damage mechanics of pantographic structures: their damaging was proven to be highly dependent on their geometry and material properties [61], and their overall resistance to damage was studied by enhanced elastic models [8]. Future prospects concerning the experimental study of pantographic structures involve the promising use of digital image correlation (DIC) techniques as a way to improve the synergy needed for model-driven design [62]. Regarding recent numerical results, pantographic structures have been modelled using finite elements at macroscopic [12] and mesoscopic scales, by using differently built meshes [45,63].

Out-of-Plane and In-Plane Compression of Pantographic Fabrics
In this article, we want to continue these lines of investigations, by attempting the modelling of in-plane compression of thick polyamide pantographic fabrics. Previous litterature in pantographic has extensively dealt with their tensile behaviour [9,38,39]. Further experimental evidence concerned their compresssion behaviour, originating locally in shear tests [8,64]. In such cases, only two layers of fibers were employed and the thickness of the pantographic fabrics was not predominant with respect to other total specimen dimensions. In those cases, cross-section of fibers constituting pantographic fabrics had a very small out-of-plane bending rigidity with respect to the in-plane one. This was due to their low (out-of-plane) thickness. Therefore, pure in-plane compression could not be observed.
By increasing the out-of-plane thickness of the constituting fibers, we avoid out-of-plane buckling, and are able to perform for the first time a purely in-plane compression test. Those experimental results are described in the present work, and a preliminary numerical investigation is reported.

Plan of the Present Work
The plan of the paper is the following. In Section 2, the main features of the used strain gradient continuum model are recalled, as well as its numerical modelling using finite elements. More specifically, a linear model for small perturbations is presented and adpated to higher deformations, by introducing several parameters whose values are to be determined by a numerical-experimental comparison.
In Section 3, the available experimental specimena are presented along the data that are considered to fit the theoretical model. In Section 4, the different steps of the model fitting are described. Using force-displacement graphs and boundary shapes to quantitatively and qualitatively compare the numerical and experimental results, the model parameters are identified and the final model is discussed.
Section 5 summarises the work and gathers the relevant results, before ending with future propspects on how to further the presented work.

Mathematical and Numerical Modelling of Pantographic Structures
In this section the general theoretical framework that will be used to model the studied pantographic structures is presented. The elastic behaviour of the studied pantographic sheet is modeled as a strain gradient continuum. This model does not describe each fiber separately, but considers the global homogenised behaviour instead. We first recall the key features of that model, before presenting its numerical modelling by the finite elements method.

Geometrically Non Linear Second Gradient Elasticity Model
The reference configuration of the studied system is represented by an open domain B 0 ⊂ R 2 , its boundary ∂B 0 (edges) and the boudary of its boundary ∂∂B 0 (wedges). More specifically, B 0 is a rectangle, due to the considered geometry (see Section 3). We consider an orthonormal basis e 1 , e 2 of B 0 . Any generic element of B 0 will be noted X = (X 1 , X 2 ) (only tensor notations are used, but they may also represent the associated vector, linear map or matrix in the basis e 1 , e 2 ).

Kinematics
For all t ∈ ] 0, T], we define the placement χ as: and the respective current configurations B t = χ(B 0 ; t). We also define the displacement u(X; t) = χ(X; t) − X. The placement is assumed to be of class C 1 and locally affine. Its gradient will be noted F( · ; t) = ∇χ( · ; t), so that the Green-Lagrange strain tensor may be expressed as G = 1 2 F t ·F − I . For the sake of notation simplicity and unless specified otherwise, we will now fix a time parameter t ∈] 0, T] and will no longer write it explicitly.

Potential Energy
We define a potential energy functional E taking the displacement field u as its argument. It accounts for both internal and external effects; formally: External Work W ext is the external work. Introducing the surface body force b ext , the linear contact force t ext and double force τ ext , and the vertex forces f ext i , the external work may be rewritten as:

Strain Energy-Linear Elasticity
E int is the elastic deformation energy. It depends on the strain tensor G and its gradient ∇G. We define its density U as follows: The studied pantographic structure is made of two families of beams. It is thus orthotropic, the two privileged directions being those of each family. Moreover, the peculiar geometry of the model makes it invariant for a π/2 rotation, as well as for a mirror reflection. That is to say, it is invariant for the D 4 symmetry group.
Under orthotropy and D 4 symmetry hypotheses, [65] gives the following expression for the linearised strain gradient deformation energy (here, [·] denotes square matrices, · row matrices and {·} column matrices): where the row matrix ε T gathers the values of the linearised strain tensor ε = ε ij e i ⊗ e j : ε T = ε 11 , ε 22 , √ 2 ε 12 η T = ε 11,1 , ε 22,1 , √ 2 ε 12,2 , ε 22,2 , ε 11,2 , √ 2 ε 12,1 and the matrices are defined as: For the linearised case, Placidi et al. [32] derived the values of the [C] and [A] matrices using analytical solutions. They obtained the form: where a lin and c lin ("lin" for "linear case") depend on the geometry and stiffness of the fibers, as well as the density of the fabric (see [32] for the complete derivation): E m being the constitutive material's Young modulus, A m the cross-section area, I m the second moment of area of the cross-section, and d m the distance between two adjacent fibers of the same family.

Strain Energy-Higher Deformations
Although the previous results were derived for the linear case, the same form for the energy will be assumed in the present work. Therefore, the deformation energy to be considered is expressed as: where, this time, the strain terms in ε and η are not linearised: and the coefficients c lin , a lin are substituted with parameters c num , a num whose values should be determined numerically, in order to fit the model with experimental results. Furthermore, the pivots are not considered perfect, and their torsional deformation energy is taken into account by an additional parameter c 33 : Thus, the aim is to build an elastic model for large displacements, starting from a small perturbations model and identifying the defined parameters by fitting them from experimental data.

Variational Approach
The problem to solve is derived using a variational approach. The previously defined potential energy E(u) is assumed to be minimised for the equilibrium displacement field u. Thus, the weak form is obtained as the stationarity condition of the potential energy.
The boundary conditions are: • prescribed null displacement at the lower edge: u |lower = 0 • no applied load on either side: ∂ n u |sides = 0 • prescribed displacement u 0 on the upper edge: u |upper = u 0

Function Space
Since the energy to be minimised involves the second gradient of displacement, the function space should be in H 2 (B 0 ). Thus, the chosen element type is cubic Hermite functions, in order to have C 1 approximations. Moreover, the boundary conditions are taken into account as constraints rather than in the function space.

Modelling with Software
The software Comsol is used for the numerical modelling. A rectangular geometry is defined and meshed (see Section 4). The compresion test is then modelled as a parametric study on the imposed displacement u 0 . The equilibrium is computed for each value of u 0 from the results of the its previous value.
The values of c num , a num and c 33 and their fitting is presented in Section 4, and come from the comparison with actual compression experiments presented in Section 3.

Studied Class of Pantographic Structure
The class of pantographic structures under study is a lattice made of two families of fibers. Every fiber is parallel to the other fibers of the same family, but orthogonal to those of the other family. Each family is planar, and the two considered families are parallel but not in the same plane, so that there is no actual intersection between the two families. However, each family is connected to the other one by what may be called pivots, located wherever the distance between two fibers from different families is minimal. Furthermore, the density of fibers is constant, so that the distance between two adjacent fibers of the same family is constant.

Modelled Specimen: Production and Characterisation
Some key elements are summarised, regarding the experimental specimen to be modelled. The machines and procedures involved are mostly similar to those presented in [51], to which the reader may refer for more details.
The specimen were printed using Powder Bed Fusion technology, more specifically Selective Laser Sintering. A Formiga P100 printer by EOS was used from Warsaw University of Technology, with EOS PA2200 polyamide powder. Every specimen was printed in horizontal configuration, i.e., the pantographic sheets were printed parallel to the support. Each layer was 0.1 mm thick.
The fibers constituting the two printed structures had the following nominal dimensions: The values that are numerically implemented will be presented in the next section (Section 4).

Experimental Procedure
The compression tests were carried out using an MTS Bionix machine in Warsaw University of Technology. One side of the specimen was fixed while the other one's displacement was imposed with a speed of 4 mm/min. The length of the specimen being 21 cm, the speed is around 2% of the length per minute. The measurement sensitivities are 1 N for the force and 0.1 mm for the displacement.
Before the compression tests themselves, each specimen underwent 200 cycles alternating between elongation and compression, in order to eliminate the powder remaining after printing. Those cycles were made with a 1 Hz frequency and 20 mm amplitude.
The presented specimen are visible in Figure 1 in their initial state, and in Figure 2 under compression. reader may refer for more details.
The specimen were printed using Powder Bed Fusion technology, more specifically Selective Laser Sintering. A Formiga P100 printer by EOS was used from Warsaw University of Technology, with EOS PA2200 polyamide powder. Every specimen was printed in horizontal configuration, i.e., the pantographic sheets were printed parallel to the support. Each layer was 0.1 mm thick.
The fibers constituting the two printed structures had the following nominal dimensions: The values that are numerically implemented will be presented in the next section (Section 4).

Experimental Procedure
The compression tests were carried out using an MTS Bionix machine in Warsaw University of Technology. One side of the specimen was fixed while the other one's displacement was imposed with a speed of 4 mm/min. The length of the specimen being 21 cm, the speed is around 2% of the length per minute. The measurement sensitivities are 1 N for the force and 0.1 mm for the displacement.
Before the compression tests themselves, each specimen underwent 200 cycles alternating between elongation and compression, in order to eliminate the powder remaining after printing. Those cycles were made with a 1 Hz frequency and 20 mm amplitude.
The presented specimen are visible in Figure 1 in their initial state, and in

Outputs of the Experiments
The experimental results were used to fit the theoretical model, each specimen having a different purpose.
Specimen 1 was used for quantitative comparison. During its compression, the reaction force was measured. Thus, the outputs of the first experiment were the evolution of the reaction and the imposed displacement. Those results were then exploited to fit the numerical parameters presented in Section 2.
Specimen 2 was used for qualitative comparison. More specifically, it showed a peculiar shape under compression (see Figure 2b), presenting two swellings, while Specimen 1 only had one swelling. Thus, the objective of the second experiment was to verify whether the fitted theoretical model was able to catch the shape of the specimen.

Measurement Uncertainty
Force measurement sensitivity is ∆F = 1 N. Assuming a uniform probability distribution, the type B uncertainty is evaluated as the standard deviation σ F = 1/ √ 3 N.

Outputs of the Experiments
The experimental results were used to fit the theoretical model, each specimen having a different purpose.
Specimen 1 was used for quantitative comparison. During its compression, the reaction force was measured. Thus, the outputs of the first experiment were the evolution of the reaction and the imposed displacement. Those results were then exploited to fit the numerical parameters presented in Section 2.
Specimen 2 was used for qualitative comparison. More specifically, it showed a peculiar shape under compression (see Figure 2b), presenting two swellings, while Specimen 1 only had one swelling. Thus, the objective of the second experiment was to verify whether the fitted theoretical model was able to catch the shape of the specimen.

Measurement Uncertainty
Force measurement sensitivity is ∆F = 1 N. Assuming a uniform probability distribution, the type B uncertainty is evaluated as the standard deviation σ F = 1/ √ 3 N.
Thus, for a measured value F m , the force lies in [F m − ∆F, F m + ∆F], and has 58% probability to be in [F m − σ F , F m + σ F ].

Model Fitting and Numerical Validation
The model presented in Section 2 is now implemented using the software Comsol. We recall that the aim of this preliminary work is to verify whether the presented model can actually be adapated to high deformations. Thus, the following study is rather qualitative, as a more precise quantitative study will be carried out in future works.

Model Fitting and Numerical Validation
The model presented in Section 2 is now implemented using the software Comsol. We recall that the aim of this preliminary work is to verify whether the presented model can actually be adapated to high deformations. Thus, the following study is rather qualitative, as a more precise quantitative study will be carried out in future works.

Fixed Parameters
The numerical compression tests involved two lattices made of the same fibers having the same geometry and made from the same constitutive material (polyamide). The common dimensions and material properties are gathered in Table 1. Both resulting pantographic sheets had 21 cm length and 7 cm width. Table 1. Geometry and material properties of the fibers constituting both tested structures. The only difference between those two structures was the distance between adjacent fibers of the same family. Thus, only the density of fibers characterized each specimen.

Parameters to Be Fitted
The theoretical model involves three parameters c num , a num and c 33 in need of identification in order to fit the model to the experimental results. This identification follows several steps, and uses the compression results of the first experimental specimen: Furthermore, once the identification is done, a qualitative comparison is made with the second specimen's compression, in order to verify if the right boundary shape can be rendered by the model.

First Experiment-Single Swelling
We first perform the parameter identification, using the compression of Specimen 1 (Figure 1a).

Without Torsional Energy
We begin with c num and a num only, leaving a null torsional energy parameter c 33 .

Parameter Values
Those coefficients are derived analytically in the linear case [32]:  Tables 1 and 2). Table 2. Distance d m between two adjacent fibers of the same family.

Force-Displacement Graph Comparison
The graph resulting from the simulation is plotted in Figure 3, as well as the experimental measurements. Measured force error-type is represented by the grey area around the experimental graph.

Force-Displacement Graph Comparison
The graph resulting from the simulation is plotted in Figure 3, as well as the experimental measurements.
Measured force error-type is represented by the grey area around the experimental graph. It is clear that the considered parameter values are not enough to predict the experimental results. The next step is to verify whether changing the value of the torsional parameter c 33 may suffice to fit the model.
One may also notice an apparent softening phenomenon in the numerical model, around an It is clear that the considered parameter values are not enough to predict the experimental results. The next step is to verify whether changing the value of the torsional parameter c 33 may suffice to fit the model.
One may also notice an apparent softening phenomenon in the numerical model, around an imposed displacement of u 0 = 50 mm. This softening is also to be investigated in the next step.

Fitting the Torsional Parameter
Several values for c 33 between 0 and 200 N/m were tested (Figure 4). The value for which the numerical force-displacement graph seemed to be the closest to experimental data was with the corresponding graph in Figure 5. The geometry and mesh of the numerical model is depicted in Figure 6, in the initial state and for an imposed displacement of u 0 = 7 cm. A graphical representation of the fibers was also defined in the numerical model, and can be seen in Figure 7.      The resulting model seems to agree with the experimental results until an imposed displacement of 5 cm, that is to say a 24% deformation. At that point, there is once again an apparent softening, beyond which the theoretical model is clearly non valid.
The actual shape of the experimental specimen at the apparent softening was investigated The resulting model seems to agree with the experimental results until an imposed displacement of 5 cm, that is to say a 24% deformation. At that point, there is once again an apparent softening, beyond which the theoretical model is clearly non valid.
The actual shape of the experimental specimen at the apparent softening was investigated to see whether a contact between fibers may have caused a difference between experimental and  The resulting model seems to agree with the experimental results until an imposed displacement of 5 cm, that is to say a 24% deformation. At that point, there is once again an apparent softening, beyond which the theoretical model is clearly non valid.
The actual shape of the experimental specimen at the apparent softening was investigated to see whether a contact between fibers may have caused a difference between experimental and numerical behaviours. However, no such contact was seen on the experimental specimen, for any value of prescribed displacement. The numerical model was also verified for the same displacement (Figure 8), and no such contact was seen either. Thus, the origin of this apparent softening must still be investigated. Now that the model is fitted, the next step is to verify whether the presence of two swellings instead of a single one may be predicted by the model.

Second Experiment-Two Swellings
The specimen to be modelled is depicted in Figure 1b, and its numerical counterpart in Figure 9 (with a representation of the fibers). The mesh remains the same as previously.

Second Experiment-Two Swellings
The specimen to be modelled is depicted in Figure 1b, and its numerical counterpart in Figure 9 (with a representation of the fibers). The mesh remains the same as previously.

Second Experiment-Two Swellings
The specimen to be modelled is depicted in Figure 1b, and its numerical counterpart in Figure 9 (with a representation of the fibers). The mesh remains the same as previously.

Parameter Values
Since the distance d m between fibers is the only quantity differing Specimen 1 and Specimen 2, c num and a num are recalculated with the right value of d m (see Table 2), while c 33 remains the same.

Qualitative Comparison
The numerical resulting shape is depicted in Figure 10, for an imposed displacement u 0 = 7 cm. Only one swelling is visible, while the experimental results showed two of them (Figure 2b). Thus, for the considered parameter identification, the model does not predict the right boundary shape.
Symmetry 2020, xx, 5 13 of 18 Parameter Values Since the distance dm between fibers is the only quantity differing Specimen 1 and Specimen 2, cnum and anum are recalculated with the right value of dm (see Table 2), while c 33 remains the same.

Qualitative Comparison
The numerical resulting shape is depicted in Figure 10, for an imposed displacement u 0 = 7 cm. Only one swelling is visible, while the experimental results showed two of them (Figure 2b). Thus, for the considered parameter identification, the model does not predict the right boundary shape. However, the double swelling can be observed for the same model, when other arbitrary values are given to the parameters ( Figure 11).
Consequently, while the model is able to render both single and double swellings, the parameter identification still needs to be improved. This discrepancy may be due to the hypotheses on which the parameter values were derived in [32]. More specifically, each fiber is assumed to behave as an Euler beam; this assumption should be investigated to determine whether it may be too restritive to model the compression of the studied lattice structure.

Work Done
The elastic response of a pantographic sheet to a compressive load was modelled. This behaviour was then numerically simulated, and the parameters of the model were fitted to experimental data.
A continuous strain gradient elasticity model was used in the presented work, assuming a quasi-static regime and a bidimensional geometry presenting a D 4 dihedral symmetry. It was derived through a variational approach, by assuming that the equilibrium point minimises the total potential Consequently, while the model is able to render both single and double swellings, the parameter identification still needs to be improved. This discrepancy may be due to the hypotheses on which the parameter values were derived in [32]. More specifically, each fiber is assumed to behave as an Euler beam; this assumption should be investigated to determine whether it may be too restritive to model the compression of the studied lattice structure.

Work Done
The elastic response of a pantographic sheet to a compressive load was modelled. This behaviour was then numerically simulated, and the parameters of the model were fitted to experimental data.
A continuous strain gradient elasticity model was used in the presented work, assuming a quasi-static regime and a bidimensional geometry presenting a D 4 dihedral symmetry. It was derived through a variational approach, by assuming that the equilibrium point minimises the total potential energy. This model has already been used in previous works, more specifically from [32], in a small displacement and deformation framework. It was here adapted to higher deformation cases. Indeed, the deformation energy was supposed to keep the same form for high deformation as for small perturbations. It was then parametrised by coefficients to be fitted numerically to experimental compression results.
Said fitting first involved finite elements modelling: using cubic Hermite elements to keep C 1 functions, the stationary condition of the potential energy was discretised using the software Comsol. The considered boundary conditions were imposed displacements, which were applied as constraints. This numerical model was then simulated, the outputs being the evolution of the reaction force against the imposed displacement, as well as the geometric shape of the boundary.
Regarding the experimental data, compression tests of two specimen were considered. Both were 3D-printed polyamide pantographic sheets, differing only by the density of fibers. Available results were, quantitatively, the evolution of the measured reaction force and imposed displacement, and qualitatively the evolution of the general shape of each specimen.
The experimental and numerical results were then compared in order to identify the parameters of the theoretical model.

Results
First, quantitative comparisons between the numerical and experimental force-displacement graphs led to the identification of the model's parameters. The two parameters c and a which were derived in [32] kept the values they had for small perturbations. The last parameter c 33 was adjusted to fit the numerical and experimental reaction force. The resulting model agrees with the experimental data up to an imposed displacement of 5 cm, i.e., a 24% deformation. An apparent softening phenomenon is then observed in the theoretical model, but does not correspond to any experimental data. In particular, there was no contact between fibers during the compression tests at that displacement. This discrepancy may be due to the hypotheses of the theoretical model. More specifically, the analytically-derived parameters a and c were calculated by considering each fiber to behave as an Euler beam, which may be too restrictive an assumption to properly model the lattice.
However, since this is a preliminary work, we recall that the interest of the previous results lies not in the values themselves, but rather in the possibility of applying the presented method to adapt a linear model to high deformation cases.
Second, a qualitative comparison of the boundary shape revealed a varying number of swellings depending on the considered specimen or parameter values. Indeed, the denser experimental specimen developped two swellings during its compression, while the other only had one. However, the corresponding models showed only one swelling. Nonetheless, the theoretical model does capture the double-swelling phenomenon for other, arbitrary parameter values. It means that, while the parameter identification should be investigated in order to model properly the experimental specimen, the model itself is still able to predict the double-swelling phenomenon.

Possible Future Work
This preliminary work should be followed by more precise studies with respect to both the experimental and numerical values.
Regarding the model, it may be interesting first to look for an analytical explanation of the varying number of swellings, as well as why a softening appears. Second, and as said above, the parameters' identification was made on an Euler beam hypothesis. Studying more specifically the influence of this assumption may reveal more clearly its limits, and enable an exploration of other modelling possibilities.
Concerning the experimental aspect, the accuracy of the data may be improved by other measurement methods, for example by using digital image correlation techniques [62]. Moreover, a better quantification of measurement uncertainty may improve the precision of the experimental data, thus yielding more relevant results for the numerical-experimental comparison.