Malo Valmalle

Benjamin Smaniotto

Mario Spagnuolo
email: mario.spagnuolo@unica.it

Alessandro Ciallella

François Hild

Mesoscale DVC analyses and parameter calibration for pantographic block in 3-point flexure

Keywords:

In the present work, the mechanical response of a fiber network metamaterial was studied when subjected to 3-point flexure. To account for the initial deformation of the sample due to the fabrication process, backtracking (i.e., the nominal mesh was accurately repositioned in the reference configuration) was implemented in the Digital Volume Correlation (DVC) procedure. Kinematic fields were then measured via DVC at two different scales. Finite element based DVC was validated thanks to a series of initialization steps. The displacement fields were successfully measured and a peculiar deformation mechanism was observed. These measurements were employed to calibrate the constitutive parameters of a second gradient continuum model introduced for modeling this type of metamaterial. Numerical simulations are shown to be in good agreement with experimental measurements at the macroscopic scale.

Introduction

The advancement of additive manufacturing techniques has made it possible to develop metamaterials whose inner architecture may be very complex [START_REF] Schaedler | Architected cellular materials[END_REF][START_REF] Pfaff | Resource analysis model and validation for selective laser melting, constituting the potential of lightweight design for material efficiency[END_REF]. This complexity determines their mechanical properties, which can be designed to be highly-performing and not achievable with standard materials. Such remarkable properties may persist when new mesostructures are embedded in a matrix. The obtained composite material, which can be regarded as a fiber reinforced matrix [START_REF] Franciosi | Uniformity of the Green operator and Eshelby tensor for hyperboloidal domains in infinite media[END_REF][START_REF] Schulte | Isogeometric analysis of fiber reinforced composites using Kirchhoff-Love shell elements[END_REF], still maintains some of the peculiarities of the underlying metamaterial architecture.

Whether considered as an (embedded) phase of a composite material or as a material per se, the mechanics of metamaterials is often modeled in the context of generalized continua models [START_REF] Alibert | Truss modular beams with deformation energy depending on higher displacement gradients[END_REF]. In the class of mechanical metamaterials, a very sophisticated one is represented by fibrous arrangements (i.e., architectures composed of fiber layers variously interconnected [START_REF] Ciallella | A rate-independent internal friction to describe the hysteretic behavior of pantographic structures under cyclic loads[END_REF][START_REF] Spagnuolo | Contact interactions in complex fibrous metamaterials[END_REF]). When designing fibrous metamaterials, a fundamental property to exploit is the fact that the morphology of their mesostructure significantly influences their macroscopic response [START_REF] Yang | Material characterization and computations of a polymeric metamaterial with a pantographic substructure[END_REF]. These properties are controlled to a lesser extent by the microscopic properties of the constituent material [START_REF] Spagnuolo | Phenomenological aspects of quasi-perfect pivots in metallic pantographic structures[END_REF]. Therefore, by working only on the geometric arrangement of basic structural elements, and without changing the constitutive material, it is possible to obtain very different macroscopic responses [START_REF] Giorgio | Numerical identification procedure between a micro-Cauchy model and a macro-second gradient model for planar pantographic structures[END_REF][START_REF] Boutin | Linear pantographic sheets: asymptotic micro-macro models identification[END_REF].

There exist natural materials in which the aforementioned principle is exploited. Remarkable examples of natural fibrous architectures are found in bone tissues or woods. In the case of bone tissues, the architecture of the fiber network is obtained through the process of topological optimization [START_REF] Desmorat | Topology optimization in damage governed low cycle fatigue[END_REF][START_REF] Giorgio | Lattice shells composed of two families of curved Kirchhoff rods: an archetypal example, topology optimization of a cycloidal metamaterial[END_REF] induced by natural evolution. In this way, the resulting bone microstructure is the best possible for finding a material having optimal ratios between resistance and weight [START_REF] George | A multiphysics stimulus for continuum mechanics bone remodeling[END_REF][START_REF] Giorgio | On mechanically driven biological stimulus for bone remodeling as a diffusive phenomenon[END_REF].

In this paper, results are presented concerning a 3-point flexural test performed on a specimen made of a fibrous architecture known as pantographic metamaterial. The simplest example of such metamaterial is a sheet composed of two mutually orthogonal fiber layers mechanically interconnected by cylindrical joints. These joints are generally referred to as pivots or hinges in the literature [START_REF] Alibert | Truss modular beams with deformation energy depending on higher displacement gradients[END_REF][START_REF] Dell'isola | Large deformations of planar extensible beams and pantographic lattices: heuristic homogenization, experimental and numerical examples of equilibrium[END_REF][START_REF] Giorgio | Two layers pantographs: a 2D continuum model accounting for the beams' offset and relative rotations as averages in SO (3) Lie groups[END_REF]. They either have some torsional rigidity, and thus add a deformation mechanism to the fibers, or behave as perfect hinges, thereby only ensuring that the families of fibers remain mutually bonded, but exerting no constraint on the relative rotations between them [START_REF] Valmalle | Local-global DVC analyses confirm theoretical predictions for deformation and damage onset in torsion of pantographic metamaterial[END_REF]. The relevance of such joints is that, in the macroscopic response, they are associated with macroshear contributions to the strain energy to be used in the continuum model, and were proven to determine the overall response of the metamaterial [START_REF] Spagnuolo | Phenomenological aspects of quasi-perfect pivots in metallic pantographic structures[END_REF][START_REF] Spagnuolo | Mesoscale modeling and experimental analyses for pantographic cells: effect of hinge deformation[END_REF][START_REF] Spagnuolo | Out-of-plane deformation reduction via inelastic hinges in fibrous metamaterials and simplified damage approach[END_REF]. It is worth noting that macroscale models used for describing pantographic sheets are based on second gradient 2D continua [START_REF] Mindlin | Second gradient of strain and surface-tension in linear elasticity[END_REF][START_REF] Pideri | A second gradient material resulting from the homogenization of an heterogeneous linear elastic medium[END_REF][START_REF] Luscher | A second gradient theoretical framework for hierarchical multiscale modeling of materials[END_REF][START_REF] Turco | Hencky-type discrete model for pantographic structures: numerical comparison with second gradient continuum models[END_REF][START_REF] Germain | The method of virtual power in the mechanics of continuous media, I: Second-gradient theory[END_REF].

The specimen considered herein is called a pantographic block, namely, it consists of 𝑛 pantographic sheets interconnected with hinges identical to those used inside each sheet. The design of the interconnecting hinges has been conducted in order to minimize any mechanical coupling phenomena inducing deformations of the pantographic block crosssections (see Figure 1). Therefore, in the considered test, one may use 2D second gradient continua models accounting for the thickness of the block via a suitable increase of model stiffnesses. Special attention is to be devoted to the deformation mechanisms of such blocks when undergoing flexure. The study of deformations at the level of the mesostructure of the pantographic metamaterial is crucial. From the point of view of continuum models describing their mechanical response, the related flexural energy determines the part of macrostrain energy depending on the second gradient of macroplacement fields [START_REF] Yang | Micromechanics based second gradient continuum theory for shear band modeling in cohesive granular materials following damage elasticity[END_REF][START_REF] Auffray | Analytical continuum mechanics à la Hamilton-Piola least action principle for second gradient continua and capillary fluids[END_REF][START_REF] Placidi | Identification of twodimensional pantographic structure via a linear D4 orthotropic second gradient elastic model[END_REF][START_REF] Barchiesi | Pantographic beam: a complete second gradient 1D-continuum in plane[END_REF]. The aim is then to observe experimentally the effects due to second gradient energy terms, in a mechanical test that reduces the strains relative to the elongation of the fibers (contributing to first gradient strain energy). In the present study, an in situ 3-point flexural test was monitored via X-ray tomography [START_REF] Maire | Quantitative X-ray tomography[END_REF], and Digital Volume Correlation was utilized to measure displacement fields. As for Digital Image Correlation [START_REF] Sutton | Computer vision-based, noncontacting deformation measurements in mechanics: A generational transformation[END_REF], there are essentially two types of approaches to measure volumetric displacements, namely, local DVC in which small and independent subvolumes are registered [START_REF] Bay | Digital volume correlation: three-dimensional strain mapping using X-ray tomography[END_REF][START_REF] Bay | Methods and applications of digital volume correlation[END_REF] and global (or FE-based) DVC [START_REF] Roux | Three dimensional image correlation from X-Ray computed tomography of solid foam[END_REF] for which the whole region of interest is meshed and correlated in a single analysis (Buljac, Jailin, Mendoza, Neggers, Taillandier-Thomas, Bouterf, Smaniotto, [START_REF] Buljac | Digital Volume Correlation: Review of Progress and Challenges[END_REF]. In the present case, FEbased DVC was utilized and a unique mesh was constructed for the whole metamaterial. FE-based DVC was already used to quantify the deformation of pantographic metamaterials subjected to torsion [START_REF] Auger | Poynting Effects in Pantographic Metamaterial Captured via Multiscale DVC[END_REF][START_REF] Valmalle | Local-global DVC analyses confirm theoretical predictions for deformation and damage onset in torsion of pantographic metamaterial[END_REF]. However, to the authors' best knowledge, it was not yet applied to any flexural test on such metamaterials. Moreover, for the first time, the experimental evidence shown in the present work has made it possible to observe phenomena that one could only expect on a theoretical level. These phenomena were not previously observed; in fact, thanks to tomographic measurements and subsequent DVC analyses, it was possible to directly observe the deformations of the connecting hinges, which are normally hidden by the fibers. The observation of these deformations is crucial. In fact, it allows to validate the hypothesis that the pantographic block can be regarded as a superposition of pantographic sheets, which can be described by means of a 2D continuum model.

In fact, the models used in the description of this type of metamaterials involve terms that depend on second order gradient of the displacement field in the strain energy. For this reason, it is useful to have precise experimental information on such second gradient terms. These second gradient terms are only related to fiber bending. This is the reason why it is useful to isolate this deformation mechanism in order to obtain a study free from quantities linked to other deformation mechanisms. Typically, in a biased extension test (i.e., biased with respect to the fiber directions, which are tilted by ±45 • with respect to the horizontal direction) other types of deformation modes emerge. We refer, for example, to the elongation of fibers or the torsion of the connecting elements (hinges). Much greater energy contributions correspond to these mechanisms than to those of the bending of the fibers (typically, this is the same ratio as between extensional and flexural rigidities of a beam). It is therefore essential to devise methods to isolate the contribution of fiber bending. It was shown how to eliminate torsion contributions of hinges by replacing them with perfect pivots [START_REF] Spagnuolo | Circuit analogies in the search for new metamaterials: Phenomenology of a mechanical diode[END_REF]. In the present work, a 3point flexural test was designed to eliminate, or minimize, the contributions of fiber elongation by only inducing fiber bending.

The outline of the paper is as follows. First, the in situ test is presented. Then, the different steps of FE-based DVC are introduced. The measured displacement fields, corresponding strain and gray level residual fields are analyzed to obtain quantitative information on the various deformation mechanisms at play in such test. Last, the modeling choices are described in which the assumption of nearly undeformable sections was probed.

In situ 3-point flexural test

The sample used herein, see Figures 1 and2 The 3-point flexural test was monitored using computed microtomography [START_REF] Maire | Quantitative X-ray tomography[END_REF]. Such imaging technique is non-intrusive and provides 3D images of samples tested in situ [START_REF] Buffière | In Situ Experiments with X ray Tomography: an Attractive Tool for Experimental Mechanics[END_REF]. Special 3D printed supports, made of ABS, were utilized (Figure 2(a)) to enable the total length and height of the sample to be imaged (i.e., they had virtually the same X-ray absorption as the tested sample). Nylon screws were used to fix the supports to the aluminum alloy platens (Figure 3(a)). The load was applied by the longitudinal actuator moving the bottom platen of the in situ testing machine, which was controlled in position (with steps of 5 mm). The test was conducted up to the stroke limit (i.e., 20 mm) of the machine. Cardboard was added to avoid friction and locking of the beam ends with respect to the supports (Figure 3). Six tomographic scans were acquired, namely, two in the reference configuration for uncertainty quantification, and four in deformed configurations. The hardware parameters of the in situ setup are gathered in Table 1. Once cropped, the reconstructed volumes covered 50 × 50 × 100 mm 3 with a 100 µm/vx resolution.

Series of 1,500 radiographs per scan were acquired to reconstruct 3D images of the deformed sample (Figure 4). This large number of radiographs per scan was needed since most of the sample occupied the whole length of the detector (Table 1). Overall, the quality of the reconstructed volumes was deemed very satisfactory for DVC calculations. Figure 5 shows the plot of the reaction force versus prescribed displacement. This curve appears to be essentially linear. Four "discontinuity zones" are also observed, which correspond to the four tomographic acquisitions. They reflect small relaxations of the material that occurred as a result of test interruption and the effects of minute misalignment due to the rotation of the turntable. The reaction force plot is of crucial significance in this study. As briefly addressed in the introduction, pantographic metamaterials are describable as generalized (i.e., second gradient) continua. In the latter ones, three contributions are introduced to the strain energy of the metamaterial, namely, one related to the elongation of the fibers, another one to the flexure of the fibers, and a third one (if any) to the torsion of the hinges. Specifically, for modeling the contribution related to flexure of the fibers, one has to introduce the second gradient of the placement field. This term is not present in theories based on classical Continuum Mechanics and, therefore, it is crucial to study it in order to fully characterize such metamaterials [START_REF] Greco | In plane shear and bending for first gradient inextensible pantographic sheets: numerical study of deformed shapes and global constraint reactions[END_REF]. In the present case, for the maximum prescribed displacement, the energy contribution associated with the elongation of the fibers was expected to be negligible [START_REF] Spagnuolo | Circuit analogies in the search for new metamaterials: Phenomenology of a mechanical diode[END_REF], while the contributions associated with fiber flexure and torsion of the hinges were present. The measured reaction force was therefore the combined signature of both deformation mechanisms. Since the measured reaction force was essentially linear with respect to the prescribed displacement, it was expected that nonlinear deformation mechanisms (i.e., plasticity and damage) remained very limited. This point will be further discussed when local strain fields will be reported.

By comparing with theoretical results [START_REF] Spagnuolo | Circuit analogies in the search for new metamaterials: Phenomenology of a mechanical diode[END_REF], it was observed that the measured reaction force was significantly higher than expected for this type of material. Two substantial differences are highlighted with respect to the investigated case: i. in the present case, the metamaterial was not a single pantographic sheet, but consisted of a multilayered network, and therefore one must take into account the combined actions of various layers of bent fibers;

ii. the hinges were not perfect but deformed, and as widely discussed in the literature, most of the deformation energy of pantographic metamaterials is due to the deformation of hinges [START_REF] Giorgio | Numerical identification procedure between a micro-Cauchy model and a macro-second gradient model for planar pantographic structures[END_REF][START_REF] Spagnuolo | Phenomenological aspects of quasi-perfect pivots in metallic pantographic structures[END_REF].

Different DVC steps

In the following analyses, two different scales were considered. Macroscale DVC dealt with a mesh for which all the morphological details of the metamaterial were not accounted for. It was useful to initialize mesoscale analyses (in which the actual architecture of the metamaterial was meshed) given the large displacement amplitudes between each acquisition (Figure 4). The macroscale measurements were subsequently compared to model predictions at the same scale. The mesoscale enabled deformation mechanisms to be quantified for the hinges. Various DVC steps had thus to be tailored to the present case.

Mesh in the reference configuration

FE-based DVC [START_REF] Roux | Three dimensional image correlation from X-Ray computed tomography of solid foam[END_REF][START_REF] Hild | Toward 4D Mechanical Correlation[END_REF] was used in the present analyses. Such approach requires an FE mesh to be constructed. In the present case, the starting point was the STL model of the to-be-printed pantographic block. This initial information corresponded to the nominal configuration of the digital twins (Figure 6). A 3D mesh made of 4noded (T4) tetrahedra was constructed from the STL model with the PDE toolbox of Matlab. Once generated, the mesh was cropped in order to measure displacement fields only between the two outer supports, and for the mesh to remain within the monitored volume when deformed. The voxelized twin, which was needed for DVC purposes, was a mask in which each voxel lying within the mesh was assigned a nonzero value. It was also constructed from the STL model of the nominal configuration (Figure 6). The STL model was cropped (with Avizo Fire 8.0) to be consistent with the DVC ROI, and to coincide with the size of the used mesh. With this cropped STL file, a ray tracing algorithm [START_REF] Patil | Voxel-based representation, display and thickness analysis of intricate shapes[END_REF] was used to create the voxelized twin. The mesh initially contained 62,270 nodes and 177,067 T4 elements whose mean size (measured as the cube root of the average elementary volume) was 6 vx. Another issue appeared, namely, the PDE mesher generated few very small elements (with volumes of the order of one voxel) that would have resulted in a very poor conditioning of the Hessian matrices used in DVC calculations. The choice was made to remove the elements whose size was less than 2 vx. This slightly modified mesh contained 62,261 nodes and 171,755 T4 elements whose mean length still was 6 vx. All hinges were explicitly meshed. Conversely, the geometric details of each beam close to the connections were not as fine (Figure 7(b-d)). This choice resulted from compromises between mesh fineness, voxel resolution and size of the DVC (Hessian) matrix.

Because of small errors due to printing, mispositioning of the sample and the application of a preload, the experimental reference configuration did not coincide with that of the digital twin. Therefore, a backtracking procedure [START_REF] Auger | Poynting Effects in Pantographic Metamaterial Captured via Multiscale DVC[END_REF] was needed to fit the mesh constructed with the digital twin onto the printed configuration. Another route could have been to resort to so-called shape correlation [START_REF] Kleinendorst | Mechanical Shape Correlation: A novel integrated digital image correlation approach[END_REF]. However, in the present case, standard DVC was sufficient with no need for mechanical regularization. The procedure consisted in registering the printed configuration with the voxelized digital twin (Figure 6). An auxiliary mesh encompassing the volume in the printed configuration was considered and DVC was run. It is worth noting that this mesh was also used in the DVC analyses at the macroscale.

Figure 9 shows the displacement field that had to be applied between the printed and nominal configurations. The rigid body translations were taken out since they were due to the fact that the reconstruction frame and that of the digital twin did not coincide. The displacement fields along the vertical and longitudinal directions were consistent with a flexural preload that was applied. The displacement

Uncertainty quantification

The reconstructed volumes were registered using the Correli 3.0 framework (Leclerc, Neggers, Mathieu, Hild and Roux, 2015) in which Hencky regularization (dell'Isola, Seppecher, Spagnuolo, Barchiesi, Hild, Lekszycki, Giorgio, Placidi, Andreaus, Cuomo, Eugster, Pfaff, Hoschke, Langkemper, Turco, Sarikaya, Misra, De Angelo, D'Annibale, Bouterf, Pinelli, Misra, Desmorat, Pawlikowski, Dupuy, Scerrato, Peyre, Laudato, Manzari, Göransson, Hesch, Hesch, Franciosi, Dirrenberger, Maurin, Vangelatos, Grigoropoulos, Melissinaki, Farsari, Muller, Abali, Liebold, Ganzosch, Harrison, Drobnicki, Igumnov, Alzahrani and Hayat, 2019b) was implemented for the mesoscale analyses (Tab. 2), since virtually no contrast existed within the printed parts (except few porosities), which prevented analyses to be performed at even lower scales. In the present case, the measured displacement field 𝐮 was parameterized with the nodal displacements associated with the finite element discretization based on 4noded (T4) tetrahedra [START_REF] Hild | Toward 4D Mechanical Correlation[END_REF]. At this stage, the regularization length was not yet selected. Therefore, the uncertainties were studied for different regularization lengths 𝓁 𝑚 with the same mesh. The noise floor levels were assessed using the two scans of the reference configuration. DVC analyses were run with the backtracked mesh (Figures 8). Rigid body motions were subtracted from the measured displacement fields, the standard deviations of each nodal displacement component was estimated and are displayed in Figure 10(a). The measurement uncertainties were very close for the three spatial directions for a given regularization length. The (classical) compromise between measurement uncertainty and regularization length was found [START_REF] Leclerc | Digital volume correlation: What are the limits to the spatial resolution[END_REF]). However, it was observed that the uncertainty level was not solely controlled by acquisition noise since the power -3∕2 dependence of the displacement uncertainty was not followed. This result proves that other sources of error also occurred, which was expected when dealing with tomographic acquisitions [START_REF] Buljac | Digital Volume Correlation: Review of Progress and Challenges[END_REF]. For large regularization lengths, standard displacement uncertainties as low as 10 -2 vx were reached thanks to the regularization strategy used herein.

The maximum principal strain 𝜖 1 and the von Mises equivalent strain 𝜖 𝑣𝑀 will be reported in the sequel. The latter one was considered to assess the level of plastic strains in the hinges and beams. Their standard uncertainties were therefore evaluated as well. They correspond, for each strain descriptor, to the standard deviation of the elementary quantities assessed in the whole mesoscale mesh. The larger the regularization length, the lower the standard strain uncertainties (Figure 10(b)). Apart from the second principal strain (𝜖 2 ), all three strain components (𝜖 1 , 𝜖 3 and 𝜖 𝑣𝑀 ) had virtually identical uncertainties. When large regularization lengths were selected, the strain uncertainties were compatible with the measurement of elastic strains. Conversely, when less regularized, DVC analyses will not allow elastic strains to be measured accurately. Conversely, plastic strains could be quantified even for very small regularization lengths (for levels greater than 10 -2 ). Last, let us note that the power law interpolation with an exponent of -5∕2 was not followed for the same reasons as those discussed for the displacement uncertainties. Last, the uncertainties were evaluated for the 1,935 hinges that were explicitly meshed (Figure 7(b-d)). The axial direction of the hinges corresponded to the 𝑦-direction in the nominal configuration. There are two nodal planes connecting each hinge to the beams (i.e., bottom and top sections). For each of them, the mean displacements and rotations were computed. The displacement jumps were obtained by subtracting the mean displacement of the bottom section of the hinge to that of the top section (Figure 7(d)). Similarly, the rotation jumps were assessed by subtracting the mean rotation of the bottom section to that of the top section. In the following, the axial components of the rotation ([[𝜃 𝑎 ]]) and displacement ([[𝑢 𝑎 ]]) jumps were computed in addition to the corresponding transverse components (i.e., the norm of the two components in the plane normal to the axial direction). Figure 11 shows their standard uncertainty levels as functions of the regularization length 𝓁 𝑚 . For the displacement jumps, the uncertainties were lower than those of the nodal displacements (Figure 10(a)). For the rotation jumps, the corresponding uncertainties were very close to those of the minor and major principal strains (Figure 10(b)).

Initialization

Due to the large deformations between consecutive scans (Figure 4), the DVC calculations had to be initialized. A digital image correlation (DIC) algorithm was used on each nodal cross-section of the macromesh encompassing the sample (i.e., the auxiliary mesh of the backtracking procedure, see Figure 8) seeking in-plane rigid body motions (i.e., assuming that there was no out-of-plane motion. Because the top grip was motionless in this experiment, the middle section of the pantograph was expected to be essentially motionless. Accordingly, the DIC calculations were The deformed configurations show that rather large displacements occurred in the (longitudinal) 𝑥-direction (Figure 13) as well. They were induced by large rotations that took place in the vicinity of the outer supports. Such displacements were not pre-corrected by the initialization procedure, contrary to those in the 𝑧-direction.

Direct DVC calculations

The displacement fields measured on the macroscopic mesh (Figures 12 and13) were used to initialize DVC analyses at the mesoscale. They were run directly for the 4 deformed scans using the backtracked mesh (Figure 8). The DVC calculations were stopped when the L2-norm of displacement corrections were less than 10 -2 vx. All calculations converged quickly. Had the DVC analyses not been initialized, the convergence would have been very slow, if possible.

Since regularized DVC was selected at the mesoscale, the measured displacement field 𝐮 was determined by minimizing the sum of squared gray level differences

Φ 2 𝑐 ({𝜐 𝜐 𝜐}) = 1 |ROI| ∑ ROI ( 𝐼 0 (𝐱) -𝐼 𝑡 (𝐱 + 𝐮(𝐱, {𝜐 𝜐 𝜐})) ) 2 (1)
where the gray level residual 𝜑 𝑐 at each voxel 𝐱 reads

𝜑 𝑐 = 𝐼 0 (𝐱) -𝐼 𝑡 (𝐱 + 𝐮(𝐱)) (2)
penalized by the equilibrium gap

Φ 2 𝑚 ({𝜐 𝜐 𝜐}) = {𝜐 𝜐 𝜐} ⊤ [𝐊] ⊤ [𝐊]{𝜐 𝜐 𝜐} (3) such that {𝜐 𝜐 𝜐} meas = arg min {𝜐 𝜐 𝜐} ( Φ 2 𝑐 ({𝜐 𝜐 𝜐}) + 𝑤 𝑚 Φ 2 𝑚 ({𝜕𝜐 𝜐 𝜐}) ) ( 4
)
where {𝜕𝜐 𝜐 𝜐} denotes the column vector of incremental nodal displacements from one analysis to the next, [𝐊] the rectangular stiffness matrix associated with bulk and free surface nodes, 𝐼 0 the volume in the reference configuration, 𝐼 𝑡 the volume in the deformed configuration, and ROI the selected region of interest. The regularization weight 𝑤 𝑚 is proportional to the regularization length 𝓁 𝑚 raised to the power 4 (Taillandier-Thomas, [START_REF] Taillandier-Thomas | Localized strain field measurement on laminography data with mechanical regularization[END_REF]. It is worth mentioning that since regularization was performed at the mesoscale, the stiffness matrix is associated with a standard (i.e., Cauchy) continuum formulation. Figure 14(a) shows the change of the root mean square (RMS) gray level residual Φ 𝑐 as a function of the equilibrium gap Φ 𝑚 for five different regularization lengths. In addition to the four acquisitions under different deflections, a fifth result was added, which corresponds to that of the uncertainty quantification for which a second scan was acquired in the reference configuration. This first case always led to the lowest gray level residuals. The fact that the gray level residuals remained rather close to these minimum levels when the sample was deformed is an indication that the registrations were successful. When the regularization length was increased, the equilibrium gap Φ 𝑚 decreased as more weight was put on the penalty term. Conversely, the gray level residuals increased as the sample was more deformed for large regularization lengths. This trend indicates that the regularization was too strong for the largest regularization lengths (i.e., Hencky elasticity was no longer sufficient to globally account for the underlying behavior). One of the criteria used to select the regularization weight consists in looking for the maximum curvature (i.e., L-curve criterion [START_REF] Hansen | The L-Curve and its Use in the Numerical Treatment of Inverse Problems[END_REF]) in the previous plot for a given deformation step (e.g., for the maximum deflection, see Figure 14(b)). According to that criterion, the regularization length 𝓁 𝑚 = 25 vx provided the best balance between Φ 𝑐 and Φ 𝑚 . This result applied to all the considered deflections and 𝓁 𝑚 = 25 vx was therefore selected (Table 2).

Figure 15 shows the gray level residual fields at the end of direct calculations for each loading step. The gray levels were the lowest close to the middle section for each scan as it did not move much (because the middle support was motionless). Higher levels were observed at the two ends of the pantograph because of very large displacements (Figures 16 and17) making the calculations more difficult. The top and bottom surface of the pantograph also concentrated higher residual levels as there was corrugated cardboard (undulating polypropylene sheet commonly known as coroplast) on the top and bottom parts of the sample to ensure smooth contact between the pantograph and the supports. Overall, the residuals remained sufficiently low for the DVC results to be deemed trustworthy.

Analysis of the results

Displacement fields

Figure 16 shows the measured vertical displacement fields on the corresponding deformed configurations. Their overall levels were similar to those measured at the macroscale (Figure 12). This observation shows that the macroscale analysis provided good initial guesses for the global characterization of the kinematics of the flexural test.

Figure 17 displays the measured longitudinal displacement fields on the corresponding deformed configurations. The displacements essentially characterize the section rotations (as was already observed on the macroscale analyses, see Figure 13), and their amplitudes increased linearly with the deflection. Very large levels were observed, which could be captured thanks to the various initialization steps that were followed herein.

In Figure 18, the measured displacement fields along the transverse direction are shown for the four deflections. Their amplitudes were significantly lower than those observed for the other two displacement components (Figure 16 and17). In some 3-point flexural tests [START_REF] Maati | Influence of the material microstructural properties on a 3-point bending test[END_REF], the emergence of Poisson effects in one direction perpendicular to that of the prescribed displacement were observed. This effect, which generally occurs in bias-extension or bias-compression tests on pantographic metamaterials as studied in this work was not detected experimentally or was almost absent, as shown in Figures 161718. As pointed out previously, a major role in the mechanics of pantographic metamaterials was played by the hinges connecting the fiber layers. The emergence of the Poisson's effect is related to the aspect ratio of the cylinders constituting the hinges. With slender hinges, as in the present case, this effect was hindered. Instead, such effect would be observed with hinges with a low aspect ratio. Further investigations are needed to clarify these questions.

Strain fields

Von Mises and maximum principal strain fields are shown in Figure 19. Both fields were equivalent in overall distributions. There was one remarkable phenomenon that was the opposite behavior of the two mutually orthogonal families of beams. Near the outer supports, the family of beams that was oriented toward the middle of the pantograph concentrated most of the strains (and thus most of the stresses). This observation is unusual for a standard flexural test to have the highest strain levels near the outer supports and not around the mid-section. Since the equivalent von Mises strain was greater than 1% in some areas at the end of the test, it was concluded that inelasticity had occurred even though the load/displacement response was still linear (Figure 5). This conclusion was confirmed by the fact that when large regularization lengths were considered, higher RMS gray level residuals were observed (Figure 14).

To further investigate this phenomenon, the maximum principal strain field is shown separately for the two families of beams in Figure 20. Apart from the first loading step, the strain distribution was no longer symmetric with respect to the mid-section plane. The deformation concentrated in the beams oriented toward the central section and located at both ends of the metamaterial.

Hinge kinematics

The analysis now focuses on the 1,935 hinges. The first set of results is related to the deformation mechanism From the previous results, it was easy to conclude directly on the deformation mode of the studied surface since it was fully dominated by the flexural mode in the longitudinal direction. Consequently, a parabolic fit was performed and the curvatures were assessed at the center of the considered surface (Figure 22). The major principal curvature was linearly related to the prescribed deflection. This is the signature of a flexural test. The minor curvature amplitude was significantly lower than that of the major curvature (as expected from Figure 21). Yet, it was negative, which indicated that the deformation mode was slightly anticlastic (i.e., negative Gaussian curvature). Interestingly, for the largest deflection step, the Gaussian curvature decreased significantly, which further indicated that nonlinear effects may have occurred. The second set of results is related to the displacement and rotation jumps for each hinge. As in the uncertainty quantification (Figure 11), the axial and transverse components of the rotation and displacement jumps are reported. Let us note that the amplitudes were at least one order of magnitude higher than their standard uncertainties. Consequently, their levels were deemed trustworthy. Figure 23(a) shows the axial displacement jump for each hinge for the last loading step. The hinges in the central part of the sample underwent rather modest displacement jumps (i.e., less than 0.5 vx in absolute value). Conversely, higher absolute levels were observed closer to the external supports.

The tangential displacement jumps were higher in the central part of the sample in comparison to their levels closer to the external supports (Figure23(b)). Similar trends were observed for the axial rotation jumps (Figure23(c)). However, their sign was alternating from one layer to the next. This was due to the fact that the two families of beams did not deform in a similar way (Figure 20). For the transverse rotation jumps (Figure23(d)), they were also alternating from one layer to the next. Their amplitudes remained smaller than those of the axial rotation jumps. The distribution of the four previous quantities are displayed in Figure 24 in terms of cumulative plots for the four deformation steps. For all four jumps, the larger the deflection, the wider their distributions. For the axial displacement jumps (Figure 24(a)), they were mostly negative (i.e., shortening of the hinges). The amplitude of the transverse displacement jumps (Figure 24(b)) was significantly higher than that in the axial direction for all deflections. The distribution of axial rotation jumps (Figure 24(c)) was virtually symmetric with respect to the 0 rad angle, which was due to the different deformation of the two families of beams (Figure 20). The amplitudes of the transverse rotation jumps (Figure 24(d)) were lower than those of the axial rotation jumps for any loading step. These observations show that the two dominant deformation modes of the hinges were shear and torsion. Conversely, axial contraction and flexure remained less active.

Very similar results were obtained on the lower row of hinges. For the sake of conciseness, they are not further discussed. Last, given the fact that the transverse displacements remained very small (Figure 18) in comparison to their vertical and horizontal components (Figures 16 and17), the two lateral surfaces of the pantographic block did not experience significant deformations.

Prediction of deformation patterns in 3-point flexure of pantographic blocks with elastic hinges

Modeling overview

The design of pantographic sheets involved two mutually orthogonal fiber arrays that were mechanically interacting thanks to continuously distributed and equidistant cylindrical joints. Their mechanical and geometric properties greatly influence the overall response of pantographic sheets. When these joints are very stiff the resulting behavior after homogenization is governed by an orthotropic first gradient strain energy [START_REF] Boutin | Linear pantographic sheets: asymptotic micro-macro models identification[END_REF]. Conversely, when they have small torsional rigidity, suitably rescaled with respect to the inter-joint distance, the response corresponds to that of a second gradient medium (dell'Isola, Seppecher, Alibert, Lekszycki, Grygoruk, Pawlikowski, Steigmann, Giorgio, Andreaus, Turco, Gołaszewski et al., 2019a;dell'Isola et al., 2019b). Elastic joints induce a microdeformation mechanism that allows for relative rotations of the fibers. In the limit case when they behave as perfect hinges, the flexural energy of the pantographic architecture becomes dominant and at the macroscale second gradient effects are observed.

Assuming that the cylindrical hinges were small and stiff enough, the two families of fibers can be described in terms of the same macrodisplacement field, as the hinges exert elastic or any constraint on their relative rotations. In macroscopic descriptions, the joints control the shear deformations and their torsional energy is equal to the macroshear energy. It was proven that small variations of their mechanical and geometric properties caused large changes in macromechanical properties, including a transition from first to second gradient macrostrain energies (Cuomo, dell'Isola, [START_REF] Cuomo | First versus second gradient energies for planar sheets with two families of inextensible fibres: investigation on deformation boundary layers, discontinuities and geometrical instabilities[END_REF]. Therefore, the models to be used in the description of pantographic sheets must be associated with second gradient 2D continua [START_REF] Dell'isola | Large deformations of planar extensible beams and pantographic lattices: heuristic homogenization, experimental and numerical examples of equilibrium[END_REF].

Pantographic blocks

The mechanical properties of the interconnecting hinges were designed to reduce all mechanical couplings that may cause deformations of cross-sections of the pantographic block. This choice was made in order to allow for the use of already available 2D second gradient continuum models. The three-dimensionality of the block was simply accounted for by using a suitable increase of 2D pantographic sheet stiffnesses (i.e., strain energy multiplied by the number of pantographic sheets composing block).

In 3-point flexure, as shown above, fiber flexure became dominant at the mesolevel. As a consequence, in the macroscopic mechanical response, the second gradient strain energy dominates. The only surviving first gradient energies are those relative to shear deformations, which are tuned by suitably choosing the geometric properties of cylindrical hinges [START_REF] Spagnuolo | Phenomenological aspects of quasi-perfect pivots in metallic pantographic structures[END_REF][START_REF] Giorgio | Numerical identification procedure between a micro-Cauchy model and a macro-second gradient model for planar pantographic structures[END_REF].

Pantographic sheet model

In the previous section, it was observed that, notwithstanding the fact that a 3D pantographic block was studied, a 2D model may be adopted because all deformations were essentially planar. The fact that multiple pantographic sheets were printed will be taken into account by suitably correcting the constitutive parameters. Let us consider a reference configuration of a pantographic lattice (Figure 25). A 2D second gradient continuum model for properly describing the pantographic sheet was introduced by dell'Isola et al. [START_REF] Dell'isola | Large deformations of planar extensible beams and pantographic lattices: heuristic homogenization, experimental and numerical examples of equilibrium[END_REF] and generalized by Spagnuolo et al. [START_REF] Spagnuolo | Qualitative pivot damage analysis in aluminum printed pantographic sheets: Numerics and experiments[END_REF]. In the following, the original model [START_REF] Dell'isola | Large deformations of planar extensible beams and pantographic lattices: heuristic homogenization, experimental and numerical examples of equilibrium[END_REF] was selected for the sake of simplicity. Let (𝐃 1 , 𝐃 2 ) be an orthogonal basis for the reference configuration, where both vectors are aligned with the fiber directions (i.e., orthogonal in the reference configuration). A 2D continuum is considered whose reference shape is given by a rectangular domain Ω = [0, L] × [0, 𝓁] ⊂ ℝ 2 , where L and 𝓁 are identified as the lengths of the sides of the continuum medium (Figure 25) and L = 3𝓁. The current shape of the rectangle is mathematically described by a regular placement function 𝝌 ∶ Ω → ℝ 2 . In Ref. [START_REF] Dell'isola | Large deformations of planar extensible beams and pantographic lattices: heuristic homogenization, experimental and numerical examples of equilibrium[END_REF], it was shown how, by means of asymptotic homogenization [START_REF] Dell'isola | Large deformations of planar extensible beams and pantographic lattices: heuristic homogenization, experimental and numerical examples of equilibrium[END_REF] and assuming that 𝝌(⋅) is at least twice differentiable, the macroscopic strain energy for the pantographic metamaterial is obtained

 (𝝌(⋅)) = ∫ Ω ∑ 𝛼 𝐾 𝑒 2 (||𝑭 D 𝛼 || -1) 2 𝑑Ω + ∫ Ω ∑ 𝛼 𝐾 𝑏 2 [ ∇ ∇ ∇𝑭 |D 𝛼 ⊗ D 𝛼 ⋅ ∇ ∇ ∇𝑭 |D 𝛼 ⊗ D 𝛼 ||𝑭 D 𝛼 || 2 - ( 𝑭 D 𝛼 ||𝑭 D 𝛼 || ⋅ ∇ ∇ ∇𝑭 |D 𝛼 ⊗ D 𝛼 ||𝑭 D 𝛼 || ) 2 ] 𝑑Ω (5)
where 𝑭 denotes the deformation gradient ∇ ∇ ∇𝝌, and

( ∇ ∇ ∇𝑭 |D 𝛼 ⊗ D 𝛼 ) 𝛽 = 𝐹 𝛽 𝛼,𝛼 = 𝜒 𝛽 ,𝛼𝛼
, with 𝛼, 𝛽 = {1, 2}. The parameters 𝐾 𝑒 and 𝐾 𝑏 represent, respectively, the elongation and bending stiffnesses. The first integral in Equation ( 5) is related to fiber elongation while the second one refers to fiber flexure. It is worth noting that the flexural energy is written in terms of the gradient of the deformation tensor ∇ ∇ ∇𝑭 , which corresponds to the second gradient of the placement function ∇ ∇ ∇ 2 𝝌. This observation is consistent with the fact that the second gradient term is representative of fiber flexure.

If microstructural joints cannot be considered as perfect hinges, as in the present study, one has to deal with their deformation adding further terms to the strain energy (Equation ( 5)). The hinges are modeled as elastic cylinders. Consequently, one should take into account their torsion and shear deformations. These two deformation mechanisms are sometimes negligible. In the present case, the diameter of the pivots was 0.90 mm and their height 2 mm. They were considered slender and their microshear should be taken into account. Yet, the experiment was designed to minimize all deformations except fiber flexure. For this reason, microshear was neglected.

A simple form for the microtorsion (corresponding to macroshear) energy term was assumed to depend on the angle between the interconnected fibers from the pivot raised to the power 𝛾 [START_REF] Dell'isola | Large deformations of planar extensible beams and pantographic lattices: heuristic homogenization, experimental and numerical examples of equilibrium[END_REF]. The latter may be obtained from the fit of experimental data, and may depend on the type of material with which the pantographic structure was additively manufactured. The proposed form for the shear energy reads

𝔘 𝑠 (𝝌(⋅))= = ∫ Ω 𝐾 𝑠 2 | | | | | cos -1 ( 𝑭 D 1 ||𝑭 D 1 || ⋅ 𝑭 D 2 ||𝑭 D 2 || ) - 𝜋 2 | | | | | 𝛾 𝑑Ω (6)
where 𝐾 𝑠 is the macro-shear stiffness, which is related to the torsional stiffness of the cylinders modeling the connecting joints. The term cos -1 (

𝑭 D 1 ||𝑭 D 1 || ⋅ 𝑭 D 2 ||𝑭 D 2 || )
-𝜋 2 represents the angular variation between two fibers in correspondence with a chosen pivot. For the numerical simulations presented hereafter, it was assumed that 𝛾 = 2.

HF2MS: info on Comsol, model implementation, etc.

Validation of modeling hypotheses

The main assumption of the bidimensional model described in the previous subsection concerned the flatness of transverse cross-sections. For a standard isotropic 3D continuum, Saint Venant solutions prove the existence of so-called anticlastic deformations of cross-sections [START_REF] Bellow | Anticlastic behavior of flat plates[END_REF]. The pantographic blocks were designed in such a way that neither anticlastic nor any other deformation would occur for transverse cross-sections. Such hypothesis is in accordance with the DVC measurements shown in Figures 2122. It was concluded that the deformation of the transverse surface of the pantographic block was weakly anticlastic. Figure 22 reports that curvature 𝜅 2 , which would vanish for synclastic deformation, was slightly negative, thereby producing a Gaussian curvature that was negative as well. This negative value was higher in modulus for the last deflection step, showing that nonlinear effects started playing a more important role. Yet, the curvature 𝜅 2 had a very low value and that this non-zero value was not due to out-of-plane deformations but to microshear of connecting hinges, which allowed for relative sliding at fiber levels that was at the basis of anticlastic deformation. It has to be noted that this effect may be considered negligible in the case under examination and for this reason the 2D model employed for numerical simulations is still valid.

A further validation of the theoretical assumptions is obtained by comparing the displacement fields displayed in Figures 16-17 with the transverse displacement in Figure 18, which was considerably lower than the other two. This result means that cross-sections did not deform much, thus the 2D model may be adopted in the description of the pantographic block.

Comparison between predicted and measured deformed shapes

Via numerical simulations based on the 2D continuum model described in Equations (5-6), the deformed shapes of a bidimensional continuum corresponding to the pantographic block were computed for the designed 3-point flexural test.

The numerical simulations were carried out using the commercial finite element code Comsol Multiphysics. This software allows a mesh to be produced in an automated manner. The simulations were conducted by implementing the equations underlying the model in weak form. Specifically, the Weak Form module was used, which allows for the minimization of the strain energy written in terms of displacement field variables and their derivatives. The finite elements used for the discretization are of the 5th order (Argyris type), while quadratic Lagrange elements were used to prescribe boundary conditions.

The constitutive parameters were calibrated using the force measurements and taking into account (i) the total number of the pantographic sheets, (ii) the fact that 3D printed materials may have mechanical properties different from the bulk material. The calibrated values for the constitutive parameters 𝐾 𝑒 , 𝐾 𝑏 and 𝐾 𝑠 are as follows 𝐾 1 𝑒 = 𝐾 2 𝑒 = 6.9 kN/m, 𝐾 1 𝑏 = 𝐾 2 𝑏 = 4.7 Nmm, 𝐾 𝑠 = 360 N/m. ( 7)

Figure 26 shows a good agreement between the experimental measurements and the predicted reaction forces once the model was calibrated. The calibration of the parameters was performed in a qualitative manner. Of course, better parameter identification systems can be implemented and in the case of discrete models this has been done. In the case of the continuum 2D model, the complexity of its implementation implies considerable effort and the creation of a code external to Comsol (e.g. in Matlab) to calibrate the parameters in order to obtain a force-displacement curve in accurate agreement with the experimental one. However, it has been observed that calibrating the parameters only on the force-displacement curve subsequently results in deformed shapes that do not coincide with the experimental ones. This suggests that the quantitative calibration technique must be specially designed and must take multiple factors into account.

In Figures 27-28, the predicted vertical and longitudinal displacement fields are shown. These fields can be compared with those displayed in Figures 1617, which were measured via DVC. They are in good qualitative agreement. It was concluded that the 2D model was consistent with macroscopic data (i.e., force-deflection curve) and with experimentally measured displacement fields. The adopted 2D model cannot be used for getting predictions analogous to the results of Figure 18, where transverse displacement fields are displayed. For the latter one, the amplitudes were low in comparison to the other two components. The observation of Figure 31 shows that the shear strain 𝜖 𝑥𝑧 was lower in magnitude than the other two components. It was thus concluded that the shear energy played a minor role in this type of test. This conclusion was also drawn by plotting the relative angle change between the two fiber families (Figure 32). As stated in Equation ( 6), the shear energy was written in terms of angular change. At the mesostructural level, the previous observation corresponds to the fact that the magnitude of joint torsion remained low. In Figure 32, the rotation fields (corresponding to relative rotation of fibers at the mesostructural level) are shown and only two portions of the specimen exhibited non negligible values. The relative angles between the two fiber families remained almost everywhere equal to 90°. Two portions of the mesostructure, instead, exhibited nonnegligible changes of angles. Specifically, there are two portions where changes in the angles are observed. These portions are located in correspondence of the vertical symmetry axis, passing through the support. Due to the fact that the concavity induced by flexure was oriented upward, then the angle between the top fibers increased, while the angle between the bottom fibers decreased. These zones were very localized and a large part of the fibers remained orthogonal even in the deformed configuration. For this reason, the mechanical behavior of the pantographic metamaterial subjected to 3-point flexure was dominated by the flexure of fibers.

Conclusion

Very few studies were conducted to quantify the deformation of pantographic metamaterials in flexure. In the present case, an in situ flexural test was performed. It allowed the metamaterial mesostructure to be imaged for different applied deflections, and kinematic details to be quantified (in particular at the beam and hinge levels) thanks to digital volume correlation (DVC). Such analyses were made possible by using FE-based DVC for which a sub-millimeter mesh, originally constructed on the nominal twin, was fitted to the underlying mesostructure thanks to a backtracking procedure. Given the large displacement amplitudes experienced by the metamaterial, a multi-step initialization was implemented in order to allow mesoscale DVC analyses to converge even under such challenging situation. Since virtually no contrast existed in the solid phase, mechanical regularization was utilized. In the present case, Hencky elasticity was assumed to penalize the DVC cost function with the equilibrium gap cost function. An L-curve criterion was followed to choose the regularization length, which led to the best balance between both cost functions. The analysis of the correlation residual fields showed that convergence was achieved with the selected parameters.

The analysis of the strain fields for the two families of beams showed that they did not deform the same way during the test. Moreover, very high levels were observed close to the external supports, which was a signature of higher order gradient continua. At the hinge level, the displacement and rotation jumps were assessed. It was concluded that the two dominant deformation mechanisms were torsion and shear in comparison to axial contraction and flexure. They led to weakly anticlastic deformations at the macroscale. As evidenced in the present study, a topic of interest for further research is related to the analysis of clasticity of pantographic sheets and blocks. Due to the fact that in this type of metamaterial large deformations and large rotations occur, one should refer to nonlinear shell theories for correctly modeling the mechanical problem [START_REF] Eremeyev | On constitutive inequalities in nonlinear theory of elastic shells[END_REF][START_REF] Eremeyev | Existence theorems in the linear theory of micropolar shells[END_REF][START_REF] Eremeyev | Basics of mechanics of micropolar shells[END_REF].

The interest of studying pantographic metamaterials lies in the fact that it requires second gradient homogenized continuum models for its description. It was shown that microstructures call for higher gradient theories in continuum descriptions [START_REF] Germain | The method of virtual power in continuum mechanics. Part 2: Microstructure[END_REF][START_REF] Germain | The method of virtual power in the mechanics of continuous media, I: Second-gradient theory[END_REF]. Specifically, a 2D planar model for pantographic sheets [START_REF] Dell'isola | Large deformations of planar extensible beams and pantographic lattices: heuristic homogenization, experimental and numerical examples of equilibrium[END_REF] was employed herein for performing numerical simulations validated by experimental measurements. This was possible for a pantographic block due to the fact that the deformations measured via DVC analyses were essentially planar. The numerical simulations were in good agreement with the experiment in terms of reaction force vs. prescribed displacement. Numerical evaluations of displacement fields were shown to be consistent with DVC results.

The study presented herein calls for further experiments. It is necessary to check the emergence of Poisson effects in relation to the aspect ratio of the hinges (i.e., clasticity). Moreover, it may be interesting to study 3-point flexural tests in the transverse direction with respect to that adopted herein. Significantly higher Poisson effects are expected to occur.

Acknowledgments

This work has been partially supported by the French "Agence Nationale de la Recherche" through the "Investissements d'avenir" program (ANR-10-EQPX-37 MATMECA Grant). MS was supported by P.O.R. Sardegna F.S.E. 2014-2020 Asse III: Istruzione e Formazione, Obiettivo Tematico: 10, Obiettivo Specifico: 10.5, Azione dell'accordo di Partenariato:10.5.12 -Avviso di chiamata per il finanziamento di Progetti di ricerca, Anno 2017 (CUP F24J17000190009).

Figure 1 :

 1 Figure 1: Example of pantographic block cross-section

Figure 2 :

 2 Figure 2: (a) 3D rendering of the studied metamaterial (length: 203.5 mm, height: 43.5 mm) with a 100 µm / vx resolution. (b) Side view of hinges (0.90 mm in diameter) and beams

Figure 3 :

 3 Figure 3: (a) Sample in the reference configuration in the in situ testing machine. (b) Sample in the deformed configuration after the first loading step

Figure 4 :

 4 Figure 4: 3D renderings of the 4 deformed configurations

Figure 5 :

 5 Figure 5: Force vs. deflection of the in situ 3-point flexural test

Figure 6 :

 6 Figure 6: Digital twins of the nominal configuration of the studied metamaterial

Figure 7 :

 7 Figure 7: (a) Mesoscale mesh of the nominal configuration with its reference frame. Side (b) and front (c) views of hinges. (d) Detail of one individual hinge ans its reference frame (axial and tangential directions). The red and blue nodes depict the top respectively the bottom planes used to measure displacement and rotation jumps

Figure 8 :

 8 Figure 8: Reference (printed) configuration to be registered with the nominal configuration via DVC using an auxiliary (i.e., macroscopic) mesh encompassing the reference configuration. The nominal mesh is then backtracked

Figure 9 :

 9 Figure 9: Displacement fields in the 𝑥 (i.e., longitudinal) (a), 𝑦 (i.e., transverse) (b) and 𝑧 (i.e., vertical) (c) directions obtained by the backtracking procedure. The fields are shown on the deformed configuration and expressed in voxels (1 vx ≡ 0.1 mm)

Figure 10 :

 10 Figure 10: Standard displacement (a) and strain (b) uncertainties as functions of the regularization length 𝓁 𝑚 . The physical size of one voxel is 0.1 mm. The cyan dotted lines show power law interpolations with exponents -3∕2 (a) and -5∕2 (b)

Figure 11 :

 11 Figure 11: Standard displacement (a) and rotation (b) jump uncertainties as functions of the regularization length 𝓁 𝑚 . The physical size of one voxel is 0.1 mm

Figure 12 :

 12 Figure 12: Initialization results in terms of displacements (expressed in voxels, 1 vx ≡ 0.1 mm) along the 𝑧 (i.e., vertical) direction for the different deformed configurations. (a) 5 mm, (b) 10 mm, (c) 15 mm and (d) 20 mm deflections

Figure 13 :

 13 Figure 13: Initialization results in terms of displacements (expressed in voxels, 1 vx ≡ 100 µm) along the 𝑥 (i.e., longitudinal) direction for the different deformed configurations. (a) 5 mm, (b) 10 mm, (c) 15 mm and (d) 20 mm deflections

Figure 14 :

 14 Figure 14: (a) Equilibrium gap Φ 𝑚 vs. gray level residual Φ 𝑐 for the five analyzed scans and for different regularization lengths (expressed in voxels, 1 vx ≡ 0.1 mm). (b) Corresponding plot for the last deformation step. The digitization of the registered volumes was 8 bits

Figure 15 :

 15 Figure 15: Gray level residuals for the different deformed configurations. (a) 5 mm, (b) 10 mm, (c) 15 mm and (d) 20 mm deflections. The digitization of the registered volumes was 8 bits

Figure 16 :Figure 17 :

 1617 Figure 16: Displacement fields along the 𝑧 (i.e., vertical) direction for the different deformed configurations. (a) 5 mm, (b) 10 mm, (c) 15 mm and (d) 20 mm deflections. The displacements are expressed in voxels (1 vx ≡ 0.1 mm)

Figure 18 :

 18 Figure 18: Displacement fields along the 𝑦 (i.e., transverse) direction for the different deformed configurations. (a) 5 mm, (b) 10 mm, (c) 15 mm and (d) 20 mm deflections. The displacements are expressed in voxels (1 vx ≡ 0.1 mm)

Figure 19 :

 19 Figure 19: Von Mises equivalent strain (left) and maximum principal strain (right) fields for the 4 different deformed configurations. (a,b) 5 mm, (c,d) 10 mm, (e,f) 15 mm and (g,h) 20 mm deflections

Figure 20 :Figure 21 :

 2021 Figure 20: Maximum principal strain fields shown for the 4 different deformed configurations on the whole mesh (left), first family of beams (center) and second family of beams (right). (a) 5 mm, (b) 10 mm, (c) 15 mm and (d) 20 mm deflections

Figure 22 :

 22 Figure 22: Principal curvatures 𝜅 1 and 𝜅 2 , corresponding mean and Gaussian curvatures for the center of the upper surface of the pantographic block (Figure 21)

Figure 23 :

 23 Figure 23: Axial (a) and transverse (b) displacement jumps (expressed in voxels, 1 vx ≡ 0.1 mm) for each hinge. Corresponding axial (c) and transverse (d) rotation jumps (expressed in rad) for each hinge. The results are shown for the last deformed configuration

Figure 24 :

 24 Figure 24: Cumulative plots of axial (a) and transverse (b) displacement jumps for the four applied deflections. The displacements are expressed in voxels (1 vx ≡ 0.1 mm). Corresponding plots for the axial (c) and transverse (d) rotation jumps

Figure 25 :

 25 Figure 25: configuration of a pantographic block

Figure 26 :

 26 Figure 26: Comparison between measured and numerically computed reaction forces vs. prescribed deflection

Figure 27 :

 27 Figure 27: Computed transverse displacement fields on deformed shapes of the pantographic block for different prescribed deflections: (a) 5 mm, (b) 10 mm, (c) 15 mm and (d) 20mm. The displacements are expressed in mm.

Figure 28 :Figure 29 :Figure 30 :

 282930 Figure 28: Computed longitudinal displacement fields on the deformed shapes of the pantographic block for different prescribed deflections: (a) 5 mm, (b) 10 mm, (c) 15 mm and (d) 20 mm. The displacements are expressed in mm.

Figure 31 :

 31 Figure 31: Computed 𝜖 𝑥𝑧 strain field on deformed shapes of the pantographic block for different prescribed displacement values: (a) 5 mm, (b) 10 mm, (c) 15 mm and (d) 20 mm

Figure 32 :

 32 Figure 32: Computed relative angles (in degrees) between the two fiber families for two deflections. (up) 10mm, (down) 20mm.

Appendix: DVC hardware and analyses parameters


Jean-Baptiste Bonnet

Guillaume Nicolet

Laurence Papinaud

Antoine Avignon

Claire Duflos

Ariane Sultan
email: a-sultan@chu-montpellier.fr

Global healthcare pathway of people living with diabetes prior to wounding is associated with a decreased risk of amputation

Keywords:

Aims:

Diabetic foot ulcer (DFU) has a poor prognosis and high amputation rate. We previously used the French National Health Data System (Syst`eme National des Donn´ees de Sant´e: SNDS) to analyze the impact of deprivation and healthcare access on DFU prognosis. The purpose of this ancillary study was to explore the relationship between the global care pathway (care consumption) the year before and after DFU and the risk of amputation.

Methods:

We conducted a study based on a cohort from the SNDS. The data came from a region of France and subjects living with a complication of DFU. We looked at care consumption one year before and one year after wound onset. Risk of amputation was calculated one year after DFU onset.

Results:

Data were extracted for 6,642 patients. Subjects with DFU had a better prognosis regarding amputation risk when they are taking cardiovascular, antibiotic, neurological, drugs. A consultation with an ophthalmologist was also linked to a better prognosis: HR = 0.71 IC95 (0.499-0.995) (p = 0.04).

Conclusion:

People with the best diabetes follow-up, even with several comorbidities, appear to have a better prognosis for their DFU. This highlights the importance of global healthcare and the care pathway in this chronic disease.

Introduction

Diabetes mellitus is complicated by diabetic foot ulcer (DFU) in one in five to one in three people living with diabetes during their lifetimes [START_REF] Armstrong | Diabetic Foot Ulcers and Their Recurrence[END_REF]. In addition to this high incidence, DFU has a poor prognosis, with a 43.8% amputation rate and a 51.7% mortality rate 6 years after hospitalization for the ulcer [START_REF] Ghanassia | Long-Term Outcome and Disability of Diabetic Patients Hospitalized for Diabetic Foot Ulcers[END_REF]. Moreover, the cost of DFU is estimated to be more than 17 billion dollars in the USA [START_REF] Boulton | The global burden of diabetic foot disease[END_REF]. Yet, few studies have considered the quality of the healthcare pathway before DFU as a predictive factor of healing.

We previously used the French National Health Data System (Système National des Données de Santé: SNDS) to analyze the impact of social deprivation and healthcare access on DFU prognosis [START_REF] Bonnet | Effects of social deprivation and healthcare access on major amputation following a diabetic foot ulcer in a French administrative area: Analysis using the French claim data[END_REF]. The SNDS database records health data for over 99% of the French population and all insurance reimbursements in the French universal social healthcare system [START_REF] Tuppin | Value of a national administrative database to guide public decisions: From the syst`eme national d'information interr´egimes de l'Assurance Maladie (SNIIRAM) to the syst`eme national des donn´ees de sant´e (SNDS) in France[END_REF]. We highlighted that four out of five DFU patients had not been hospitalized and therefore were not part of the large hospital cohorts already studied.

Although comorbidities, such as chronic kidney disease [START_REF] Bonnet | Narrative Review of the relationship between chronic kidney disease and diabetic foot ulcer[END_REF], cardiovascular disease [START_REF] Meloni | A multidisciplinary team reduces in-hospital clinical complications and mortality in patients with diabetic foot ulcers[END_REF], diabetic retinopathy [START_REF] Ramsey | Keeping an eye on the diabetic foot: The connection between diabetic eye disease and wound healing in the lower extremity[END_REF] have been described as factors of DFU poor prognosis with an increased in amputation risk, it is unknown whether management of these comorbidities could impact on the DFU prognosis. Moreover, assess glycemic status with HbA1c monitoring is also recommended [START_REF] Elsayed | Glycemic Targets: Standards of Care in Diabetes-2023[END_REF] although the impact of glycemic control on the risk of amputation is described but debated [START_REF] Wang | Diabetic foot ulcers: Classification, risk factors and management[END_REF][START_REF] Fernando | Intensive versus conventional glycaemic control for treating diabetic foot ulcers[END_REF].

Regarding healthcare professionals, we previously highlighted the link between the presence of a nurse in the living area of the patient and the reduced risk of amputation [START_REF] Bonnet | Effects of social deprivation and healthcare access on major amputation following a diabetic foot ulcer in a French administrative area: Analysis using the French claim data[END_REF]. The international working group for diabetic foot proposed a strategy for management of DFU. [START_REF] Schaper | Practical Guidelines on the prevention and management of diabetic foot disease (IWGDF 2019 update)[END_REF]. Podiatrists play a critical role in the management.

All these contacts with the healthcare system (consultation, drug consumption, biology) are referred to as "global healthcare pathway". The association between this global healthcare pathway before or after the wound and DFU prognosis is under-studied.

Therefore, the purpose of the present study was to explore the relationship between healthcare consumption (drug consumption, nurse or physician visits) the year before and the year after DFU onset and DFU prognosis in terms of amputation risk.

Methods

Study design

We conducted a study based on a retrospective cohort. The data came from a region in the south of France, eastern Occitanie (formerly Languedoc-Roussillon), which is an administrative area with a population of 2.727 million people. All the medical data that were extracted from the SNDS were anonymized. Patients were included from January 2016 to December 2017.

Population

The population was defined as subjects living with any type of diabetes with a complication of DFU. This study used an algorithm for the inclusion of people with DFU through the SNDS [START_REF] Bonnet | Effects of social deprivation and healthcare access on major amputation following a diabetic foot ulcer in a French administrative area: Analysis using the French claim data[END_REF].

Outcome

Variables

Statistical analysis

Sample size was not calculated, as the entire database of the region was analyzed. In addition, a descriptive analysis is provided. Standardized incidences were obtained using direct standardization, with the entire included population as the reference population.

Kaplan-Meier curves are displayed according to the health consumption data, and univariate hazard ratios were computed using Cox models after assessment of the linearity hypothesis and the proportional hazard hypothesis. Statistical analysis was performed using SAS 9.4.

Ethic

This study was part of the usual public health assessment of the Regional Medical Department of the Caisse Nationale d'Assurance Maladie (CNAM), the national health insurance institution in charge of the SNDS. Therefore, there was no need to obtain authorization from a Research Review Board.

Results

Population

Data were extracted on 6,642 patients. The mean age was 71.8 (SD 13.7) years, and 54.5% were male. The number of amputations was 2.03 per 100 patient-years (Table 1). Men had a higher risk of amputation (Fig. 1). The analysis showed that 97% of the population had pharmacologic treatment for diabetes, and 42% were reimbursed for insulin. Subjects with an insulin reimbursement had a worse prognosis for amputation risk: HR = 2.18 IC95 (0.50-0.99) (p = 0.047). The average duration for bandage reimbursements was 105 days (SD: 167) but the median duration was only about 59 days, which is concordant with the wide distribution of healing times (Table 1).

Impact of comorbidities on amputation risk

68% of subjects had at least two reimbursements of HbA1c dosing in the year before DFU without impact on the DFU prognosis. Regarding comorbidities treatment, a better prognosis for amputation was observed for patients taking cardiovascular drugs, antibiotics or neurological drugs. (Fig. 1).

Impact of healthcare professionals on amputation risk

Ninety-eight percent of the subjects had contact with their GP and 21% with an endocrinologist in the year before DFU occurrence (Table 1). A large proportion had a previous contact with a private-practice nurse (82.1%). Visits to these professionals were not statistically linked to the amputation rate. However, a consultation with an ophthalmologist in the year before DFU, 47% of the population, was associated with a better prognosis: HR = 0.71 IC95 (0.50-1.00) (p = 0.04).

Among subjects with a podiatrist reimbursement in the year before the wound (25.7%), the mean number of consultations was 2.47 (Fig. 1).

Reimbursement for offloading shoes was related to a worse prognosis, whether the reimbursement was before (15.5% of the population concerned) or after (20.8% of the population concerned) the beginning of the wound (Table 1).

Discussion

In this study, we highlighted that patients who were reimbursed for neurological, cardiovascular, or antibiotic drugs, or who had a consultation with an ophthalmologist in the year prior to the wound, had reduced amputation risk. We also confirmed some known factors of a poor progression of DFU such as male sex [START_REF] Sen | Meta-analysis of risk factors for amputation in diabetic foot infections[END_REF], insulin need [START_REF] Cox | Locality deprivation and Type 2 diabetes incidence: A local test of relative inequalities[END_REF], and previous DFU revealed by the recent purchase of offloading shoes [START_REF] Sen | Meta-analysis of risk factors for amputation in diabetic foot infections[END_REF].

Care for comorbidities

Receiving several drug reimbursements was related to a reduced risk of amputation, although paradoxically this finding should have revealed those patients with more comorbidities. First, we assume that not all the patients who needed these drugs were granted an appropriate prescription. This assumption was previously highlighted for the achievement of LDL targets in high-risk people living with diabetes [START_REF] Breuker | Nonachievement of LDL-cholesterol targets in patients with diabetes at very-high cardiovascular risk receiving statin treatment: Incidence and risk factors[END_REF]. Second, we can assume that good management of other diabetic complications, such as cardiovascular, might be a benefit for reducing amputation risk after DFU Assuming that all people living with diabetes, and even more so with a complication of diabetes, have a very high cardiovascular risk, the probability of not achieving the LDL-cholesterol target without treatment target is significant. This result is therefore probably in favor of reducing LDL cholesterol levels to reduce the risk of amputation. It was not within the scope of our study to determine the exposure time required to reduce this risk. This is a matter for further study. Thus, early management of these complications may prevent the development of more serious comorbidities. Therefore, comprehensive care of the person living with diabetes appears to have a positive impact on amputation risk after DFU, but further studies are needed.

Global care pathway

Nearly half our population had the annual ophthalmologist check-up and had reduced amputation risk. However, although ophthalmological complications are related to DFU [START_REF] Gazzaruso | Impact of microvascular complications on the outcomes of diabetic foot in type 2 diabetic patients with documented peripheral artery disease[END_REF], ophthalmologists are not generally involved in DFU prevention or management. Yet patients with regular follow-up might be better integrated into a healthcare pathway and thus treated earlier. To our knowledge, this is the first result of this kind for DFU prognosis. Interestingly, we did not have a similar result for a consultation with an endocrinologist (who in France are diabetologists), although only one patient in five consulted one. The GP was not tested because almost all the population had had a consultation. An evaluation of the number of GP visits and their regularity could be a next objective.

Podiatrist and offloading shoes

The few number of podiatrist consultations in a population specifically made up of patients with DFU raises questions. Efforts to systematize these consultations should help to reduce the incidence of DFU and, perhaps, a certain number of patients should not have ended up in our cohort if they had been provided with these consultations. Elsewhere, it is quite disturbing to not observe an increase in offloading shoes refunds after DFU onset when offloading is known as one of the most effective treatment for DFU. It is right that offloading shoes price is known to be a hindrance [START_REF] Bus | Guidelines on the prevention of foot ulcers in persons with diabetes (IWGDF 2019 update)[END_REF]. The link between amputation risk and an offloading shoes refund is more likely to be related to more serious DFU. Furthermore, the introduction of podiatric care has been previously linked to a decrease in the number of amputations [START_REF] Schmidt | Podiatry impact on high-low amputation ratio characteristics: A 16-year retrospective study[END_REF] and is a keystone of the International Working Group for Diabetic Foot (IWGDF) guidelines [START_REF] Schaper | Practical Guidelines on the prevention and management of diabetic foot disease (IWGDF 2019 update)[END_REF].

Strengths and limitations

The first strength of our study lies in its power and exhaustiveness. To our knowledge, this is the first broad epidemiological study in France that focuses on patients with DFU, whether or not it led to hospitalization, and healthcare consumption.

The use of the SNDS database offers the advantage of exhaustiveness, with this study having one of the largest published cohorts of subjects with DFU. It should also be noted that the lack of information about real drug consumption is compensated by exhaustive data about drug and bandage reimbursements.

The first limitation of our study is that it relies on an algorithm to define DFU that depends on the SNDS database [START_REF] Bonnet | Effects of social deprivation and healthcare access on major amputation following a diabetic foot ulcer in a French administrative area: Analysis using the French claim data[END_REF]. Furthermore, outcomes concerning the one-year amputation rate might be considered as quite restrictive clinically. A study with 3 or 5 years of follow-up should be considered, but it is important to keep in mind that the SNDS provides universal data only since 2014.

We performed only descriptive and univariate analyses. The complete assessment of the risk factors of amputation would need multivariable analyses. However, this was beyond the scope of this study.

Clinical implications and future directions

We used the SNDS database, for the first time, to explore the relationship between the global healthcare pathway, comorbidities and amputation risk after DFU in a region of France. To our knowledge, the SNDS has never been used to explore care consumption among people with DFU. We previously published our findings on DFU prognosis and healthcare access using the same methodology [START_REF] Bonnet | Effects of social deprivation and healthcare access on major amputation following a diabetic foot ulcer in a French administrative area: Analysis using the French claim data[END_REF]. We concluded that including only hospitalized patients with a DFU creates a bias, as they are the most severe cases. Thus, healthcare organizations will now have to manage these unseen DFUs in the calibration of the medical response. In addition, diabetologists also need to be informed that they do not see four out of five patients with DFU, whereas these patients have almost all seen a GP. Therefore, GPs and nurses should be included in the elaboration of new healthcare pathways for the management of DFU. Moreover, they need to know about the nearest DFU reference centers, and dialogue between all these professionals should be encouraged.

At last, the management of co-morbidities appears to be effective in reducing the risk of amputation after DFU. Preventing cardiovascular risk by consuming drugs targeted against this risk is, for instance, a good example of our study.

Conclusion

We highlight the link between the global healthcare pathway and the risk of amputation after DFU. People with the best diabetes follow-up, even with several comorbidities, appear to have a better prognosis of their ulcer. This highlights the importance of the care pathway in this chronic disease. These results can, inspire health policies to promote access to care for people living with diabetes. Reducing the risk of amputation secondary to diabetes requires a global health effort.

Data availability

Acknowledgments

The authors would like to acknowledge the support of Bernard Ledesert from the Centre Régional d'Études, d'Actions et d'Information, Observatoire Régional de la Santé for assistance in identifying relevant datasets of social deprivation.

Fun**ding

This research received no specific grant from any funding agency in the public, commercial or not-for-profit sectors.

8.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.