Micromorphic Computational Homogenization for Mechanical Metamaterials with Patterning Fluctuation Fields

This paper presents a homogenization framework for elastomeric metamaterials exhibiting long-range correlated fluctuation fields. Based on full-scale numerical simulations on a class of such materials, an ansatz is proposed that allows to decompose the kinematics into three parts, i.e. a smooth mean displacement field, a long-range correlated fluctuating field, and a local microfluctuation part. With this decomposition, a homogenized solution is defined by ensemble averaging the solutions obtained from a family of translated microstructural realizations. Minimizing the resulting homogenized energy, a micromorphic continuum emerges in terms of the average displacement and the amplitude of the patterning long-range microstructural fluctuation fields. Since full integration of the ensemble averaged global energy (and hence also the corresponding Euler--Lagrange equations) is computationally prohibitive, a more efficient approximative computational framework is developed. The framework relies on local energy density approximations in the neighbourhood of the considered Gauss integration points, while taking into account the smoothness properties of the effective fields and periodicity of the microfluctuation pattern. Finally, the implementation of the proposed methodology is briefly outlined and its performance is demonstrated by comparing its predictions against full scale simulations of a representative example.


Introduction
Mechanical metamaterials exhibiting exotic mechanical properties, such as auxetic behaviour, serve dedicated applications, most notably in soft robotics (Whitesides, 2018;Mirzaali et al., 2018;Zhang et al., 2018;Mark et al., 2016;Jiang and Wang, 2016;Yang et al., 2015). The development of advanced manufacturing techniques, like 3D printing, has enabled a significant increase of the design space for making materials with certain specific architected microstructures (Ren et al., 2018;Truby and Lewis, 2016;Wang et al., 2015). In this contribution, we in particular concentrate on elastomeric materials with periodic microstructures consisting of a regular grid of circular holes, which undergo a pattern transformation triggered by externally applied compression upon reaching a critical compressive strain (Bertoldi et al., 2010). The pattern-transformed material has been shown to reveal a significantly different effective behaviour compared to its initial untransformed state .
The overall effective response of this class of materials depends on the number of holes (which may be intractable) as well as their relative (unknown) position with respect to the specimen's geometry. This is a challenging problem because, depending on the loading and boundary conditions applied to the specimen, the patterning may occur or it may be restricted. Adjacent cells kinematically communicate in this respect and the effective response is intrinsically non-local. A typical example is shown in Fig. 1a, where an elastomeric sheet of size L × L with circular holes subjected to combined compression and shear is depicted. The resulting deformation exhibits a transformed pattern, a shear band, and stiff boundary layers. The boundary layers are present due to the kinematic constraints close to the physical boundaries (the two horizontal edges), where the regular buckling pattern cannot fully develop. The non-transformed material in such constrained boundaries has a higher stiffness, leading directly to stiff boundary layers. These non-uniformities in the deformation pattern give rise to a size effect in the overall mechanical behaviour of the considered system, i.e. the overall response depends on the scale ratio. The scale ratio is, in this particular case, defined as the total height or width of the entire specimen L with respect to the size of one primitive cell (containing only one hole), i.e. as L/ . In order to obtain the effective macroscopic response, use can be made of the periodicity of the microstructure. For the example shown, a microstructural periodic cell consisting of 2 × 2 holes can be evaluated to obtain the homogenized response, as shown experimentally by , and computationally by Ameen et al. (2018) using ensemble averaging. A fruitful starting point for computing the effective response of a heterogeneous material is provided by computational homogenization. Conventional first order computational homogenization (Kouznetsova et al., 2001;Miehe and Koch, 2002;Matsui et al., 2004), can provide accurate results for elastomeric metamaterials only for very large scale ratios (Ameen et al., 2018). Here, a microscale boundary value problem is solved at every macroscopic point to retrieve the macroscopic constitutive response, which is used in a concurrent macroscopic simulation. For smaller scale ratios, however, the long-range fluctuation field, resulting from the kinematic interaction between neighbouring microstructural periodic cells, needs to be captured by the homogenization framework. Various higher-order homogenization frameworks have been proposed in the literature, which help to capture the non-local behaviour. Second-order computational homogenization (Kouznetsova et al., 2004b), for instance, incorporates the gradient of the macroscopic deformation gradient tensor into the kinematic macro-micro scale transition. It entails, however, additional complexity through the solution of a higher-order equilibrium equation at the macroscale, for which the length scale equals the size of the Representative Volume Element (RVE), see e.g Kouznetsova et al. (2004a) for further details.
Another class of extended continua has been proposed via generalized micromorphic theories by Eringen (1968) and Forest (2009). Here, additional kinematic fields are introduced in order to incorporate underlying microfluctuation fields which are not captured by standard kinematic variables. The Cosserat continuum (Cosserat and Cosserat, 1909) is a well-known example, where the additional kinematic field corresponds to local rotations. Hütter (2017) proposed a methodology for homogenization towards a fully micromorphic continuum and derived the relevant generalized stresses. Biswas and Poh (2017) proposed a computational homogenization scheme for matrix-inclusion composites, in which an additional kinematic field characterizing the average strain in the inclusions is introduced, also resulting in a micromorphic continuum. It was shown that this framework is capable of capturing the homogenized response for such materials without a clear separation of length scales. Other relevant contributions can be found in Jänicke and Steeb (2012) or Forest (2002), and in the references therein.
The present paper proposes a novel computational homogenization framework for materials exhibiting long-range correlated fluctuation fields. Unlike the above-mentioned references, the adopted approach directly incorporates the dominant underlying microstructural kinematical fluctuations through a characteristic fluctuation mode. This mode uniquely relates to the morphology of the microstructure, and is not necessarily limited to rotation or higher-order deformation only. Any relevant kinematical interaction between adjacent periodic cells can be captured. To this end, first the displacement field is decomposed into three sub-fields: a mean effective field, a long-range spatially correlated field, and a local microfluctuation field, cf. Figs. 1b and 1c. In the next step, an entire family of translated microstructures relative to the specimen geometry is considered, the average of which defines the effective behaviour. The entire problem is then formulated variationally, which upon minimization yields the associated Euler-Lagrange equations of the effective micromorphic continuum in terms of the mean displacement and the amplitude of the correlated part of the fluctuation, along with the associated boundary conditions. To avoid having to solve for the microfluctuation field in the entire domain, which effectively makes the solution computa-tionally as expensive as full-scale simulations, two approximations are considered. The first one neglects the uncorrelated microfluctuation field completely. This works well in terms of kinematic quantities, but it turns out that the conjugate quantities are very inaccurate. A second approach is therefore elaborated, in which an approximate microfluctuation field is computed only locally, inside a microstructural periodic cell. To render the computation more efficient, approximations of the local average energy density in the neighbourhood of the considered macroscopic Gauss integration points are further considered (with respect to the microstructural translations), while taking into account the smoothness properties of the effective fields and periodicity of the microfluctuation field. The performance of the proposed computational homogenization framework is evaluated using a representative example, by comparing against the Direct Numerical Simulations (DNS) published by Ameen et al. (2018). The effective kinematic quantities are thereby captured with adequate accuracy. Moreover, the homogenized nominal stresses now also closely predict the size effect exhibited by the DNS solution, achieving a maximum relative error of less than 10 % for the smallest scale ratio considered.
The paper is divided into five sections as follows. The next section, Section 2, describes the full-scale problem considered, its corresponding DNS solution, and the methodology used to define the mean (reference) solution using ensemble averaging. Section 3 then presents the decomposition of the kinematics into the three sub-fields and provides the derivation of the full generalized homogenization framework, which leads to the emergence of an effective micromorphic continuum. Because the resulting problem is still computationally expensive, Section 4 introduces a simplified approach in which the microfluctuation field is neglected and considered as a truncation error. The results obtained with this approach suggest the need to incorporate the fast fluctuating field, since it contributes significantly in terms of energy. Section 5 therefore elaborates on a computationally integrated formulation, in which the microfluctuation field is computed only locally, inside periodic cells associated with individual Gauss integration points. This formulation is dual to the computational homogenization solution of the micromorphic continuum. The performance of the proposed methodology is compared against full-scale simulations in Section 6, and the paper closes with a summary and conclusions in Section 7.
Throughout the paper, the following notational conventions will be used scalars a, vectors a, second-order tensors A, fourth-order tensors 4 C, matrices A and column matrices a, derivatives of scalar functions with respect to second-order tensors

Problem Statement
A representative example consisting of an elastomeric metamaterial exhibiting a patterned deformation under compressive loading is described in Section 2.1, which is used as the reference problem. Although the developments made throughout this paper are only demonstrated on this simple example, the entire methodology proposed is equally applicable to more complex cases. The full-scale problem for a low scale separation regime, along with the DNS solutions are briefly discussed in Section 2.2. Ensemble averaging and the adopted definition for the effective reference solution are presented in Section 2.3.

Reference Problem Definition
The reference problem geometry is sketched in Fig. 2a. A semi-infinite specimen Ω S = (−∞, ∞) × [−L/2, L/2] is considered. Since the specimen geometry is infinitely wide (along the horizontal axis e 1 ) and the applied deformation is uniform in this direction, a model domain Ω = [− , ] × [−L/2, L/2], with Periodic Boundary Conditions (PBC) along the two vertical edges (i.e. along AD and BC), is considered. This is justified experimentally and computationally, since the deformation field typically exhibits a 2 × 2 repetitive pattern, see e.g.  or Ameen et al. (2018). The specimen is subjected to uniaxial compression by prescribing vertical displacements along the two horizontal edges Γ D = AB ∪ CD ⊂ ∂Ω, as u D (on the top edge AB) and zero (on the bottom edge CD), while constraining the horizontal displacements on both these edges. The considered domain consists of a number of square single-hole primitive cells, the size of which is the characteristic microstructural length . The macroscopic length, on the other hand, is specified by the height of the specimen L, giving rise to a scale ratio L/ . In what follows, only positive integer scale ratios larger than 3 will be considered, i.e. L/ ∈ N >3 .
The constitutive behaviour of the elastomeric base material is assumed to be hyperelastic, described by the strain energy density whereas the individual holes are captured through the topology of the microstructure. In Eq.
(1), X ∈ Ω, X = X 1 e 1 + X 2 e 2 is the position vector in the reference configuration, F = I + ( ∇ u) T is the deformation gradient tensor, ∇ indicates the gradient with respect to the reference configuration, u( X) is the displacement field, J = det F , I is the second-order identity tensor, and I 1 = tr C and I 2 = 1 2 (tr C) 2 − tr (C 2 ) are the invariants of the right Cauchy-Green deformation tensor C = F T · F .
The constitutive and geometric parameters adopted throughout this contribution are listed in Tab. 1, based on the experimental characterization by . Note that the same parameters have been used for generating the full-scale DNS solutions reported in Ameen et al. (2018).

Full-scale Problem
For the above hyperelastic problem, the solution (displacement field) u( X) is obtained by minimizing the total potential energy of the entire system, i.e. u( X) ∈ arg min u( X)∈U (Ω) E( u( X)), where χ is the indicator function of the base material, i.e. χ = 0 inside holes and χ = 1 outside (implying also that all internal holes remain stress-free), and U (Ω) denotes the space of kinematically admissible displacement fields over the domain Ω (i.e. the vector space of all vector functions that satisfy the boundary conditions applied on the two horizontal edges, while being periodic in the horizontal direction with the period 2 ). In anticipation of ensemble averaging in Section 2.3, where knowledge of the displacement field u( X) at each position X ∈ Ω (including holes) is required, the entire rectangle Ω = [− , ] × [−L/2, L/2] is considered as the domain of u, i.e. u( X), X ∈ Ω. To this end, displacements and strains inside holes are interpolated using an ultra-soft linear-elastic material, whereas all stress components are set to zero there. As in Eq. (2), throughout this manuscript admissible test functions are denoted with hats ( •), whereas minimizers are free of hats. The inclusion sign ∈ in Eq. (2) emphasizes the non-uniqueness resulting from non-convexity of the total potential energy E, which allows for multiple solutions, instabilities, and buckling.
The typical behaviour of the solutions for the considered problem is shown in Fig. 2. Here, an approximately bi-linear nominal stress-strain response is observed in Fig. 2c, the individual linear regimes being separated by a critical strain of about 3 %. At this strain, the microstructure buckles and nucleates a regular deformation pattern, as shown in Fig. 2b, which results in a significant change in nominal stiffness. Due to kinematic constraints close to the physical boundaries, the regular buckling pattern cannot fully develop. This implies that such constrained boundaries directly lead to more stiff boundary layers, as the non-transformed material has a higher stiffness. Clearly, considering more complex loading scenarios would entail more complicated mechanisms, including boundary layers, localization bands (recall Fig. 1) and percolation paths. The presence of such non-uniformities in the deformation pattern furthermore gives rise to a size effect in the overall mechanical behaviour of the system, i.e. a dependence of this overall response on the scale ratio L/ -particularly for ratios L/ < 10. Because most of these phenomena result from the non-linearities and kinematic interactions between individual neighbouring periodic cells, the considered problem is challenging from a homogenization point of view.

Ensemble Averaging
One might not be, in general, able to control the positioning of the boundary relative to the microstructure-and the other way around. Hence, all possible realizations are considered and the reference mean solution is defined as their ensemble average. In order to establish an effective (homogenized) solution, from which all the fluctuations due to the microstructure are eliminated, a family of minimization problems with translated microstructures similar to Eq. (2) are therefore considered. The translation of the microstructure relative to the specimen's boundaries is denoted by ζ ∈ Q, ζ = ζ 1 e 1 + ζ 2 e 2 , where Q = [− , ] × [− , ] is the corresponding 2 × 2 periodic cell. The periodic cell should not be confused with the primitive cell (defined as [0, ] × [0, ]), which reflects the geometrical periodicity in the undeformed state. The periodic cell has the period of the fluctuating pattern, which is important for the homogenization, and which relies on experimental and numerical evidence. In general cases, Bloch analysis (cf. e.g. Geymonat et al., 1993;Singamaneni et al., 2008) can be used to estimate the wavelength of the pattern. The minimization problems associated with the individual translated microstructures are specified as By considering a fixed point X in the energy density Ψ( X, ζ, F ) of Eq.
(3) we see that this point is occupied by different material points for different translations ζ. Consequently, for certain translation of the microstructure, X will be positioned inside a hole, whereas for other translations in the base material of the elastomeric matrix. Following the argumentation of Smyshlyaev and Cherednichenko (2000) and Cherednichenko and Smyshlyaev (2004), the effective solution is obtained as which corresponds to pointwise ensemble averaging (i.e. for each fixed X ∈ Ω) over all the translated microstructures with a uniform probability density over Q.
The effective solution u can be obtained by considering brute force calculations via DNS and discretisation of the integral in Eq. (4), as reported in detail by Ameen et al. (2018) for various scale ratios L/ . Fig. 3a shows only the components u 2 (X 2 ) of the displacement fields u( X, ζ) for three vertical translations ζ ∈ { 0, 2 e 2 , 3 4 e 2 }, and the corresponding component of the effective field u. Due to e 1 -periodicity, the effective solution u is constant along the horizontal direction in this particular case. Recall that inside the holes, an ultra-soft linear-elastic material has been considered for interpolation purposes.
The local magnitude of the fluctuations for the ensemble u( X, ζ) can be estimated pointwise as a quantity which will be referred to as standard deviation in what follows, shown as the red smooth curve in Fig. 3b. This field is again constant along the e 1 direction. The fluctuating part, approximately determined as [ u( X, ζ) − u( X)]/σ( X) for each translation ζ, is plotted in Fig. 3c for ζ = 0, and essentially corresponds to a microstructural buckling mode. The irregularities in the vicinity of the two horizontal edges originate from the fact that the individual fields u( X, ζ)− u( X) as well as corresponding standard deviation σ( X) are close to zero there, resulting in numerical inaccuracies when dividing them. Since individual periodic as well as primitive cells of the microstructure mutually interact through the pattern, the amplitude of the mode (estimated as σ( X)) and its spatial gradient are important for the effective behaviour of the system. Accordingly, this mode magnitude will be considered as part of the effective solution in addition to u( X) in the homogenization method developed in what follows.

Homogenization Towards a Micromorphic Continuum
Based on the insights obtained from the DNS results presented, a homogenization approach is developed in this section, which consists of three steps. First, a suitable ansatz that accurately decomposes the kinematic fields corresponding to all translated microstructures is developed in Section 3.1. It serves as a basis to build an approximation space over which the effective solutions are sought. In the next step, in Section 3.2, the ensemble averaged energy is introduced, the minimum of which corresponds to the average over all energies with translated microstructures considered at their solutions u( X, ζ), i.e. at their minima. This average energy functional can be minimized directly with respect to the kinematic ansatz introduced in the first step, providing the closest possible approximation to the effective solution u( X) within the kinematic subspace considered. Using the standard arguments of the calculus of variations, a set of governing Euler-Lagrange equations along with the associated boundary conditions are finally derived and discussed in the third step in Section 3.3.

Kinematic Ansatz
The displacement field, common to all translated microstructures according to the definition in Eq. (3), can be generally decomposed as Here, for fixed ζ, u( X, ζ) is the displacement field corresponding to a particular microstructure translated by ζ (i.e. the minimizer of the energy E in Eq. (3)). For a fixed position in the specimen X, we may similarly interpret u( X, ζ) as the displacement at that particular point as a function of the translation vector ζ. The vector fields ϕ i ( X, ζ), i = 1, . . . , n, are long-range spatially correlated modes characteristic of the problem morphology and deformed underlying microstructural pattern (translated by ζ), which are in addition scaled by their magnitudes v i . Individual modes ϕ i ( X, ζ) are periodic and translate along with the microstructure, i.e. ϕ i ( X, ζ) = ϕ i ( X + ζ, 0), and have zero mean. As a result of this, and of the zero mean of w( X, ζ) (see below), v 0 ( X) corresponds to the mean effective displacement field. This decomposition enables the fields v i to 'switch off' individual modes near the boundaries, for example, and thus serve as 'amplitude modulators' for the individual modes ϕ i ( X, ζ). For a motivation of such a general decomposition see for instance Okumura et al. (2004), where multiple bifurcations in foams have been studied.
For the reference problem considered in this manuscript, only one spatially correlated mode is adopted, i.e. n = 1, which is defined as see Fig. 3d. In Eq. (7), C 1 = 1 |Q| Q ϕ 1 ( X, 0) 2 d X is a normalizing constant ensuring that the standard deviation of ϕ 1 ( X, ζ) equals one, recall Eq. (5). One then expects v 1 to approximately recover the shape of the standard deviation shown in Fig. 3b, while v 0 is assumed to capture the trend of the mean, as depicted by the thick red line in Fig. 3a. The component v 1 ϕ 1 essentially establishes the kinematical interactions between individual neighbouring periodic cells through the mode ϕ 1 . For the present case these interactions can be interpreted as a micro rotation field, cf. Figs. 3c and 3d. The mode amplitude is therefore the highest in the bulk of the specimen, and due to the imposed constraints, equals zero at the top and bottom boundaries. The mode is moreover active only after the specimen undergoes instability, whereas in the pre-bifurcation regime it is expected to be inactive (i.e. v 1 vanishes).
Finally, w( X, ζ) is the fast and local microfluctuation field complementary to the sum of the smooth and spatially correlated parts v 0 ( X) + v 1 ( X) ϕ 1 ( X, ζ). The w component is expected to capture fluctuations in the linear regime (in which v 1 vanishes), and to compensate any inaccuracies of the mode ϕ 1 in the non-linear post-bifurcation regime. Because the split in Eq. (6) is not unique, additional conditions are required. Considering a fixed point in space X, both effective fields v 0 ( X) and v 1 ( X) are constant as a function of ζ in Eq. (6), meaning that w should be orthogonal to a constant function in ζ (represented by v 0 ), and to ϕ 1 (which is multiplied by a constant arbitrary magnitude v 1 ). Expressed mathematically, this means that where 1 = e 1 + e 2 is a constant vector in ζ. An extension to multiple modes ϕ i , i = 1, . . . , n > 1, requires orthogonality with respect to the entire vector space spanned by the considered modes.

Energy Considerations
The energy associated with one fixed realization ζ ∈ Q of the translated microstructure reads whereas the ensemble averaged energy may be written as (see also Smyshlyaev and Cherednichenko, 2000;Cherednichenko and Smyshlyaev, 2004;Peerlings and Fleck, 2004) 1 Note that in Eq. (10), the minimization with respect to u( X, ζ) in the first equality is considered for all X and only one particular fixed ζ, whereas in the second equality ζ is considered as a continuous variable equally to X; this means effectively that u( X, ζ) is a vector function over a four-dimensional space (in the case of Ω ⊂ R 2 and Q ⊂ R 2 as considered here). In the second equality of Eq. (10), the average energy functional has also been introduced as The decomposition introduced in Eq. (6) is next used as a basis for the minimization of the variational homogenization problem by considering the set of admissible functions in the form u( X, ζ) = v 0 ( X) + v 1 ( X) ϕ 1 ( X, ζ) + w( X, ζ), where the present problem is limited to one mode (n = 1) only. The problem in Eq. (10) is then effectively replaced with min u( X, ζ)∈U (Ω) E( u( X, ζ)), from which the kinematical functions v 0 ( X), v 1 ( X), and microfluctuation functions w( X, ζ) can readily be obtained as minimizers, i.e.
The minimizer v 0 results in as close an approximation to u as possible within the test space U (Ω) considered, i.e. within the decomposition according to Eq. (12). So far the exact solution for all ζ can be represented because w still allows to cover the entire space U (Ω). Later, however, this space will be restricted by adopting certain constraints on w, which introduces an approximation.

Euler-Lagrange Equations
Following the standard argumentation of variational calculus, we take in the next step the Gâteaux derivative of E and require it to vanish, which provides where the fully-expanded definition of the displacement gradient for the considered kinematic field of Eq. (12) reads and, accordingly, Before proceeding, let us note that one could, at this stage, make use of the divergence theorem in the last row of Eq. (15) to obtain 1 |Q| Ω Q (− ∇ · P T ) · δ u d ζd X + boundary terms = 0, which requires that the balance law ∇ · P T = 0 holds for each realization ζ. This condition is satisfied only if w is unconstrained (i.e. the decomposition of Eq. (12) covers the entire kinematic space U (Ω)), meaning that Eqs. (15) and (18) reduce to a classical continuum. This would require computing the exact solution everywhere in Ω for each translation ζ. Below we will restrict the space in which w is considered and require only that the ensemble average vanishes, meaning that only the ensemble average of the last row in Eq. (15) is zero and that individual realizations may not be in equilibrium under constraints on w imposed. In order to establish ensemble averaged balance laws for the homogenized quantities, i.e. stress quantities averaged over Q and tested by variations of the effective fields δ v 0 and δv 1 , the last row of Eq. (15) is rewritten as In compact form, Eq. (19) reads Ω P 1 ( X) : ∇δ v 0 ( X) + P 3 ( X) · ∇δv 1 ( X) + P 2 ( X)δv 1 ( X) which, with the help of the divergence theorem applied to ∇δ v 0 , ∇δv 1 , and ∇δ w, can be rewritten as where Γ N denotes the free part of the domain boundary ∂Ω, and N the unit outer normal to ∂Ω in the reference configuration. Recall, however, that prescribed tractions have been omitted in the definition of the ensemble averaged energy of Eq. (11), and that only essential boundary conditions on the Γ D part of the boundary along with periodic boundary conditions have been considered. Making use of the localization argument in space as well as in translations (because w is a space-translation quantity), Eq. (21) gives the following set of governing Euler-Lagrange equations δ v 0 : ∇ · P T 1 = 0, in Ω, P 1 · N = 0, on Γ N , δv 1 : ∇ · P 3 − P 2 = 0, in Ω, δ w : ∇ · P T = 0, in Ω, where the averaged stress quantities P 1 ( X), P 2 ( X), and P 3 ( X), have been introduced as Although the stress quantity P 1 is a simple ensemble average and its balance law in Eq. (22a) has the standard form, the balance law is not standard as P 1 depends also on v 1 and ∇ v 1 . Through these two terms the balance law in Eq. (22a) couples to the non-standard equilibrium Eq. (22b), expressed in terms of the two non-standard stress quantities P 2 and P 3 . Since all homogenized stresses depend on (gradients) of v 0 and v 1 , a micromorphic continuum emerges. So far, all stresses depend on ∇ w as well. As Eq. (22) has been derived for the ensemble averaged quantities v 0 and v 1 , it depends on the coordinate vector X only. In contrast to that, w and its variation δ w are also functions of ζ, and hence δ w cannot be taken outside the integral over Q in the derivation of Eq. (21). This effectively means that w needs to be solved in Eq. (23) for each translation ζ separately.
Unlike the natural boundary conditions expressed in terms of P 1 , P 2 , and obtained from the minimization procedure of the ensemble averaged energy (cf. Eqs. (21) and (22)), the essential boundary conditions for all coarse fields v 0 and, in general, v i can be derived directly from the considerations of individual translations, recall Section 2.2 and Fig. 3. This may generally not be as straightforward for other homogenization techniques based on, e.g., volume averaging. In the case of the developed micromorphic homogenization technique it is clear, however, that when the displacement field is prescribed on a part of the domain boundary Γ D ⊂ ∂Ω, i.e. u( X, ζ) = u D ( X) on Γ D , all fluctuations v i ( X) ϕ i ( X, ζ) need to vanish there for all translations ζ ∈ Q. This effectively means that all coarse modulating functions vanish at this part of the boundary, i.e. v i ( X) = 0 on Γ D , whereas v 0 ( X) = u D ( X) on Γ D captures any prescribed displacements.
In total three unknown fields are considered in Eqs. (22) and (23), namely v 0 ( X), v 1 ( X), and w( X, ζ), which will be discretised using the finite element method. Assuming that v 0 ( X) and v 1 ( X) are slowly varying effective fields, their spatial approximation or discretisation may be relatively coarse, considered at the level of the entire specimen L. On the contrary, w( X, ζ) is a rapidly oscillating field, which needs to be resolved finely at the microstructural level for each translation ζ ∈ Q. In order not to solve for w( X, ζ) globally inside the entire domain Ω (which would effectively correspond to the DNS solutions for all translated microstructures), two approaches are presented below: • First, w( X, ζ) is neglected and considered as a truncation error in the closed-form homogenization approach detailed in Section 4. This is only a reasonable assumption if a high approximation quality can be achieved with the decomposition v 0 + v 1 ϕ 1 , i.e. when ϕ 1 captures most of the significant fluctuations present.
• Second, w( X, ζ) is computed only locally inside a periodic cell associated with each coarse-scale Gauss integration point for only one translation, reducing to a computational homogenization approach. This is discussed in more detail in Section 5.

Closed-Form Homogenization Approach
The first approach adopted neglects the microfluctuation field w (i.e. w is constrained to w( X, ζ) ≡ 0), assuming that the leading terms v 0 +v 1 ϕ in Eq. (12) approximate the macroas well as micro-kinematics sufficiently accurately. Upon such a simplification, the governing Eqs. (22) and (23) effectively reduce to Eq. (22) for the homogenized quantities only, along with the definitions of the corresponding conjugate stresses in Eq. (24). Consequently, the resulting problem can readily be discretised and solved as follows.

Discretisation of the Governing Equations
A macroscopic discretisation is introduced, and both effective fields are expressed in terms of shape functions as where v 0 and v 1 are columns storing nodal values of the respective fields, and N 0 ( X), N 1 ( X), are matrices storing the individual shape functions. Spatial derivatives are then expressed as where B 0 ( X) and B 1 ( X) are standard matrices of the shape functions' spatial derivatives. The variations δ v 0 , δv 1 , and their spatial derivatives are approximated in the same way. Making use of the first variation of the average energy functional E in Eq. (20) and the definitions of homogenized stresses in Eq. (24), one can write for the column matrix of internal forces where A denotes the assembly operator over n e elements with spatial subdomains Ω e . For the average energy itself, one can write In both Eqs. (27) and (28), the integrals over individual elements Ω e are carried out numerically in the standard way using a Gauss integration rule, as indicated in Eq. (28). Here, the total average energy is computed as a sum over all integration points, where w ig is the integration weight, J ig is associated Jacobian, and Ψ ig denotes the average energy density associated with the considered integration point i g . Note that when the microfluctuation field w is not present, all average conjugate quantities P 1 , P 2 , P 3 , and the average energy density Ψ, can be computed by simply sampling the corresponding constitutive law (no solution at the 'microscale' is required). Let us consider a fixed Gauss integration point i g , in which the coarse fields ∇ v 0 , v 1 , and ∇v 1 are known. Because also a fixed position is assumed, no further function dependencies on X need to be considered. As a consequence, the deformation gradient becomes an explicit function of ζ only, i.e.
where ϕ 1 ( ζ) is an a priori known vector field, recall Eq. (7). The deformation gradient Eq. (29) can thus directly be substituted into the adopted constitutive law and integrated over Q, i.e.
Similar expressions also hold for the remaining homogenized quantities P 2 and P 3 . In Eq. (30) we have, without any loss of generality, considered the position of the Gauss integration point as X = 0. For the numerical results presented below in Section 4.2, piece-wise linear shape functions have been used with a one-point integration scheme, for both effective fields v 0 and v 1 . Because the example considered is periodic in e 1 direction, the effective problem reduces to be one-dimensional only, in which all averaged quantities are constant along the e 1 direction. Although a standard numerical Newton solver can be used to minimize E, or in other words to equilibrate f , for convenience and for brevity a quasi-Newton solver has been used instead. This minimization technique does not require any evaluations of the Hessian of the average energy E (or in other words the tangent stiffness matrix associated with f ), but more iterations are usually required.

Results
The effective solutions v 0 (X 2 ) and v 1 (X 2 ) obtained for two scale ratios L/ = 5 and L/ = 12 and for applied nominal strains of u D /L = 0.02 and 0.075 are compared against the corresponding characteristics of the respective DNS solutions in Fig. 4. Here we can see that the achieved accuracy in terms of kinematics is encouraging, especially in the post-bifurcation regime (Figs. 4b and 4d). The boundary layers are captured accurately, scaling well with the adopted scale ratio L/ . The plateau of v 1 in the bulk region, especially for the scale ratio 12, is also well reflected. The small differences between v 1 and the standard deviation σ are explained by the fact that the two quantities are not directly comparable (recall that whereas σ includes the effect of w, v 1 does not). For both homogenized solutions in the pre-bifurcation regime, the term v 1 ϕ 1 tries to represent the effect of w, mimicking a buckled state in which v 1 = 0 (Figs. 4a and 4c); this is a result of the overly stiff (over-constrained) representation, as explained in the next paragraph.
Although the kinematic fields are relatively acceptable, especially in the post-bifurcation regime, the corresponding conjugate quantities are severely inaccurate, see Fig. 5 where a logarithmic scale has been adopted on the vertical axis to allow for comparison. Here, the (mean) vertical stress component (P 1 ) 22 which results from the applied deformation as a function of the scale ratio L/ is shown. In the homogenized solution this quantity is simply the outcome of the numerical analysis, whereas in the DNS case it has been obtained by brute force averaging over all realizations. The dependence of the results on the scale ratio L/ originates from restricting the pattern transformation in a boundary layer (thickness of which scales with ) close to the top and bottom boundaries. In the limit L/ → ∞ this influence becomes negligible and a scale-independent limit is reached. The closed-form homogenized solution is able to qualitatively capture this trend in the post-bifurcation regime, although the values differ by a factor of 15 (Fig. 5b). Pre-buckling, however, the size effect exhibited by the closed-form homogenized solution is artificial (and wrong), as almost no dependence on the scale ratio is observed for the DNS results, see Fig. 5a. This is due to the fact that actually almost all homogenized solutions buckle already even for applied strains as small as u D /L = 0.01. As mentioned above, premature bifurcation occurs to relax the energy through the long-range correlated field v 1 ϕ 1 to compensate for the neglected microfluctuation field w. This is the result of overconstraining w, which leads to an overly stiff system compared to the DNS. Note that one option to improve on this is to consider more modes ϕ i with shorter wavelengths, capable of approximating the vector field w in a better way. This approach is not considered in what follows, and a computational scheme is adopted instead.

Computational Homogenization Approach
In spite of the reasonable effective post-bifurcation kinematic fields obtained under the constraint w ≡ 0, Section 4.2 convincingly showed that the microfluctuation field w needs to be accounted for. As briefly mentioned at the end of Section 3, the difficulty with w is twofold: (i) it is defined on the entire (full-scale) problem domain Ω, meaning that one has to basically solve the full-scale problem to obtain it, and (ii) one needs to do this for every realization (translation of the microstructure). Hence, in order to avoid resorting to DNS type of simulations, additional approximations on w need to be made to simplify the multi-scale ansatz. A first approximation which we propose here is to compute (for each translation) the microfluctuation field w on a 2 × 2 periodic cell around the point of interest (e.g. a macroscopic integration point) instead of on the full scale problem domain Ω. This means that boundary conditions have to be applied on the boundary of that periodic cell instead of on the outer boundary of the problem domain. (Note that near the specimen boundary ∂Ω, part of this periodic cell may be outside of the problem domain Ω.) To achieve this, two assumptions are made: • On the periodic cell, the effective fields v 0 and v 1 vary only slowly and hence they can be approximated by a linear expansion around the point of interest (i.e. the centre of the periodic cell); • The remainder w is periodic over the periodic cell.
These conditions, together, allow us to estimate the energy density at a point X for every translation ζ.
A second approximation is made in order to avoid the computation of w for every translation ζ. Hereto, we estimate the energy density in a point X for all the translations from the energy density distribution in the neighbourhood of that point in one realization only.
All considered assumptions are elaborated formally in Section 5.1. Approximations of the average energy along with its minimization leading to the governing equations are further discussed in Section 5.2, whereas the overall computational framework is summarized in Section 5.3. For results obtained with this framework we refer to Section 6.

Assumptions and Approximations
Using a Taylor series expansion and assuming that the v 0 and v 1 fields vary spatially at a much coarser scale than ϕ 1 and w, the first assumption implies that the effective fields can be approximated in a small neighbourhood of a point X (spanned by ∆ X) as The microstructural components of the ansatz Eq. (6), on the other hand, fluctuate at a fine scale compared to v 0 and v 1 and hence, they need to be evaluated exactly as ϕ 1 ( X + ∆ X, ζ), and w( X + ∆ X, ζ).
Substituting Eqs. (31) and (32) into Eq. (12) provides an expression for the approximate displacement field in the form e 1 e 2 X ζ Ω X X + ζ ζ Figure 6: A sketch showing a point X in the reference microstructure ζ = 0 (in black) and a microstructure translated by ζ (in blue). The point X in the translated microstructure is microstructurally (or pattern-wise) at the same position as the point X + ζ in the reference microstructure, although these two points experience slightly different effective fields v 0 and v 1 .
with corresponding deformation gradient The differentiation with respect to a (fine scale) variable ∆ X in Eq. (34), i.e. ∇ m = ∂/∂∆X i , is valid because within the small neighbourhood considered, the global position X is kept fixed.
In the second step, two realizations as sketched in Fig. 6 are considered. The reference microstructure is defined as ζ = 0 (in black), whereas the translated microstructure is defined as ζ = 0 (in blue). A point X in the microstructure translated by ζ has exactly the same relative microstructural position as the point X + ζ in the reference microstructure, as indicated in the right part of Fig. 6. This means that these two points are microstructurally (or pattern-wise) identical, although they are experiencing slightly different underlying effective fields v 0 and v 1 . Using this second assumption, the constitutive law of a point X in a translated realization exactly corresponds to that of a point X + ζ in the reference realization, i.e. Ψ( X, ζ, F ) = Ψ( X + ζ, 0, F ), where an arbitrary deformation gradient F has been considered. Next, a similar procedure is applied to the deformation gradient at a point X for the translated microstructure, i.e. F ( u( X, ζ)), which is approximated by the deformation gradient in the corresponding position X + ζ of the reference microstructure. To this end, use of the Taylor-expanded deformation gradient in Eq. (34) is made, which provides F ( u( X, ζ)) ≈ F ( u( X, 0, ζ)) ≈ F ( u( X, ζ, 0)).
In the first approximation of Eq. (36), the exact deformation gradient at a position X and microstructure translated by ζ has been replaced by the Taylor-expanded one from Eq. (34) (as the point X is considered, ∆ X = 0 while ζ = 0), whereas in the second approximation of Eq. (36) the corresponding point X + ζ in the reference microstructure is used (i.e. the translation vector ζ effectively becomes a position vector relative to X, replacing ∆ X in Eq. (34)). Translating the constitutive law (Eq. (35)) while employing the kinematic approximation of the deformation gradient (Eqs. (34) and (36)), one can finally estimate the energy density at a point X in the microstructure translated by ζ as Ψ( X, ζ, F ( u( X, ζ))) = Ψ( X + ζ, 0, F ( u( X, ζ))) ≈ Ψ( X + ζ, 0, F ( u( X, ζ, 0))).
The rightmost form of the energy density in Eq. (37) clearly needs to be evaluated for the reference microstructure only (i.e. ζ = 0), over a small spatial neighbourhood of a point X (in fact, over the periodic cell Q).
The assumption of periodicity of the microfluctuation component requires that w is periodic, superposed on the approximated v 0 + v 1 ϕ 1 field (recall Eq. (33)), meaning that holds, whereas the four corner points P i of the periodic cell Q are fixed, i.e. w(P i ) = 0, i = 1, . . . , 4. In Eq. (38), Γ i , i = 1, . . . , 4, denote boundary segments of the periodic cell Q, cf. Fig. 8c. In order to alleviate uniqueness concerns (recall Eq. (8) and the discussion therein), the orthogonality condition ∇ w, ϕ 1 2 = 0 still needs to be enforced in addition to the periodicity of w, although it is now considered in space rather than in translations.

Euler-Lagrange Equations
It is important to realize at this point that by the three assumptions made above, the integration over all translations ζ effectively reduces to a spatial integration over some local neighbourhood in space in one particular realization. To indicate this more clearly, a change of variables X m = ζ is employed hereafter. The integration over Q then reduces to a spatial integration over Ω m , which again spans the 2 × 2 periodic cell, and the differentiation ∇ m with respect to the local variable ∆ X in Eq. (34) reduces to ∇ m = ∂/∂X m,i . For further convenience, and in order to establish a link to computational homogenization, Ω m will be referred to as RVE, which is translated such that its centre point is located at a 'macroscopic' position X. Furthermore, X m will be referred to as the 'microscopic' coordinate.
Taking all considerations into account, the approximate average energy takes the form in analogy to the exact average energy functional defined in Eq. (11). Setting the first variation of E equal to zero provides which differs from the variation of the exact average energy functional in Eqs. (15) and (16) in only one term, namely in X m · ∇δv 1 ∇ m ϕ 1 , and in the fact that δ w is periodic. Following therefore exactly the same procedure as in Section 3.3, the set of governing equations for v 0 and v 1 is derived, cf. Eq. (22), where the definitions of the homogenized stresses now read In Eq. (41), the stress quantities P 1 and P 2 are analogous to their exact counterparts in Eq. (24), whereas P 3 has an extra term due to the extra contribution X m · ∇δv 1 ∇ m ϕ 1 in δE, cf. Eq. (40). The spatial dependence of individual quantities in Eq. (41) is to be understood as follows. Each macroscopic point X has assigned its own RVE, Ω m , over which the local stress P m ( X, X m ) is defined; hence the dependence of P m on both variables X and X m . On the other hand, the fluctuation field ϕ 1 is considered relative to the RVE's centre point, which is assumed to be 0 for convenience (recall that all RVEs are translated in space such that their centre points coincide with X); hence ϕ 1 eventually only depends on X m . The same explanation also holds for w. Note that the adopted approach may reflect situations in which individual RVEs (and hence also their fields ϕ 1 and w), corresponding to a pair of closely positioned macroscopic points, overlap.
At the microscale, the following equation holds where the orthogonality constraint w, ϕ 1 2 = 0 and periodicity of w over Ω m need to be enforced in addition (recall Eqs. (8) and (38), and the discussions therein). The orthogonality condition is now considered in space over Ω m instead of in translations ζ, as microstructural translations ζ have been reduced to changes in spatial position X m relative to X in Section 5.1. The periodicity of w and δ w further entails that the boundary term, similar to the one in Eq. (23) and now expressed as P m · N m = 0 on Γ m , disappears. From Eq. (42) it is furthermore clear that all microscopic quantities are considered for one realization only. a Solution of the RVE problem requires an independent (Newton, or other path-following algorithmic) loop with assembly of the microscopic gradients g m (and Hessians H m ), a result of contributions from all microscopic Gauss integration points. In addition, orthogonality (i.e. w, ϕ 1 2 = 0) and periodicity constraints over Ω m are enforced for w.

Computational Homogenization Framework
The outline of the solution scheme adopted for the multiscale computational homogenization is given in Algorithm 1 and the overall procedure is sketched in Fig. 7. At each macroscopic Gauss integration point, three quantities are sampled at the macroscale and prescribed to the microscale, namely ∇ v 0 , v 1 , and ∇v 1 (note that v 0 is irrelevant as it does not affect the energy). At the microscale, i.e. inside Ω m , the mode ϕ 1 is introduced, and the underlying smooth field constructed. Considering the necessary orthogonality conditions, the microfluctuation field w is computed, which allows for the computation of local stresses P m . The homogenized stresses are subsequently evaluated according to Eq. (41), and passed on to the macroscale. Here, the global governing equations are assembled and solved.
Note that exactly the same set of multiscale governing equations can be obtained by using arguments of computational homogenization, see e.g. (Kouznetsova et al., 2004b). Starting directly from a two-scale decomposition of the form u( X m , X) ≈ v 0 ( X) + X m · ∇ v 0 ( X) + v 1 ( X) + X m · ∇v 1 ( X) ϕ 1 ( X m , 0) + w( X m , 0), (43) where ∇ v 0 , v 1 , and ∇v 1 are macroscopic quantities assumed to be constant within each RVE, one can arrive to Eqs. (22), (41), and (42), by following standard procedures of computational homogenization. We remark, however, that the derivation based on microstructural translations adopted in this paper is (until making approximations on w) more rigorous, and that alternative assumptions may lead to implementations which are not comparable with computational homogenization and thus leave room for improvement. Such considerations are, nevertheless, outside the scope of this contribution.

Results and Comparison with Full Scale Simulations
Following the same procedure as the one outlined in Section 4.1, both coarse quantities v 0 and v 1 are discretised by one-dimensional piecewise affine elements in order to implement the computational homogenization approach. Two discretisation schemes along the specimen height are adopted, as sketched in Fig. 8. The first one is uniform, with element size for scale ratios larger than 10, whereas 10 elements are used for all scale ratios smaller than 10 (cf. Fig. 8a). Such a discretisation is certainly excessively fine, and hence serves as a reference solution. The second discretisation, referred to as coarse, uses prior knowledge on the thickness and location of the two boundary layers, while keeping the size of the smallest macroscopic element to be 2 (except for small odd scale ratios), see Fig. 8b. Finally, in Fig. 8c, the RVE discretisation of the 2 × 2 periodic cell Ω m is shown, which makes use of quadratic isoparametric triangular elements. This discretisation has the same typical element size as the one used to obtain the DNS results discussed in Section 2.2. All constraints for w, i.e. periodicity in Eq. (38) and orthogonality w, ϕ 1 2 = 0, are treated as classical equality constraints. All microfluctuation fields w are condensed out for each macroscopic integration point i g by means of an equality constrained Newton-based minimization algorithm, cf. e.g. Bonnans et al. (2006) or Nocedal and Wright (2006), whereas the macroscopic governing equations are solved using again a quasi-Newton solver.
Computed results in terms of the coarse quantities v 0 and v 1 are compared against the DNS data u and σ in Fig. 9. An adequate approximation is now achieved in both the linear and in the post-bifurcation regimes. The boundary layers are captured properly, even for the coarse discretisation. In the linear regime, σ = 0 again as a consequence of the discrepancy in its definition compared to that of v 1 . The homogenized field v 1 does vanish here as expected, but unlike in the closed-form homogenized solutions of Section 4. The computed responses of the two uniform discretisations corresponding to scale ratios L/ = 5 and 12 and an applied strain of u D /L = 0.075 are further shown with the corresponding deformed RVEs in a couple of integration points in Fig. 10. At the centre, where v 1 is constant, the cells deform practically periodically, according to the pattern given by ϕ 1 . But closer to the boundaries v 1 approaches zero and hence the pattern becomes less prominent. Since v 1 varies significantly on the size of the RVE, 2 , the RVE deforms non-periodically.
Finally, the homogenized stresses are reported in Fig. 11, where a significant improvement compared to the closed-form homogenization approach presented in Section 4 can be observed. In the linear regime, the modulation function v 1 = 0, whereas w is periodic on top of an underlying affine deformation v 0 . This effectively results in an FE 2 type of approximation, with a relative error below 2 %, and practically no size effect, cf. Fig. 11a. In the post-bifurcation regime, v 1 = 0, and individual RVEs buckle in a synchronized manner coordinated through the term v 1 ϕ 1 . The microfluctuation field w still compensates any local fluctuations on top of the approximate quantity v 0 + v 1 ϕ 1 , significantly relaxing the system. This results in nominal stresses (P 1 ) 22 that are within 10 % of relative error in the case of the regular macroscopic discretisation, which is systematically more compliant compared to the DNS results. Note further that the largest error occurs for the smallest scale ratio considered, i.e. L/ = 4. The coarse discretisation is, on the other hand, somewhat stiffer, and provides in this particular case a slightly more accurate result, cf. Fig. 11b.

Summary and Conclusions
This paper presented a micromorphic homogenization framework, capable of capturing long-range correlated fluctuation fields emerging in elastomeric mechanical metamaterials due to local buckling of their microstructure. The ensemble averaging based homogenization framework, the emergent Euler-Lagrange equations and the numerical approximation schemes have been elaborated and commented. Results have been compared against averaged full-scale numerical simulations reported by Ameen et al. (2018).
The key developments and results can be summarized as follows: 1. A kinematic decomposition of the underlying displacement field considered for an entire family of translated microstructures has been proposed. 2. The proposed pattern-based decomposition, upon making use of the variational homogenization framework and ensemble averaging, has provided an effective micromorphic continuum enriched with the magnitude of the underlying long-range correlated fluctuation field. The additional micromorphic field establishes kinematic interactions between individual cells, which is essential for accurate homogenization and capturing quantitative size effects. 3. A computational homogenization scheme has been developed, making use of assumptions on the local energy density, smoothness properties of the effective kinematic fields, and the periodicity of the microfluctuation field. This allowed for an efficient and accurate solution of the macroscopic governing equations. 4. The key aspects of the implementation of the proposed framework have been discussed, focusing in particular on the effective fields v 0 and v 1 , along with the appropriate boundary conditions and orthogonality constraints. 5. Using a representative example, the ability of the proposed methodology to capture accurate solutions has been demonstrated in terms of kinematics as well as stress quantities. Homogenized stresses were always within 10 % relative error compared to the direct numerical simulations, capturing relevant size effects.
The proposed methodology, although presented on a specific case, can easily be extended to more complex systems, containing multiple long-range correlated fluctuation fields and exhibiting complex size effects.