Harvesting Deformation Modes for Micromorphic Homogenization from Experiments on Mechanical Metamaterials

A micromorphic computational homogenization framework has recently been developed to deal with materials showing long-range correlated interactions, i.e. displaying patterning modes. Typical examples of such materials are elastomeric mechanical metamaterials, in which patterning emerges from local buckling of the underlying microstructure. Because pattern transformations significantly influence the resulting effective behaviour, it is vital to distinguish them from the overall deformation. To this end, the following kinematic decomposition into three parts was introduced in the micromorphic scheme: (i) a smooth mean displacement field, corresponding to the slowly varying deformation at the macro-scale, (ii) a long-range correlated fluctuation field, related to the buckling pattern at the meso-scale, and (iii) the remaining uncorrelated local microfluctuation field at the micro-scale. The micromorphic framework has proven to be capable of predicting relevant mechanical behaviour, including size effects and spatial as well as temporal mixing of patterns in elastomeric metamaterials, making it a powerful tool to design metamaterials for engineering applications. The long-range correlated fluctuation fields need to be, however, provided a priori as input parameters. The main goal of this study is experimental identification of the decomposed kinematics in cellular metamaterials based on the three-part ansatz. To this end, a full-field micromorphic Integrated Digital Image Correlation (IDIC) technique has been developed. The methodology is formulated for finite-size cellular elastomeric metamaterial specimens deformed in (i) virtually generated images and (ii) experimental images attained during in-situ compression of specimens with millimetre sized microstructure using optical microscopy.


Introduction
Mechanical metamaterials have recently gained substantial attention in the literature due to their potential for various engineering applications such as strain tunable photonic crystals (Krishnan and Johnson, 2009), tunable acoustic properties (Zhao et al., 2015;Bilal et al., 2017;Wormser et al., 2017;Lin et al., 2023), soft robotics (Yang et al., 2015;Mirzaali et al., 2018b;Yu et al., 2018;Shen et al., 2023), programmable mechanical metamaterials (Florijn et al., 2014;Bertoldi et al., 2017;Kolken and Zadpoor, 2017;Li et al., 2023), biomedical prosthetics (Yavari et al., 2015;Zadpoor, 2017), or functionally graded metamaterials (Francesconi et al., 2019;Farzaneh et al., 2022).Various aspects of such metamaterials with 2D (Bertoldi et al., 2008) and 3D (Javid et al., 2016;Hedayati et al., 2019) microstructures have been investigated.A number of studies focused on the effect of the hole geometry and stacking (Bertoldi et al., 2010;Overvelde et al., 2012;Overvelde and Bertoldi, 2014;Hu et al., 2013).Several studies reported in the literature made use of full-field displacement measurement techniques to investigate cellular metamaterials, among which Slann et al. (2015) used local Digital Image Correlation (DIC) to validate their numerical study on auxetic metamaterials, Xu and Pasini (2016) used DIC in the design of thermally expanding tunable 3D metamaterials, Easey et al. (2019) studied different designs of auxetic dome-shaped structures using stereo DIC, and Tang et al. (2015) compared DIC and finite element analysis results on highly stretchable and reconfigurable metamaterials.Apart from displacement and strain assessment, Mullin et al. (2007) made a qualitative comparison between numerically and experimentally assessed stresses in cellular elastomers.Besides the geometrical and morphological effects, also the type of material is important.Mirzaali et al. (2018a) used 3D printing to manufacture and study multi-material mechanical metamaterials, Wu et al. (2014) investigated different pattern transformations as a result of swelling in cellular hydrogel membranes, whereas Hu et al. (2018) investigated the effect of friction on buckling pattern formation in cellular structures.A large number of works also focused on experimental investigation of pantographic metamaterials, e.g., (Auger et al., 2020;Hild et al., 2020;Barchiesi et al., 2021).
The common denominator to the listed applications are exotic properties of metamaterials, typically emerging from a coordinated motion of their microstructure, i.e., from the so-called patterning.The patterning usually occurs upon reaching a critical level of compressive load, resulting in a local buckling of the underlying microstructure and a significant change in the overall mechanical properties.Fig. 1a shows an example of a cellular elastomeric metamaterial specimen with a regular grid of circular holes, where the local buckling results in a pattern of alternating elliptical holes (Fig. 1b).Here it can be observed that the patterning is locally constrained by external boundaries, and depends on the size of the holes relative to the size of the entire specimen (i.e., on the scale ratio), studied computationally in detail by Ameen et al. (2018).Since such irregularities in the pattern directly influence the overall behaviour of the specimen-even on a larger scale-, it is important from the engineering perspective to be able to identify these distributions, independently of the global deformation, for a proper assessment of and design for certain mechanical properties.
Apart from classical homogenization approaches mostly limited to linear elastic behaviour, such as, e.g., (Toupin, 1962;Mindlin, 1964Mindlin, , 1965)), an effective tool for predicting mechanical behaviour of materials with complex microstructures on an engineering scale and operating in a non-linear regime is first-order computational homogenization (Kouznetsova et al., 2001;Miehe and Koch, 2002).To determine effective constitutive behaviour, computational homogenization solves a microscopic boundary value problem associated with each macroscopic integration point.Obtained results are used to compute effective stresses, which are equilibrated at the macroscale.Because separation of scales is assumed, such a scheme is able to properly model only macroscopically local continuum.For better approximation, extensions to second-order homogenization schemes were proposed (Kouznetsova et al., 2004;dell'Isola et al., 2019;Giorgio et al., 2023), which use gradient of the deformation gradient to account for non-local effects.Non-locality can also be addressed using generalized micromorphic theories, e.g., (Forest, 2009;Hütter, 2017;Biswas and Poh, 2017;Rokoš et al., 2019;Rokoš et al., 2020), in which additional kinematic fields are introduced to account for the underlying microfluctuations.For a recent overview and further approaches see, e.g., Findeisen et al. (2020).
Of particular interest in this contribution is the micromorphic computational homogenization scheme based on recent theoretical and numerical studies (Rokoš et al., 2019;Rokoš et al., 2020;van Bree et al., 2020), which relies on a general decomposition of the displace-ment field according to the following ansatz (1) Here, the total displacement u(x), which is a vector field expressed as a function of the material coordinate in the undeformed configuration x, is divided into three parts: (i) the mean smooth displacement field, v 0 (x), corresponding to the slow variations at the macro-scale; (ii) long-range correlated fluctuations at the meso-scale, represented by the vector fields φ i (x), i = 1, . . ., n, scaled by their spatial magnitudes v i (x); and (iii) a remaining microfluctuation field, w(x), representing the uncorrelated fluctuations at the micro-scale.By means of the patterning modes, φ i , the micromorhpic scheme is capable of capturing spatial nonlocal effects in the microstructure, present typically because of the lack of separation of scales and emerging due to effects such as boundary layers or temporal and spatial mixing of modes.Non-local effects have been studied numerically as well as experimentally within the context of metamaterial structures in, e.g., (Ameen et al., 2018;Maraghechi et al., 2020).Using a variational formulation towards a numerical approximation, (cf.van Bree et al., 2020), an effective potential energy is minimized with respect to the kinematic decomposition of Eq. ( 1), resulting in a micromorphic type of effective continuum, i.e., not of the standard Cauchy-Born type.Although the standard balance of momentum equation ∇ 0 • P = 0 governs the evolution of the v 0 (x) field, where P is the homogenized first Piola-Kirchhoff stress tensor energetically conjugated with the deformation gradient tensor F = I+(∇ 0 v 0 (x)) T , an additional coupled micromorphic balance equation emerges.Herein, ∇ 0 denotes the gradient operator with respect to the reference configuration.This equation governs the evolution of v i (x), which is in the form of ∇ 0 • Λ − Π i = 0, where Λ is energetically conjugated with ∇ 0 v i (x) while Π i is conjugated with v i (x).Since v i (x) corresponds to the magnitude of the patterning fluctuation field φ i (x), its kinematic interpretation is straightforward and essential boundary conditions can be identified as v i (x) = 0 along any fixed boundaries.On the contrary, Λ • N = 0 on the free boundaries may be adopted (where N denotes the outer normal).We refer the interested reader for further details and derivations to (Rokoš et al., 2020, Section 3.3).Note that unlike second-order computational homogenization (e.g., Kouznetsova et al., 2004), which uses continuity of the gradient of the deformation gradient tensor field F to account for non-locality and size effects, micromorphic computational homogenization introduces an additional independent kinematical field with a governing equation complementing the standard balance of linear momentum.The homogenization scheme can thus facilitate optimized designs of engineering systems, including deeper understanding of their behaviour and development of new metamaterials.Although the kinematic fields v 0 , v i and w are an outcome of a numerical analysis, the microfluctuation patterning fields φ i are input parameters characteristic to the geometry of the underlying microstructure that need to be provided a priori, in analogy to material constants of an underlying constitutive model.The patterning modes can be estimated through an educated guess (cf.Rokoš et al., 2019), or a numerical Bloch-type analysis (cf.Bertoldi et al., 2008).It would, however, be much better to measure the patterning modes directly from experiments, as this would include the real material behaviour, all micro-fabrication flaws, and errors in the correction of spatial distortions in the geometrical images.Therefore, the goal of this contribution is to develop a novel full-field Integrated Digital Image Correlation (IDIC) framework for quantitative experimental identification of the kinematic fields decomposed according to Eq. ( 1), including the meso-scale long-range correlated fluctuation fields φ i for the micromorphic homogenization scheme.To this end, highresolution images captured at different stages of deformation during an experimental test of a finite-sized specimen are analysed with this micromorphic IDIC framework.IDIC is a powerful full-field displacement measurement technique based on regularization of the kinematics according to specific knowledge about the problem at hand, resulting in identification of the optimal input model parameters minimizing the difference between experimentally measured data and model predictions (Roux and Hild, 2006;Neggers et al., 2015;Blaysat et al., 2015).Knowledge of the nature of the problem at hand can be integrated in the DIC scheme in the form of a parametric analytical description, which is utilized here in the form of the kinematic fields of interest.Specifically, the total displacement field is regularized based on the ansatz of Eq. ( 1) by parametrizing all relevant kinematic components, providing upon correlation the smooth mean field, v 0 , the long-range correlated fluctuation modes carried by the microstructure, φ i , as well as their spatial distribution magnitudes, v i .
The proposed methodology is tested on virtually generated and deformed images, as well as optical microscopy images taken during an in-situ compression test on a finite-size cellular elastomeric metamaterial with rectangular stacking of circular holes, cf.Fig. 1.A spectral density analysis is shown to be an effective tool for estimation, parametrization, and initialization of the fluctuation mode for a specific metamaterial.It is demonstrated that the micromorphic IDIC scheme leads to a proper decomposition of the kinematic fields of interest, both before and after the emergence of the fluctuation pattern.Considering geometrical arguments, the patterning mode is shown to be independent of the considered periodic cell size, although effective mechanical behaviour of such a material is still size dependent.Moreover, it is furthermore shown that the patterning mode is only weakly dependent on the hole diameter versus unit cell size ratio, which enables the design of metamaterials with graded microstructures.
Although a square stacking of holes is used throughout to guide and validate our developments, the proposed methodology is applicable for a general class of cellular metamaterials including hexagonal stacking of holes and chiral or re-entrant microstructures, as well as graded microstructures in terms of geometry (stacking of holes) and material properties (e.g., stiffness).The reason for focusing on the square stacking of holes herein is twofold.First, our detailed analysis requires testing inside an in-situ micro-tensile stage underneath an optical microscope, which means that the total specimen size needs to fit within a small, rectangular space.Therefore, square stacking of holes is preferred as it results in a thin, constant-thickness boundary layer along the edge, in contrast to other microstructures that may result in a boundary layer with a varying thickness along the edge.Second, a rectangular specimen shape with square stacking of holes is easiest to manufacture with high quality and apply a constant load along the boundary edge.Finally, the remaining boundary layer effect of square-holed microstructure has been studied in detail (for precisely the same rectangular specimen shape) elsewhere (Ameen et al., 2018;Maraghechi et al., 2020), so this configuration is known to work well.But, as indicated earlier, it will be demonstrated that the herein proposed methodology is generally applicable also to various (graded) geometries.

Integrated Digital Image Correlation
Integrated digital image correlation is a full-field displacement measurement method based on the minimization of the residual field, which is the difference of the reference image, f , and the deformed image, g, probed at, respectively, reference coordinates and the deformed positions corresponding to those coordinates based on the displacement field.The problem is solved over the whole region of interest, where the residual field is defined as with ϕ M (x, a) = x + u (x, a) being the mapping function corresponding to the mechanical displacement field, u.The displacement field, in turn, is based on the regularization according to a model with specific knowledge of the kinematics at hand, thereby reducing the number of degrees of freedom (dof) contained in a, and improving the robustness of the solution (Ruybalid et al., 2016).This model can be a numerical model in which the dofs would be the material parameters or an analytical model describing the deformation field.
The optimal values of the degrees of freedom, a opt , are found as This minimization problem is then solved by an iterative scheme such as Gauss-Newton algorithm (Neggers et al., 2016).

Regularization of the Kinematics Based on the Micromorphic Kinematical Ansatz
The kinematical ansatz introduced in Eq. ( 1) is used to experimentally obtain the decomposed kinematic fields of cellular metamaterials.The local microfluctuations are the uncorrelated part of the kinematics and by definition unknown, therefore w(x) is not included in the micromorphic IDIC parametrization for cellular metamaterials (in fact, w(x) is directly related to the residual field r(x, a) in Eq. ( 2) via the image gradient).Accordingly, the displacement field for the IDIC problem is defined as where a = a v 0 a v i a φ i T is the column of degrees of freedom, with a v 0 , a v i and a φ i the sets of dofs for v 0 , v i and φ i .Since each mode φ i is periodic within an integer multiple of the unit cell size (Rokoš et al., 2019), a truncated 2D Fourier series is a natural choice for its parametrization.As will be shown below, the selection of relevant sine waves can be easily done by means of a spectral density analysis of one example for any family of cellular metamaterials and loading conditions.Note that the number of modes that are activated by different loading conditions in cellular metamaterials is often only one, e.g., Fig. 1, and always limited to a few (Rokoš et al., 2020).Each long-range correlated fluctuation mode φ i is defined such that it is non-zero where the cellular microstructure exists and decays linearly to zero inside a region of surrounding (non-cellular) bulk material over a width of half the unit cell size.The boundary between bulk and cellular material is a priori selected manually by visual inspection of the reference image of the undeformed geometry; high accuracy in boundary identification is not needed.Because v 0 and v i are slow macroscopic fields, smooth functions such as globally-supported polynomials are the best choice for their parametrization.To this end, Chebyshev polynomials of 5th and 6th order will be used hereafter to parametrize v 0 and v 1 .Since the behaviour of these polynomials outside the specimen's domain is irrelevant (because no kinematics is considered there), potentially large oscillations may be ignored in these regions.
Finally, in any in-situ test under optical visualization, in addition to the mechanical displacement field, the images are affected by optical distortions of the optical lenses.Such distortions may introduce large errors if neglected.Therefore, the residual field of Eq. ( 2) is rewritten as where ϕ S (x) = x+S(x) is the mapping function describing the spatial distortion (Maraghechi et al., 2018).The spatial distortion field, S(x), is determined a priori in a calibration step by the method introduced in (Maraghechi et al., 2019).Hence, there are no additional degrees of freedom related to spatial distortion that need to be identified in the micromorphic IDIC problem of Eq. (3).

Results and Discussion
To validate the micromorphic IDIC framework, it is applied to both virtual (Section 3.2) and real in-situ micro-compression experiments (Section 3.3) on cellular elastomeric metamaterial.Its general applicability to various types of metamaterial microstructures in terms of spectral analysis is first demonstrated in Section 3.1 for square and hexagonal stacking of holes, as well as for a chiral metamaterial, whereas considerations towards graded microstructures are discussed in Section 3.4.

Parametrization and Initial Guess of the Patterning Modes φ i
In the underlying deformed or buckled microstructure, corresponding to an infinitely large specimen or a periodic arrangement of unit cells, typically a long-range correlated pattern emerges, implying non-local behaviour.Such an order is a periodically repeating pattern which can well be represented with a few lowest-frequency terms of a truncated Fourier series.To show this, we perform a spectral density analysis (Grigoriu, 2002) on the reference displacement field, available from numerical simulations obtained for infinite periodic specimens and for which the affine part of the deformation has been subtracted.Note that in real experiments the displacement fields will also be available from regular DIC, as shown below in Section 3.3.Three typical metamaterial microstructures of square and hexagonal stacking of holes and a chiral metamaterial are considered.All three geometries are depicted in the first column of Fig. 2, their corresponding deformed states for various types of compression (uniaxial or biaxial) in the second column, whereas spectral densities in the last two columns.The spectral densities are the energy spectral densities of the displacement components in the horizontal and vertical directions, S xx = ûx (ξ)û * x (ξ) and S yy = ûy (ξ)û * y (ξ), where u x and u y are the x-and y-components of u(x), the hat indicates the Fourier transformation, and * complex conjugate.Inspection of spectral density graphs indicates that most of the energy is localized only around a few particular wave vectors, hence allowing to make a specific choice on the parametrization and initialization of the longrange correlated fluctuation field, see, e.g., Rokoš et al. (2019) and Eq. ( 6) below, and this parametrization plus initialization is then applicable to a broad range of metamaterials.Note that in the spectral density graphs there are multiple secondary peaks present, corresponding to the microfluctuation w, which are (almost) invisible in the images by the naked eye due to their small magnitudes.

Virtual Experiments
The micromorphic IDIC scheme will be first tested on virtually generated and deformed speckle images.By this means the method is tested in a case where there is no influence of image distortions and noise in the images, and there are no imperfections in the specimen and the loading, etc.The geometry of the specimen analysed virtually and later experimentally (although for a different scale ratio) is depicted in Fig. 3.It has hole diameters d = 1.5 mm and centre to centre pitch l = 1.9 mm (i.e., unit cell size l) with edges of 1.9 mm and 0.9 mm bulk material in loading and transverse directions.The length of the specimen, excluding the bulk edges, is denoted L. The scale ratio of the specimen, L/l, is equal to the number of holes in the loading direction.Details on the specimen fabrication are discussed thoroughly in the PhD thesis of Maraghechi (2019).To obtain virtual data, a plane strain finite element simulation, using quadratic isoparametric triangular elements with three-point Gauss integration rule and large-strain Total Lagrangian formulation with stability control and bifurcation, is performed for the compression of a cellular elastomeric metamaterial specimen with scale ratio of 6, i.e., 6 and 10 holes in loading and transverse directions, respectively.Dirichlet boundary conditions are applied on the two vertical edges perpendicular to the loading direction to apply overall compression to the specimen, while the two remaining edges are free.Upon reaching the critical overall strain, the imposed deformation triggers an instability in the underlying microstructure leading to patterning.The material, i.e., PDMS, is modelled by a hyper-elastic Ogden material model.The material considered, although rubber-like, is of sufficiently high compressibility (i.e., the initial Poisson's ratio < 0.5) to avoid any numerical problems with locking of elements.The compressibility of the material has been verified numerically as well as experimentally.More details on the numerical simulations are presented in PhD thesis (Maraghechi, 2019).The displacement fields obtained from this simulation are used to deform a speckle pattern for Loading Direction l (1.9 mm) L 1.9 mm 0.9 mm 1.5 mm φ Figure 3: Analysed cellular elastomeric metamaterial specimen, depicting the total length of specimen excluding the bulk side edges (L: the large dashed square), the size of the unit cells (l = 1.9 mm: small dashed square), hole diameter (d = 1.5 mm), the size of the bulk side edges and the loading direction.
two different time steps, one before microstructural buckling and the other after the emergence of the patterned fluctuation field, corresponding to 7.8% and 12.8% applied nominal strain.A total formulation within the context of IDIC was employed, in which two snapshots in a deformed state (pre-bifurcation and post-bifurcation images) have been correlated with the reference undeformed image, although an updated formulation correlating incremental states is equally possible as well.To generate deformed images, pixel intensities of the reference image are mapped by displacements, interpolated at nodal or pixel positions, by inverting the iso-parametric mappings of the underlying FE approximations.This step is performed to avoid any potential additional biases due to interpolation errors.Fig. 4a depicts the deformed configuration of the virtual speckle pattern after the buckling-induced pattern emerged, which exhibits the pattern transformation more pronounced in the centre and restricted at the edges, specifically the vertical edges.The employed speckle pattern is obtained as a thresholded two-dimensional scalar and stationary random Gaussian field, generated according to (Grigoriu, 2002, Section 5.3.1.2).The spectral density is chosen as a radially-symmetric normal probability density function with the mean zero and standard deviation 1/20.No additive image white noise has been added to such generated images.The pattern is blurred using Gaussian filter with the standard deviation of 3 pixels, resulting in the pattern shown in Fig. 4a.
As the virtual images are free of spatial distortions, ϕ S (x) = x is set for the virtual experiments.Chebyshev polynomials of 5th and 6th order are used to parametrize v 0 and v 1 , respectively.In order to attain a proper initial guess for φ i (x, a φ i ), it is first realized that only one long-range correlated fluctuation mode is triggered in the specimen, as can be easily seen in Fig. 4a, thus n = 1 is set in Eq. (1) with only φ 1 considered.It may also be clear that the triggered mode is periodic with a periodicity of two unit cells, recall Figs.2c-d where the peaks are localized at ξ 1 = ξ 2 = ±1/2.The horizontal and vertical lines in Figs.4b and  4c (not present in Figs.2c-d) correspond to the discontinuities in the mean deformation due to the edges of the specimen, while the four peaks correspond to the correlated fluctuation mode, revealing two sine functions roughly in diagonal directions of the Cartesian coordinate system for the x and y components of φ 1 .The purpose of the spectral density analysis is merely the rough identification of the type of long-range correlated fluctuation modes, so only the corresponding four peaks are considered here.Inspection of ûx (ξ) and ûy (ξ), at the frequencies corresponding to the four peaks of S xx and S yy , suggests that the two sine waves are of the same sign in the x component, and of the opposite sign in the y component of the displacement field.This inspection also suggests that the two sine waves are of the same amplitude, for both x and y components of the displacement in agreement with the symmetry of the problem.Based on these observations, φ 1 is parametrized as the sum of two sine functions for each component, i.e., where a φ 1 = [a 1 , a 2 , a 3 , a 4 , a 5 , a 6 , a 7 ] T is the column of degrees of freedom describing the mode.Note that the direction of both sine waves and their wavelengths are free and independent, which compensates for inaccuracies in the measurement of the unit cell size due to processing uncertainties.This allows for an accurate identification of the shape and magnitude of the fluctuation mode in the micromorphic IDIC analysis, whereas the spectral density analysis is only used for a proper parametrization of the fluctuation mode and for an approximate initial guess of all these parameters.Note also that it is sufficient to perform the spectral density analysis only once for a family of cellular metamaterials triggering the same pattern.In cases where more than one mode is activated, the spectral analysis should be done separately on each of the subregions containing the individual modes, cf.Section 3.1.
The dofs describing φ 1 are initialized as: a φ 1 = [1, 1, 1, 1, 1, 0, 0] T .The initial values for a 2 to a 5 are set to a wavelength of 2 unit cells and their signs and ratios are set such that the initial orientation of the two sine waves is aligned with the two diagonals of the coordinate system.Note that a 1 is only defining the ratio of the amplitude of the mode in x and y direction, while v 1 (x) defines its absolute amplitude according to Eq. ( 4).In order to assure a proper scaling of the minimization problem, φ 1 (x) is scaled such that the initial maximum value of the mode is 0.04.Normalization of the microfluctuation field φ i is an inherent choice of the employed micromorphic computational homogenization formulation.We refer to (Rokoš et al., 2020, Eq. ( 3)), for more details.The scaling max x ∥φ i (x)∥ = 0.04 is chosen to avoid poor conditioning of the resulting Hessian matrix during the solution procedure using the Gauss-Newton algorithm.Moreover, since the sensitivity functions of the mode φ 1 used in the Gauss-Newton optimization algorithm (Neggers et al., 2016) are a function of v 1 (x), the initial guess for the latter needs to be non-zero to avoid an ill-posed problem in the first iteration.To this end, coefficients of linear combination of the Chebyshev polynomials corresponding to the v 1 (x) field are initiated such that a small constant non-zero spatial field results.The applied compressive global strain results in large displacement on the edges of the specimen requiring an approximate initial guess in order to ensure convergence.Thus, the smooth mean field v 0 (x) is initialized such that the first order term in x direction approximately accounts for the applied global strain.The IDIC procedure overall identifies 77 dofs, out of which 2•21 = 42 correspond to coefficients of linear combination a v 0 of the Chebyshev polynomials of 5th order describing the field v 0 (x), 28 to coefficients of linear combination a v 1 of the Chebyshev polynomials of 6th order describing the field v 1 (x), and 7 corresponding to the degrees of freedom a φ 1 of the long-range correlated fluctuation field φ 1 (x) specified in Eq. ( 6).This is a negligible amount of dofs compared to local DIC which needs at least tens of thousands dofs to capture such locally-varying displacement fields.
Fig. 5 shows the results of the virtual experiments before (left) and after (right) the emergence of the buckling pattern.The micromorphic IDIC results in an accurate identification of the v 0 (x), v 1 (x) and φ 1 (x), considering the residual fields that are small everywhere both before and after buckling occurs, see Figs. 5g and 5o.For sensible initial guesses the convergence is consistent, but the procedure may be trapped in local minima during minimization for poor initial guesses.For each case, the mean smooth displacement field, v 0 (x), in x and y directions, and v 1 (x) are depicted in the first two rows, i.e., Figs.5a-c before and Figs.5i-k after buckling occurs.The v 0 fields show compression in x direction both before and after buckling (a linear profile).Before pattern transformation, lateral expansion is visible in the y direction, while afterwards the auxetic effect, i.e., lateral contraction, is clearly visible near the specimen's centre.As expected, the v 1 field is close to zero before buckling occurs, although its absolute value increases slightly in the corners of the region with cellular microstructure (the area of the specimen excluding the bulk edges), which is due to the fact that the fluctuation mode, φ 1 , partially captures the deformation that occurs in the corner holes where the edge effect is considerable.After the buckling occurs, the v 1 field is close to zero on the edges and increases towards the centre, as expected.Again, v 1 is not exactly After microstructural buckling zero at the edges and the corners of the metamaterial, since the fluctuation mode partially describes the deformation in the corners where the edge effect is considerable.Note that the large values in v 1 (x) are due to the scaling of the corresponding mode φ 1 (maximum value of 0.04).
The identified long-range correlated mode, φ 1 (x), is depicted in Figs.5d and 5l before and after buckling, where their amplitudes are scaled to assure the visibility of the deformed shape of the modes.The small v 1 prior to buckling entails a low sensitivity of the mode parameters.As a result, φ 1 is significantly perturbed with respect to the initial guess, see Fig. 5d where the periodicity does not match twice the unit cell size any more.After buckling, however, the mode becomes representative, and it resembles the initial guess more closely and matches the mode observed in the simulations (Fig. 5l).The long-range correlated fluctuation fields, v 1 φ 1 (x) in x and y directions, are depicted in their deformed configurations prior to buckling in Figs.5e and 5f, and beyond buckling in Figs.5m and 5n.Prior to buckling, v 1 φ 1 (x) is close to zero everywhere apart from the four corners, where the identified mode, although significantly perturbed with respect to the initial guess, partially captures the edge effect described above.Post buckling, φ 1 is identified properly, resulting in small modifications to the wavelengths and orientations of the sine waves compared to the initial guess based on spectral density analysis.
The residual fields in normalized grey scales before and after buckling are shown in Figs.5g and 5o, in the deformed configuration based on the total evaluated displacement u(x).The residual fields are consistently small everywhere over the region of interest.However, the microfluctuation fields, discarded in the micromorphic IDIC scheme which seeks to extract the correlated fields only, are by definition reflected in the residual fields.In order to evaluate the performance of the method, the microfluctuation field w(x) is recovered as the difference of the reference displacement field (here available from the simulation) and the displacement field u(x) according to Eq. ( 4).The Euclidean norm of the microfluctuation field, ||w||, is depicted in Figs.5h and 5p before and after the buckling, respectively.Note that the micromorphic IDIC scheme is constrained to identify displacement fluctuations with spatial frequencies bound to the highest spatial frequencies in the long-range correlated fluctuation modes, φ i .Thus, the microfluctuation field, w, contains all parts of the field with spatial frequencies higher than the characteristic frequency of φ 1 , as well visible far from the edges of the specimen in the post buckling case, see Fig. 5p.As expected, the microfluctuation field, w, is of a significantly smaller amplitude than the correlated fluctuation field beyond buckling (compare Fig. 5p with Figs.5m and 5n).The opposite is true prior to buckling, due to the almost negligible contribution of the long-range correlated fluctuation field (compare Fig. 5h with Figs.5e and 5f).Note the difference in scale bars for the figures of the pre-and post-bifurcation state.
In order to assess the robustness of the method, the correlations of the virtual experiments, after the onset of buckling, are repeated with different initial guess values for the fluctuation mode.The values of a φ i are perturbed by a random value up to 10% of the nominal value stated above, resulting in random changes in the wavelengths and orientations of the sine waves as well as the ratio of the x and y components of φ 1 .The perturbed correla-tions are repeated 100 times, and the error in each case is evaluated as: ϵ = 100 × ||u||−||upert|| ||u|| , where u is the total displacement field evaluated with no perturbation of initial guess, and u pert is the perturbed total displacement field.The error averaged over all the repetitions is 0.86%, entailing only two cases with an erroneous assessment of the total displacement field, i.e., ϵ > 1.5%, confirming the 98% robustness of the method to initial guess.The average error ϵ, ignoring the two outlier cases, is 0.45%, confirming also the high accuracy of the methodology.

Real Experiments
The results of an in-situ compression test on a cellular elastomeric metamaterial specimen, manufactured using a customized mould as presented in PhD thesis of Maraghechi (2019), are analysed with the micromorphic IDIC method.The real experiment includes a range of additional complexities compared to the virtual experiment, including the image distortions and image capturing noise as well as the imperfections in the specimens and loading conditions of the test.The specimen is 22.5 mm thick and has 12 and 10 holes in loading and transverse directions, corresponding to a scale ratio 12.The specimen is made of Sylgard 184, with a 10:1 mixing ratio.A two-layer speckle pattern is deposited on the top surface (white powder sprayed to make a fine-grained background and India ink sprayed using an air brush, see Fig. 6a), resulting in a speckle size between 30 and 80 µm (3 to 8 pixels).An in-situ compression test is performed with a Kammrath & Weiss micro tensile/compression stage with a 50 N load cell, while a Zeiss V20 stereo microscope is used to acquire images at each load step at 7.6× magnification and 2751 × 2207 pixels.The stereo microscope and the objective was used for proper imaging of the specimen providing highenough resolution for the fine speckle pattern applied on the specimen with an adequate level of detail.More details on the specimen processing and the in-situ test can be found in PhD thesis of Maraghechi (2019).Two images obtained during the in-situ test, before and after the emergence of the buckling pattern, corresponding to 8.2% and 12.7% applied nominal strain on the metamaterial, are correlated to the undeformed reference configuration.The deformed configuration after emergence of the buckling pattern is shown in Fig. 6a.
A spatial distortion calibration step Maraghechi et al. (2018Maraghechi et al. ( , 2019)), is performed to measure the spatial distortions of the imaging system at the desired magnification.The resulting distortion field, which is as large as 10% of the actual displacements, is used to construct the distortion mapping function S(x).
Chebyshev polynomials of 5th and 6th order are used to parametrize v 0 and v 1 .The same procedure as in the case of virtual experiments is performed for estimating the reduced regularization of φ 1 (x) and the initialization of the associated parameters.To this end, first, the post buckling kinematics is determined using local DIC, in order to attain S xx and S yy , shown in Figs.6b and 6c.Similar to the virtual experiments, vertical and horizontal lines relate to discontinuities in the smooth mean deformations due to the specimen's edges, while the two sine waves in approximately diagonal directions for each displacement component correspond to a single long-range fluctuation mode.Based on these observations, the fluctuation mode is assigned the same way as for the virtual experiments, i.e., Eq. ( 6) with a φ 1 = [1, 1, 1, 1, 1, 0, 0] T as initial guess.In order to assure a proper scaling of the minimization problem, φ 1 (x) is scaled such that the initial maximum value of the mode is 0.04.In order to avoid an ill-posed problem in the first iteration, v 1 (x) is again initiated with a small constant value in space.The applied compressive global strain results in large displacement on the edges of the specimen, which requires an approximate initial guess in order to ensure convergence.Thus, the smooth mean field v 0 (x) is initialized such that the first order term in x direction approximately accounts for the applied global strain.Similarly to the virtual experiment of Section 3.2, the IDIC procedure identifies in total 77 dofs a.
For these real in-situ compressive experiments, the micromorphic IDIC also results in a proper and accurate identification of v 0 (x), v 1 (x) and φ 1 (x), whereby the residual fields are small everywhere both before and after buckling, see Figs. 7g and 7o.Convergence is consistent.Fig. 7 shows the results of applying micromorphic IDIC on real experiments before (left) and after (right) the emergence of the buckling pattern.For each case, the mean smooth displacement field, v 0 (x), in x and y directions, and v 1 (x) are depicted in the first two rows, i.e., Figs.7a-c before and Figs.7i-k after the onset of buckling.The v 0 field shows compression in x direction both before and after buckling.Similar results as in the virtual experiments are observed in v 0 and v 1 .Note that the large values in v 1 (x) are due to the scaling of the corresponding mode φ 1 (maximum value of 0.04).The auxetic effect is more pronounced in the real experiments (Figs.7j,o) compared to the virtual experiments (Figs. 5j,o), due to the larger specimen size, i.e., the larger scale ratio.In a specimen with a larger scale ratio the buckling pattern in the centre is less affected by the two stiffening lateral edges and, thus, it is more pronounced and results in larger auxetic effect.
The long-range correlated mode, φ 1 (x), as well as the long-range correlated fluctuation fields, v 1 φ 1 (x), in x and y directions, are depicted in their deformed configurations in Figs.7d-f and Figs.7l-n, before and after buckling, respectively.In Figs.7d and 7l, φ 1 (x) is scaled to assure the visibility of the deformed shape of the modes.As expected, small values for v 1 prior to buckling entail a low sensitivity to the fluctuation mode, i.e., a slightly different mode is identified at this stage.Indeed, the long-range correlated fluctuation mode is as good as absent before the emergence of the buckling pattern, where it is not yet representative.The edge effect prior to buckling, which is most pronounced at the corner unit cells, is partially captured by the identified mode.This explains the non-zero values of v 1 and thus v 1 (x)φ 1 (x) in these areas.
The residual fields in normalized grey scales are shown in Figs.7g and 7o, in the deformed configuration based on the evaluated displacement u(x) from Eq. ( 4), before and after buckling.In order to extract the microfluctuation field w(x), the displacement fields obtained by local DIC on the same images are used as the reference.Note that the local DIC procedure is not a part of the methodology introduced here, and is used to assess the microfluctuation field (with errors from local DIC) in order to evaluate the performance of the method.Although local DIC gives the displacement field with minimum kinematical constraints, yet with significant statistical and systematic errors, it requires tens of thousands of dofs to fit the displacement field compared to the 77 dofs for the micromorphic IDIC framework.Figs.7h and 7p show the Euclidean norm of the microfluctuation field, ||w||, before and after buckling.The microfluctuation field w, is of a higher spatial frequency than the characteristic frequency of φ 1 , which is more pronounced far from edges of the specimen in the post-buckling regime (Fig. 7p).Yet, it can be observed prior to buckling as well (Fig. 7h).Beyond buckling, the microfluctuation field w has a significantly smaller amplitude than the correlated fluctuation field.This is not the case prior to buckling, considering the almost zero contribution of the long-range correlated fluctuation field.The residual field is larger in certain areas, which corresponds directly to higher values in the microfluctuation field, both before and after buckling.
All the identified fields, i.e., v 0 (x), v 1 (x) and φ 1 (x), satisfy the expected mechanical behaviour, and the correlation residual, reflected also in the evaluated microfluctuation field w(x), confirms the good performance of the micromorphic IDIC method.

Graded Microstructures and Size Independence of the Modes
As long as the constitutive law of the base elastomeric material is size independent (e.g., the hyperelastic Ogden material used in Section 3.2, usually used for constitutive modelling of metamaterials), the modes φ i are independent of the unit cell size l.To show that the modes are independent also of d/l ratio (i.e., hole diameter to cell size ratio), an infinitely large specimen with a square stacking of holes is considered, cf.Figs.2a-b.In such a case the buckling analysis carried out on a 2 × 2 periodic cell provides the reference (exact) mode-exact within the accuracy of the adopted numerical model-, against which the approximation of φ 1 given by Eq. ( 6) can be compared.Results are obtained for a larger and smaller d/l ratio of, respectively, 0.95 and 0.25, compared to the reference value of 0.79 identified by the micromorphic IDIC procedure, see Fig. 8.The deformed configuration of the reference mode is shown in black colour, whereas the deformed shape of the experimentally identified mode is shown in light grey.Note that only one quarter of the 2 × 2 periodic cell is shown due to various symmetries (see e.g.Bertoldi et al., 2008).Fig. 8a suggests that the overall shape of the mode is captured sufficiently accurately for d/l = 0.95 (recall also Fig. 5), although local discrepancies can be observed in higher-frequency components.Such components are not visible in Figs.4b-c and 6b-c, since their contributions to spectral densities are very small.More interestingly, comparable errors are observed also in Fig. 8b for d/l = 0.79, meaning that identified mode coefficients a φ 1 do not change substantially as a function of small variations of d/l.For a very different ratio, d/l = 0.25, however, larger discrepancies between the two modes are observed in Fig. 8c, suggesting that the experimental mode should be identified again.Although this insight, i.e., that the modes are independent for moderate variations in the hole diameter to cell size ratio, has been demonstrated here for a square stacking of holes, it is a general feature of cellular metamaterials.This brings us to the realization that the introduced micromorphic IDIC technique, combined with the previously developed micromorphic computational homogenization scheme, provides a general tool for the design of engineering systems at two scales.(i) Size effects and microstructural changes are communicated through the patterning modes to the macroscale, where such information can be utilized for the design of engineering systems.(ii) Microstructures with graded properties smoothly varying in space (e.g., d(x)) or material property (e.g., Young's modulus, E(x), or Poission's ratio, ν(x)), can be optimally adjusted and understood at the micro-or meso-scale level.For more details on the modelling and application of functionally graded materials the reader is referred to the extensive literature on this topic, see, e.g., (Anthoine, 2010;El-Helou and Harne, 2019;Francesconi et al., 2019;Farzaneh et al., 2022;Morris et al., 2023).It is important to realize that considerable design freedom is available even for geometrically simple microstructures, such as the herein studied square stacking of holes.Adopting different types of microstructural morphologies (e.g., hexagonal stacking of holes, chiral or re-entrant microstructures) expands the design space further.

Summary and Conclusions
Cellular elastomeric metamaterials generally exhibit a pattern transformations based on microstructural buckling at a critical level of compressive load.These pattern transformations are correlated over long ranges in the specimen and result in considerable changes in the overall mechanical response after their emergence.A micromorphic Integrated Digital Image Correlation (micromorphic IDIC) method, was here developed to experimentally identify and quantify the decomposed kinematics for cellular metamaterials.The decomposition follows a recently introduced micromorphic computational homogenization scheme based on a kinematical ansatz consisting of three parts: (i) mean smooth displacements, corresponding to the slow-scale material level; (ii) long-range correlated fluctuation fields, representing the buckling patterns; and (iii) a microfluctuation field, which captures the remaining uncorrelated part of the fluctuations.This kinematical ansatz was integrated within the digital image correlation scheme, enabling a direct identification of the smooth displacement, and the long-range correlated fluctuation pattern along with its spatial distribution field.On the basis of spectral analysis, a 2D Fourier series approximation of the long-range fluctuation mode with a proper parametrization was performed.The identified patterning modes serve as input parameters for the micromorphic computational homogenization scheme, which is capable of providing mechanical predictions of engineering systems including size effects as well as temporal and spatial mixing of patterning modes.
The methodology is validated on both virtual and real experiments for a specific case of elastomeric metamaterials with rectangular stacking of millimetre-sized circular holes.A standard spectral density analysis is used to attain a regularization with a reduced number of degrees of freedom and to initialize the associated parameters of the long-range fluctuation mode for the considered metamaterial.The method successfully decomposes the displacement field of cellular elastomers under compressive load, both before and after local microstructural buckling leading to long-range correlated pattern transformations.The performance of the method is assessed by identifying the microfluctuation fields through a comparison with local DIC data.In all cases, the microfluctuation field is of higher spatial frequency than the long-range fluctuations and of smaller amplitude after the emergence of the buckling pattern, thereby validating the kinematical ansatz as well as the micromorphic IDIC identification routine.An initial guess robustness study, performed by random perturbations up to 10% in fluctuation mode parameters, resulted in 98% robustness of the correlations and negligible errors (0.45% on average) in the results, confirming the robustness and the accuracy of the proposed method.
The micromorphic IDIC methodology may be considered novel and unique in the following aspects: (i) The methodology aims at and succeeds in identifying long-range correlated fluctuation fields, which is not existing yet in commercial DIC packages.(ii) The decomposition of the displacement field is performed in a single minimization step, identifying the long-range correlated fluctuation mode, its amplitude in space, and the mean smooth displacement field; (iii) The methodology can be readily extended to any cellular metamaterial by making proper initialization choices of the fluctuation modes; (iv) It is easy to implement based on minimal modifications to conventional global DIC formulations.The proposed micromorphic IDIC technique, in combination with the previously developed micromorphic computational homogenization scheme, thus provides an integrated experimental-computational tool for optimizing the design of metamaterial-based engineering systems, including materials with various metamaterial microstructures and/or with graded properties in space, such as the hole size, or material properties, such as the Young's modulus or Poisson's ratio.

Figure 1 :
Figure 1: Elastomeric metamaterial specimen (a) before and (b) after the onset of microstructural buckling.The emerging antisymmetric pattern is clearly visible, which is a result of a compression in the horizontal direction.
Figure 2: Reference configuration, one particular deformed configuration, and spectral densities S xx (ξ) and S yy (ξ), normalized with respect to the unit cell size l, expressed as a function of the wave vector ξ = [ξ 1 , ξ 2 ], corresponding to a square stacking of holes (a)-(d), hexagonal stacking of holes (e)-(h), and a chiral microstructure (i)-(l).The most of the pattern energies localize around a few particular wave vectors, allowing in all cases for effective approximation through a truncated Fourier series.
Figure 4: (a) Virtually generated and deformed speckle pattern, corresponding to compression in horizontal direction after the onset of pattern transformation in a cellular elastomeric metamaterial, revealing the presence of one fluctuation mode inside the specimen.Energy spectral density functions of the corresponding displacement field in (b) x and (c) y direction, revealing the principal wave vectors related to the long-range correlated fluctuation mode (four peaks approximately on the diagonals).The range of the colour bars is reduced to highlight the frequencies related to the long-range correlated fluctuation mode.

Figure 5 :
Figure 5: Results of the virtual experiments before (left) and after (right) the onset of microstructural buckling: (a, b, i, j) the smooth mean displacement fields, v 0 (x), in (a, i) x and (b, j) y directions; (c, k) the pattern amplitudes, v 1 (x), (d, l) the deformed configuration of identified patterns/modes, φ 1 (x), scaled to assure the visibility of the deformed shape, and (e, f, m, n) the resulting fields, (e, m) in x and (f, n) y directions, corresponding to the long-range correlated fluctuations.In (d, l) φ 1 (x) are scaled to assure the visibility of the deformed shape of the modes.The long-range correlated fluctuation fields, v 1 (x)φ 1 (x), are plotted in the deformed configuration corresponding to their contribution only.The residual fields in normalized gray scales are in (g, o), and (h, p) show the Euclidean norm of the microfluctuation field, w(x), plotted in deformed configuration based on the total displacement.
Figure 6: (a) Image of the deformed configuration of a cellular elastomeric metamaterial with millimetresized circular holes in a rectangular stacking, corresponding to a compressive load step after the onset of pattern transformation, acquired during an in-situ compression test.Energy spectral density functions of the displacement field, assessed by local DIC, in (b) x and (c) y direction, revealing the principal frequencies related to the long-range correlated fluctuation mode (four peaks approximately on the diagonals).The range of the colour bars is reduced to highlight the frequencies related to the long-range correlated fluctuation mode.

Figure 8 :
Figure 8: Patterning modes for an infinite microstructure with a square stacking of holes for a unit cell size l = 1.9 mm and a varying hole diameter d ∈ {1.81, 1.50, 0.48} mm.Numerically computed modes are shown in black, whereas experimentally obtained modes according to Eq. (6) in Section 3.3 are shown in light grey.Large hole diameter to cell size ratio d/l = 0.95 is shown in (a), d/l = 0.79 corresponding to the reference geometry analysed in detail in Section 3 is shown in (b), while small hole diameter d/l = 0.25 is shown in (c).