A Newton Solver for Micromorphic Computational Homogenization Enabling Multiscale Buckling Analysis of Pattern-Transforming Metamaterials

Mechanical metamaterials feature engineered microstructures designed to exhibit exotic, and often counter-intuitive, effective behaviour. Such a behaviour is often achieved through instability-induced transformations of the underlying periodic microstructure into one or multiple patterning modes. Due to a strong kinematic coupling of individual repeating microstructural cells, non-local behaviour and size effects emerge, which cannot easily be captured by classical homogenization schemes. In addition, the individual patterning modes can mutually interact in space as well as in time, while at the engineering scale the entire structure can buckle globally. For efficient numerical macroscale predictions, a micromorphic computational homogenization scheme has recently been developed. Although this framework is in principle capable of accounting for spatial and temporal interactions between individual patterning modes, its implementation relied on a gradient-based quasi-Newton solution technique. This solver is suboptimal because (i) it has sub-quadratic convergence, and (ii) the absence of Hessians does not allow for proper bifurcation analyses. Given that mechanical metamaterials often rely on controlled instabilities, these limitations are serious. To address them, a full Newton method is provided in detail in this paper. The construction of the macroscopic tangent operator is not straightforward due to specific model assumptions on the decomposition of the underlying displacement field pertinent to the micromorphic framework, involving orthogonality constraints. Analytical expressions for the first and second variation of the total potential energy are given, and the complete algorithm is listed. The developed methodology is demonstrated with two examples in which a competition between local and global buckling exists and where multiple patterning modes emerge.


Introduction
Acting like a carefully engineered structure, rather than a standard bulk material, is a common characteristic of mechanical metamaterials. Recent advances in 3D printing and additive manufacturing enable the production of such structures on a relatively small scale, allowing to treat them as a homogeneous medium. Metamaterials are typically designed to exhibit an exotic behaviour which cannot be found in nature, such as a negative compressibility (Nicolaou and Motter, 2012), negative Poisson's ratio (Kolken and Zadpoor, 2017), or a high stiffness with an ultra low density (Zheng et al., 2014). In this contribution, we focus on elastomeric mechanical metamaterials, which under compression exhibit microstructural buckling resulting in a pattern transformation. Such a transformation induces an abrupt change in effective properties including Young's modulus and Poisson's ratio, with envisioned applications in, e.g., soft robotics (see Yang et al., 2015;Mark et al., 2016;Mirzaali et al., 2018).
Because elastomeric mechanical metamaterials rely mostly on local instabilities in their microstructural morphology, large deformations, rotations, and strains occur. In particular, the microstructure undergoes a pattern transformation, due to coordinated buckling of the underlying microstructure, resulting in a strongly non-local behaviour. If the specimen is restricted, e.g. by applied essential boundary conditions, the expected pattern cannot fully develop, and in the vicinity of the restriction the so-called boundary layers are formed. Because of pattern restriction such boundary layers generally behave stiffer compared to the bulk of the (meta)material, which may significantly influence the overall response even at the engineering scale. Such a configuration is depicted in Fig. 1a, in which an example of an elastomeric metamaterial beam subjected to compressive load is shown. The buckled pattern vanishes close to the two vertical boundaries, resulting in a stiffening effect. The extent to which the boundary layers influence the effective mechanical behaviour depends on the ratio of their thickness and the overall size of the specimen, or more generally on the scale ratio defined as the ratio between the overall size of the specimen H relative to the typical size of the microstructural features , i.e. H/ .
For predictive modelling of engineering-scale applications, it is important to accurately yet efficiently predict the overall mechanical response of the structure. To this end, homog-(a) local buckling (b) global buckling u 2 /u D Figure 1: An elastomeric mechanical metamaterial containing circular holes in a square packing exhibiting a pattern transformation under a compressive load induced by clamping the right-hand side of the specimen while prescribing a horizontal displacement u D at the left-hand side. Restricted evolution of the patterning near the edges results in stiff boundary layers. Depending on the slenderness of the specimen, local (a) or global (b) buckling occurs. In the case of global buckling, the patterning mode (i.e. local buckling) is localized only in the compressive regions of the specimen. The colour indicates the pointwise magnitude of the displacement field u relative to u D . enization techniques are employed, replacing the complex microstructural behaviour with an equivalent continuum model. However, due to the non-locality, patterning, and buckling of the microstructure, homogenization of mechanical metamaterials presents a difficult challenge. First-order computational homogenization, outlined e.g. by Kouznetsova et al. (2001), can significantly reduce computing time compared to full scale Direct Numerical Simulations (DNS). However, its inherent assumptions on locality and scale separation prevent it from predicting any size effects in the microstructure. Ameen et al. (2018) showed that as a consequence relative errors up to 40% can be induced in terms of force quantities by the first-order method in the post-bifurcation regime for small scale ratios. With an increasing scale ratio, the accuracy of the predicted overall behaviour typically improves, with an exact match for H/ → ∞. The second-order computational homogenization (Kouznetsova et al., 2004), which incorporates the gradient of the macroscopic deformation gradient in the micro-to-macro transition, permits to reflect non-locality and ensuing size-effects. The effective behaviour captured by this method coincides satisfactorily with DNS results even for low scale ratios, albeit at the cost of additional complexity stemming from a higher-order continuum formulation at the macroscale level; see Sperling et al. (2020) for more details. Recently, Rokoš et al. (2019) proposed a micromorphic computational homogenization framework, specifically designed to predict the effective behaviour of mechanical metamaterials by decomposing the displacement field into three components: (i) a smooth, mean displacement field, (ii) a spatially correlated microfluctuation field, and (iii) an uncorrelated, local microfluctuation field. This decomposition ensures an adequate performance and accuracy by introducing prior knowledge on the patterning fluctuation. See Sperling et al. (2020) and Rokoš et al. (2020a) for more details.
A major practical limitation of the previously reported implementations of the micromorphic computational homogenization framework is that a quasi-Newton solution method has been employed. This method is not as efficient as a full Newton scheme, it has shown to be quite sensitive to the initialization of the (equilibrium) iteration process, and in particular to the perturbations applied to trigger buckling and, related to the latter point, it does not allow for a proper bifurcation analysis. For elastomeric metamaterials in particular, such a buckling analysis is essential for the prediction of the local patterning (resulting in a transition in the effective mechanical properties, see Fig. 1a) or global buckling (indicating a potential failure of the entire structure, see Fig. 1b). Without a proper Newton algorithm, such phenomena can hardly be captured in an accurate and reliable way.
The macroscopic instability at the level of a material point induced by a microscopic bifurcation has been studied by Saiki et al. (2002), whereas Wadee and Farsi (2015) have investigated geometrical effects on the buckling behaviour of cellular structures in which a transition between local and global buckling is observed as a function of the specimen slenderness. Experimentally and numerically, Niknam and Akbarzadeh (2018) have compared in-plane and out-of-plane buckling of various types and sizes of architected cellular structures. Specimens with a hexagonal honeycomb microstructure, which exhibit multiple buckling patterns under different compressive biaxiality ratios, were studied by Ohno et al. (2002a). Rokoš et al. (2020a) demonstrated the same multi-pattern character for hexagonally stacked cells with circular holes using the micromorphic homogenization framework. Obtaining correct patterning of the microstructure required, nevertheless, intervention of the user based on insight in the mechanics of the system, precisely because no reliable solver was available to tackle the buckling.
The main goal of this paper is to derive the tangent operator for the micromorphic computational homogenization framework, enabling a more efficient and robust solution procedure using a full Newton algorithm and allowing for bifurcation analyses (Miehe and Bayreuther, 2007). Unlike the derivation of the Hessians (i.e., the macroscopic tangents or stiffnesses) for first-order computational homogenization-see detailed explanation in (Miehe and Koch, 2002) or (Miehe, 2003)-, the Hessians for the micromorphic scheme require a non-trivial extension. Additional orthogonality constraints acting within each Representative Volume Element (RVE) need to be enforced in order to guarantee uniqueness of the adopted kinematic decomposition. Using variational calculus, the first variation of the averaged energy resulting in the microscopic and macroscopic governing equations is derived, followed by the second variation from which the micro-, and coupling macro-Hessians can be obtained. Following Rokoš et al. (2020a), the formulation involves an arbitrary number of patterning modes, and introduces a slight reformulation of the orthogonality constraints with respect to gradients of individual modes as compared to the original framework (Rokoš et al., 2019) in order to eliminate spurious oscillations observed in the resulting micromorphic fields. Employing a standard Finite Element (FE) discretization, the associated internal forces and stiffnesses are constructed, from which the local microfluctuation fields are condensed out, yielding a macroscopic Newton algorithm. We illustrate the performance of the method with two examples, one focusing on local versus global buckling of a metamaterial column, and one on local patterning of a hexagonally-voided microstructure.
The remainder of this paper is organized as follows. After recalling the kinematic decomposition pertinent to the micromorphic framework, Section 2 details the derivation of the first and second variation of the ensemble averaged energy, resulting in both macroand microscopic governing equations accompanied by the relevant macro-and microscopic Hessians. Employing standard FE procedures, Section 3 describes the discretization of the governing equations at both scales and addresses the bifurcation analysis. Section 4 illustrates the developed methodology with two examples: (i) a compressed metamaterial column with a varying slenderness ratio in which a competition between local and global buckling exists, and (ii) a microstructure with hexagonally-stacked holes subjected to biaxial compressive loading exhibiting three distinct pattern transformations. Finally, the summary and conclusions follow in Section 5.
Throughout the paper, the following notational conventions are used scalars a, vectors a, position vector in the reference configuration X = X 1 e 1 + X 2 e 2 , fourth-order tensors 4 A = A ijkl e i e j e k e l , matrices A and column matrices a, derivatives of scalar functions with respect to second-order tensors where Einstein's summation convention is adopted on repeated indices i, j, k, l, and e i , i = 1, 2, denote the basis vectors of a two-dimensional Cartesian coordinate frame.

Kinematic Decomposition
The micromorphic computational homogenization framework (Rokoš et al., 2019), depicted schematically in Fig. 2, relies on the decomposition of the kinematic field u into the mean effective displacement v 0 , long range correlated fluctuation components v i ϕ i , i = 1, . . . , n, and the remaining local microfluctuation field w, i.e. (1) The vector field ϕ i corresponds to the i-th patterning mode of the underlying microstructure, whereas the scalar field v i regulates, spatially and in time, its magnitude. Unlike the original formulation (Rokoš et al., 2019), who considered only one such mode, we consider here an arbitrary number of modes, n; see also Rokoš et al. (2020a). Because in general it may not be possible to control the positioning of the microstructure relative to the specimen's boundary, all possible microstructural translations should be taken into account via ensemble averaging (cf. Ameen et al., 2018). The micromorphic scheme avoids this costly procedure by approximating the mechanical state of a point in a translated microstructure by evaluating the mechanical state of a microstructurally equivalent point in the reference microstructure. In addition, a separation of scales into a macroscopic position vector X and a microscopic position vector X m is introduced, assuming the fields v 0 and v i to vary slowly over a close vicinity of each macroscopic point X spanned by a microscopic Representative Volume Element (RVE) with a domain Ω m . Consequently, the microfluctuation field w is computed only locally over each RVE, and is independent for RVEs associated with distinct macroscopic points X; they communicate only by means of the macroscopic fields v 0 and v i . Therefore, v 0 and v i become functions of the macroscopic position only, whereas ϕ i and w are functions of the microscopic as well as the macroscopic position. Because the patterning mode ϕ i is the same for each macroscopic point, it eventually depends on the microscopic position vector X m only. Using a first-order Taylor expansion and the above considerations, the decomposition of Eq. (1) can be approximated as For more details on the decomposition (2), the reader is referred to Rokoš et al. (2019)albeit for a single mode ϕ 1 . The fields v 0 , v i , and w, are unknown and need to be solved for, while the individual patterning modes ϕ i are characteristic for the underlying microstructural morphology and are assumed to be known a priori, computed either from a Bloch-type analysis (Bertoldi et al., 2008), estimated analytically from full-scale numerical simulations (Rokoš et al., 2019), or identified experimentally (Maraghechi et al., 2020). Because the patterning modes ϕ i are defined with respect to the reference microscopic configuration X m , the effect of the macroscopic rotations needs to be factored out from the macroscopic deformation gradients, where R M is the macroscopic rotation tensor and U M is the macroscopic stretch tensor. This effectively means that the term ∇ v 0 needs to be replaced with U M − I in Eq.
(2), and that all effective stress and stiffness quantities need to be rotated back accordingly, see e.g. (Kunc and Fritzen, 2019, Section 3.2) for more details. Such a distinction is, however, omitted hereafter to simplify all derivations, and can even be neglected in the limit of small rotations as is the case for both examples shown in Section 4. Note also that the additional micromorphic fields v i in Eqs.
(1) and (2) relate directly to the displacement field u. Such a setup contrasts with standard micromorphic formulations, in which the microdeformation χ typically relates to the gradient of u, see e.g. (Forest and Trinh, 2011, Eqs. (29)-(32)). The decomposition of Eq.
(2) further suggests that when all micromorphic Microscale boundary value problem, Eq. (13) Figure 2: A schematic representation of the micromorphic computational homogenization framework. The macroscopic displacement gradient ∇ v 0 , micromorphic fields v i , and their spatial gradients ∇v i at a macroscopic point X are sampled and sent to the microscale, where a boundary value problem for the microfluctuation field w is solved. Based on the solution w, the macroscopic stresses, Θ, Π i , Λ i , and stiffnesses are computed and passed back to the macroscale, where the macroscopic boundary value problem is assembled and solved.
The uniqueness of the decomposition (2) within an RVE is guaranteed by introducing the following additional orthogonality conditions: w( X, X m ) · ϕ i ( X m ) Ωm = 0, i = 1, . . . , n, Recall that throughout this contribution the angle brackets indicate the integration over the domain specified by the subscript. The first condition (3) requires zero mean of w over Ω m , effectively eliminating rigid body translations. The second condition establishes the uniqueness in terms of the patterning modes ϕ i themselves, because, upon assuming a homogeneous state with ∇v i = 0, for instance, the patterning can be equally well represented by the product of the micromorphic field and the patterning mode, i.e. v i ϕ i , or by the microfluctuation field w while keeping v i = 0. The last orthogonality condition has not been previously introduced in the original formulation of Rokoš et al. (2019). It follows, however, from the decomposition outlined in Eg.
(2), and eliminates non-uniqueness issues related to linearly varying magnitudes of the patterning fields ϕ i , which may cause spurious macroscopic oscillations for microstructures with a hexagonal stacking of holes. This condition acts mainly as a stabilization condition and does not significantly affect the overall mechanical response. Following the same reasoning as in first-order computational homogenization, the remaining orthogonality condition, which might arise from the non-uniqueness related to the ansatz (2) (i.e. orthogonality with respect to X m · ∇ v 0 ( X), or w( X, X m ) X m Ωm = 0 for an arbitrary ∇ v 0 ), is accounted for differently. As a well-accepted modelling choice, which provides accurate results in the first-as well as second-order computational homogenization schemes, a periodicity constraint on w is adopted ensuring this orthogonality, i.e.
where w( X, X m ) = w( X, ∂Ω + m ) − w( X, ∂Ω − m ) denotes the jump of the field w( X, X m ) on the RVE boundary split into two parts: ∂Ω m = ∂Ω + m ∪ ∂Ω − m . As an example, a 2 × 2 RVE of a square stacking of holes is schematically shown in Fig. 3 The same approach is used for polygons with multiple edgessee e.g. ahead to Fig. 10b for the case of a hexagonal RVE. The periodicity constraint in combination with the condition of Eq. (3) also eliminates all rigid body modes of w. Although periodic boundary conditions for w (and thus antiperiodic RVE boundary tractions) are adopted, such a choice might not necessarily result in the optimal performance of the proposed micromorphic scheme. See e.g. Forest and Trinh (2011), where antiperiodic microfluctuation fields w with periodic RVE boundary tractions have been observed. Yet, it will be demonstrated in the results Section 4 that the periodicity constraints (6) provide an adequate accuracy here.

Potential Energy
Because we are restricting ourselves to hyperelastic materials, the unknown parts v 0 ( X), v i ( X), and w( X, X m ) of the solution u( X, X m ) according to (2) can be found by minimizing the total potential energy of the system while accounting for the constraints introduced above, i.e.
where the displacement field u( X, X m ) is, through Eq. (2), also a function of all unknown macroscopic, v 0 ( X), v 1 ( X), . . . , v n ( X), and microscopic, w( X, X m ), quantities in addition to the two spatial variables X and X m , and where the fields µ( X), ν( X) = [ν 1 ( X), . . . , ν n ( X)], and η( X) = [ η 1 ( X), . . . , η n ( X)] collect the Lagrange multipliers pertinent to rigid body modes and orthogonality constraints with respect to individual patterning modes as their components. While the bulk constraints of Eqs.
(3)-(5) are considered for each RVE separately, i.e. µ( X), ν( X), and η( X) are a function of the macroscopic position vector X only, the periodicity constraint (6) is a continuous function on the boundary of each RVE and therefore, it is also a function of the microscopic position vector X m , i.e. λ( X, X m ). In Eq. (7), the Lagrangian L = E + C consists of the total potential energy E where the minus sign in front of the last constraint term is adopted for consistency reasons (Miehe and Bayreuther, 2007, Section 2.2.2), and the remaining terms take a plus sign to obtain positive quantities on the right-hand side of the resulting governing equation (see Eq. (13) below). Note that this slightly deviates from Rokoš et al. (2019), where the constraint term C was incomplete. A hyperelastic constitutive law specified through the energy density function Ψ( X m , F ) is adopted for the description of the material behaviour; F ( u( X, X m )) = I + ( ∇ m u( X, X m )) T is the deformation gradient and ∇ m = e i ∂/∂X m,i is the microscopic gradient operator defined in the reference configuration. Note that the explicit dependency of F on u has been dropped in Eq. (8), and will be omitted for brevity hereafter as well.

First Variation and Governing Equations
A minimizer of the total potential energy can be found by taking the Gâteaux derivative (i.e. the first variation) of the Lagrangian L and requiring it to vanish: where the local first Piola-Kirchhoff stress tensor P ( X m , F ) is defined as Making use of the decomposition introduced in Eq. (2), ∇ m δ u reads Upon substituting this expansion, application of the divergence theorem, and rearrangement of individual terms, Eq. (10) leads to a set of microscopic and macroscopic balance equations. At the microscale, only the variations of the microfluctuation field w( X, X m ) and its related Lagrange multipliers µ( X), ν i ( X), η i ( X), and λ( X, X m ) matter. Since Eq. (10) must hold for arbitrary variations of these quantities, the following set of microscale balance equations results δ w : on ∂Ω ± m , δ λ : periodicity constraint for w, Eq. (6), δ µ, δν, δη : kinematic constraints for w, Eqs. (3)-(5).
As a consequence of the constraint term introduced in Eq. (9), the right hand side of the first governing Eq. (13) involves Lagrange multipliers acting as body forces inside each RVE, anti-periodic condition for RVE boundary tractions, and an additional set of orthogonality constraints. These terms were not present in the original formulation, where the constraint equations were enforced differently. Although the body forces µ, v i ϕ i , and η i · ( ϕ i X m ), are directly associated with the orthogonality constrains of Eqs. (3)-(5) in the form of Lagrange multipliers, they can be introduced also directly at the level of governing equations, see e.g. (Yvonnet et al., 2020, Eqs. (8) and (9)). At the macroscale, only the slowly varying fields v 0 ( X) and v i ( X) are relevant, and their governing equations read with the following definitions of macroscopic stress-like quantities: The dependence of the homogenized stress quantities on the macroscopic position X originates from the dependence of the deformation gradient F ( u( X, X m )) on the macroscopic position through all smooth fields v 0 ( X) and v i ( X), cf. also Eq. (2).

Second Variation
The second variation of the Lagrangian L reads where the spatial dependencies on X and X m have been dropped for brevity, and where the local stiffness tensor 4 D m ( X m , F ) has been introduced as Positive definiteness of the modified Lagrangian, from which all constraint variables-i.e. Lagrange multipliers-have been condensed out, reflects stability of the combined multiscale system including micro-as well as macro-quantities. Upon condensing out also the microfluctuation field w, stability of the macroscopic system can be assessed. The microfluctuation field as well as all Lagrange multipliers can be eliminated for each macroscopic point X by static condensation. This is done by means of the Schur complement after the FE numerical discretization, as described in Section 3 below, whereas stability of the macroscopic system is detailed in Section 3.3. Substituting Eq. (12) into the first expression on the right hand side of Eq. (18), one obtains a set of p(p + 1)/2, p = n + 2, coupled specific stiffness quantities for the macroscopic, v 0 , v i , i = 1, . . . , n, and microscopic, w, fields. In particular, the individual terms corresponding to the macroscopic quantities read as where the subscript 0 relates to ∇δ v 0 , i B to ∇δv i , and i N to δv i terms, respectively, with i = 1, . . . , n, and j = 1, . . . , n, which are essential for numerical solution and stability assessment of the macroscopic system.

Numerical Implementation
At every macroscopic integration point, the macroscopic quantities ∇ v 0 , v i , and ∇v i , are sampled and passed down to the microscale, cf. Fig. 2, where the modes ϕ i and microfluctuation field w are defined. Taking into account the constraints (3)-(6), the microfluctuation field w is computed by solving the microscale boundary value problem defined in Eq. (13). Knowing w, the homogenized macroscopic stresses and stiffnesses are computed following Eqs. (15)-(18) and (20)-(25). All of these quantities are required for the solution of the macroscopic boundary value problem of Eq. (14) using the standard Newton algorithm, leading to a fast quadratic convergence and allowing for a bifurcation analysis. Although multiple approaches can be adopted for the solution of the multiscale problem, see e.g. (Okada et al., 2010), here we adopt the condensation method. The discretization and numerical solution of the microscopic problem is presented in Section 3.1. The macroscopic problem is subsequently elaborated upon in Section 3.2. The bifurcation analysis is detailed in Section 3.3, including a nested algorithmic scheme. A Matlab implementation of the presented Micro-Morphic homogenization for Multiscale Metamaterials framework (mm4mm) is available at git repository mm4mm.

Microscopic Problem
Using standard FE procedures, the microfluctuation field w and its gradient ∇ w are expressed over each microscopic element e within an RVE in terms of the shape functions and nodal values as where w e is a column of element nodal values of the w field, collected for all n e elements in a column matrix w, N e w is a matrix of the corresponding shape functions, and B e w a matrix of the shape function derivatives. Corresponding variations δ w and ∇δ w are discretized in the same way. The patterning modes are expressed similarly as where ϕ e i is a column storing element nodal values of ϕ i , collected over all elements in a column matrix ϕ i . Although analytical expressions for the patterning modes ϕ i are provided below in Eqs. (56) and (62), we opted here for a discretized version for convenience and generality, since alternative definitions of patterning modes may be more appropriate, e.g. based on linearized buckling analysis or true deformed shapes obtained numerically. The discretized version of the orthogonality constraint (3) takes the form where w Ω e m is the symmetric Gramian matrix of microscopic shape functions, and 1 k is a column matrix with ones at the positions of Degrees Of Freedom (DOFs) corresponding to the k-th component (i.e. either horizontal or vertical component) of w.
A denotes the standard FE assembly operator, and the integration over each microcopic element volume Ω e m is performed using a standard Gauss integration rule, for brevity not expressed explicitly as a sum. Analogously, the discrete forms of the scalar constraints (4) whereas the set of vector constraints (5) is discretized as where q i,k stores ( ϕ i · e 1 )X k and ( ϕ i · e 2 )X k at the positions of the DOFs corresponding to the first and second component. The periodicity constraint of Eq. (6) is expressed with the help of the link topology matrix C pbc , described e.g. in (Miehe and Bayreuther, 2007), as Finally, the microscopic governing equation (13) is solved iteratively using the standard Newton method (see Bonnans et al., 2006, Section 14) for the linear system or in a compact form as In Eqs. (32) and (33), is the microscopic stiffness matrix, f w is the column of microscopic nodal internal forces f w = A ne e=1 (B e w ) T P Ω e m and unbalanced equality constraints, where P is a column storing the components of the stress tensor P , and dw is a column storing the iterative change of the microscopic fluctuation field and current iterative values of Lagrange multipliers. Note that, in analogy to η and ν in Eq. (7), µ stores individual components of µ. An asterisk indicates that the matrix accounts for all the equality constraints and thus also has the corresponding extra entries. To evaluate the vector of current internal forces f w and stiffness matrix D ww , the nodal values of the entire displacement field u of Eq. (2) need to be constructed over the entire RVE from the knowledge of the current state of the macroscopic quantities ∇ v 0 , v i , ∇v i , patterning modes ϕ i , i = 1, . . . , n, and iterative state of the microfluctuation field w. Because the primary buckling modes are captured by the patterning fields ϕ i , a microscopic bifurcation analysis and stability control analogous to Section 3.3 below is not required.

Macroscopic Problem
The macroscopic fields v 0 and v i are discretized within each macroscopic element E as where N E • and B E • are macroscopic element shape functions and v E • the corresponding column matrices of DOFs, collected globally for all n E elements in column matrices v 0 and v i ; the same forms and expressions are used to discretize their variations. The internal element forces are then obtained as a sum over n g macroscopic Gauss quadrature points, expressed explicitly as and assembled over all macroscopic elements to form the global internal force vector f M = A n E E=1 f E M . In Eqs. (37) and (38), w ig are integration weights with the corresponding Jacobians J ig , and the column matrices Θ ig , Λ In order to condense out the effect of the constrained microfluctuation field w and all Lagrange multipliers, the following monolithic incremental system of equations is assembled at each macroscopic quadrature point i g of each element E, including both macroscopic as well as microscopic quantities (superscripts E and i g are dropped in Eqs. (39)-(49) for brevity):  In Eq. (39), r i , i = 0, . . . , n, are the residuals reflecting the fact that equilibrium is not satisfied at the level of an integration point, but only for the entire assembly over all quadrature points of all elements. The last row is zero, since the microscopic system of Eq. (33) has been equilibrated, implying also that f w = 0. The specific stiffness for the macroscopic quantities can be then obtained as a Schur complement via static condensation: Note that whereas for a given element E an n g microfluctuation fields w ,E,ig are computed and condensed out, only one set of macroscopic DOFs for the coarse fields v E 0 and v E i pertinent to that element are involved. The asterisk superscript to the specific stiffness components in Eqs. (39) and (40) again indicates that these matrices have extra zero entries corresponding to the Lagrange multipliers. The individual sub-matrices are defined as where D 00 , , are matrix representations of the specific stiffnesses defined in Eqs. (20)-(25), obtained upon RVE discretization simply as volume integrals using a standard Gauss integration rule. For instance, D is a 4 × 4 matrix representation of the microscopic fourth-order stiffness tensor 4 D, etc. In addition to the expressions (41)-(46), there are 3n + 1 coupled specific stiffnesses related to the microscopic variation ∇δ w and one of the macroscopic variations ∇δ v 0 , δv i , or ∇δv i , and one microscopic stiffness quantity that relates twice to ∇δ w. Because ∇δ w depends on the microscopic position X m , over which the integral in Eq. (18) is carried out, these specific stiffnesses can be written only upon RVE discretization, yielding where the integration over Ω e m is again carried out numerically. The stiffness matrix of the entire macroscopic element is obtained by summing the contributions from all quadrature points, which are eventually assembled into a global stiffness matrix K M = A n E E=1 K E M . The resulting macroscopic system is again solved using the standard Newton method with an incremental system of linear equations where f ext denotes a column of externally applied forces (acting only on v 0 ), and is an iterative increment of the global macroscopic quantities.

Bifurcation analysis
Following Miehe and Koch (2002), an equilibrated configuration v M of a system is considered to be stable if the energy of this state is lower than the energy associated with a state obtained by adding a small kinematically admissible perturbation δv M to the equilibrated configuration. That is, if where the second-order Taylor series expansion of the total energy has been used. Notice that E corresponds to the total potential energy of the entire system (Eq. (8)) from which microfluctuation fields and Lagrange multipliers, w , associated with all macroscopic Gauss points have been condensed out, and that the first-order term vanishes as a result of equilibrium. The condition (53) is equivalent to the requirement of positive definiteness of K M , i.e. to the requirement that all eigenvalues of K M are positive. If the lowest eigenvalue α 1 is non-positive, the associated configuration is unstable and the corresponding eigenvector ψ 1 determines the buckling mode. The equilibrated solution v M of the current increment is then perturbed with the eigenvector ψ 1 multiplied by a small perturbation factor τ > 0, and the system is equilibrated again. The factor τ is increased until a stable equilibrium is reached, i.e. until a possible energy barrier between the current unstable and a stable buckled configuration is overcome, and at the same time until the lowest eigenvalue α 1 of the updated macroscopic stiffness matrix K M does not become positive. If the system fails to find a stable equilibrium even for large τ , the previous increment is halved to decrease the energy barrier, and the entire procedure is repeated. An outline of the overall micromorphic computational homogenization scheme for multiple modes, including stability control, is given in Algorithm 1.

Numerical Examples
In this section, predictions made using the micromorphic computational homogenization scheme, introduced in Sections 2 and 3, are compared against Direct Numerical Simulations (DNS) for two examples. The first example represents a metamaterial column composed of a square stacking of holes subjected to compression, whereas the second example considers a specimen with a hexagonal stacking of holes subjected to a uniform compressive loading with various biaxiality ratios, buckling locally into one of multiple possible patterns.
The constitutive behaviour of the elastomer base material is modelled by a hyperelastic law with the following energy density where F = I + ( ∇ u) T is the deformation gradient tensor, where the gradient operator ∇ is defined with respect to the reference configuration, J = det F , and I 1 = tr C is the first invariant of the right Cauchy-Green deformation tensor C = F T · F . The values of the constitutive parameters employed, listed in Tab. 1, are based on the experimental characterization of Bertoldi et al. (2008). The smallest RVE domains of the size 2 are adopted in both examples for the micromorphic homogenization scheme, as shown in Figs. 4b and 10b below, which are large enough to accommodate the longest microstructural patterning modes (see e.g. Bertoldi et al. 2008 for the square and Ohno et al. 2002a for hexagonal stacking of holes). Although the chosen RVE size is sufficiently large to accommodate microstructural buckling, because of the periodicity assumption on the microfluctuation fields w (recall Eq. (6) and the discussion therein), choosing larger RVE domains may still slightly affect obtained results, especially for a vanishing separation of scales. M is non-positive, perturb the system with corresponding eigenvector τ ψ 1 , and equilibrate iteratively for an increasing perturbation factor τ until the system becomes stable. Then continue to (i) for k = k + 1. If perturbation fails, halve the load increment and proceed to (i) with current k. 3: end for

Example 1: Local Versus Global Buckling
The first example analyses a finite column of width W and height H with a microstructure consisting of a square stacking of unit cells with edge size = 9.97 mm and circular holes of diameter d = 8.67 mm, cf. Fig. 4a. The bottom and top edges of the specimen domain are displaced by ±u/2 e 2 to induce 10% overall compressive strain, defined as u/H. Depending on the slenderness ratio H/W , a competition between microstructural buckling (pattern transformation) and macrostructural buckling of the structure is expected. A similar example has been investigated numerically as well as experimentally by Coulais et al. (2015).
For the DNS solutions, the entire domain is discretized using isoparametric quadratic triangular elements of typical size h m = /10 with three Gauss integration points, as shown in Fig. 4b. For this case, a single microstructural realization adequately represents the ensemble averaged DNS solution. This is shown in Fig. 5a, where nominal stress-strain diagrams corresponding to 100 microstructural translations (all possible combinations of 10 translation steps in horizontal and 10 in vertical direction, covering together one period of the microstructure × ) are shown for a 6 × 30 specimen. The overall response is initially linear until the first bifurcation point is reached, upon which a local patterning emerges and the specimen's stiffness drops close to zero. Further increasing the compressive strain leads to the second bifurcation, corresponding to global buckling of the specimen, upon which the overall stiffness becomes negative. The maximum and minimum envelopes of all realizations deviate less than 6% from the corresponding mean, suggesting that the reference microstructure is acceptable for the representation of the effective response.
Only one local patterning mode emerges for the adopted microstructural morphology, i.e. n = 1, see Figs. 1-2 and Bertoldi et al. (2008), approximated analytically as (see Rokoš et al., 2019, Eq. (7)) were C 1 is a normalization constant ensuring that 1 |Q| Q ϕ 1 ( X) d X = 1, and where Q is the 2 × 2 periodic cell of Fig. 3. For the RVE discretization, the same type and density of elements is used as for the DNS system (Fig. 4b). To provide sufficient kinematic freedom to the macroscopic micromorphic system, a mesh convergence study is performed. A uniform macroscopic mesh of quadratic isoparametric triangular elements with three Gauss points is considered with characteristic element sizes h M ∈ {1, . . . , 30} . The same mesh is considered for both macroscopic fields v 0 and v 1 . An example of a particular mesh of an element size h M = 2 for a 4 × 8 specimen (for which local buckling is expected to occur) is shown in Fig. 4c. The obtained results in terms of nominal stress-strain diagrams are plotted in Fig. 5b. For element sizes h M ≤ 4 , the behaviour is similar to the DNS result of Fig. 5a; moreover, the results of element sizes h M = 2 and h M = 1 are indistinguishable. The element size h M = 2 is thus adopted in what follows, although from the deformed shape of Fig. 6d it may be clear that a locally refined mesh might be useful. For more details on appropriate choice of element types and associated integration rules see Rokoš et al. (2020b). Depending on the slenderness ratio H/W , two basic and mutually interacting deformation mechanisms occur, as shown in Fig. 6 by the deformed configurations for 6 ×14 and 6 ×34 specimens. The first mechanism is local patterning (Fig. 6a), emerging for low slenderness ratios upon reaching a critical compressive strain of approximately 3%. The cells fold in a typical pattern of alternating ellipsoidal holes and an auxetic effect is observed along with boundary layers where the local buckling is restricted. The second mechanism, occurring for higher slenderness ratios, is global buckling (Fig. 6c), which is triggered upon reaching the critical Euler buckling stress. The corresponding buckling strain can be estimated as In Figs. 6b and 6d, the micromorphic field normalized by its maximum value considered over space and a parametrization pseudo-time t, i.e. v 1 = v 1 / v 1 (t, X) ∞ , is shown in colour. Comparing Fig. 6a with 6b, and Fig. 6c with 6d, we conclude that the micromorphic homogenization scheme is capable of accurately reconstructing the overall kinematic response, correctly capturing the auxetic effect for lower slenderness ratios, and accurately indicating regions of localized patterning reflected by the magnitude of the micromorphic field v 1 for larger slenderness ratios. The pattering regions localize in the compressive parts of the bent domain, situated close to the specimen's centre and near the supports at both ends. The overall deformed shape for the large slenderness ratio, i.e. the v 0 field, is captured with good accuracy as well. Nominal stress-strain diagrams for W = 6 and slenderness ratios H/W ∈ [1, 30] are shown in Fig. 7. Here, the two mechanisms and their mutual interactions are visible more clearly. For the DNS results (Fig. 7a) and the applied strain range u/H ∈ [0, 0.1] the slenderness ratios up to H/W = 2.33 buckle only locally, ratios H/W ∈ [2.67, 6.67] buckle first locally and then globally, the ratio H/W = 7 buckles locally and globally at the same time, whereas ratios H/W ∈ [7.33, 30] buckle first globally and then locally. Bilinear stress-strain responses typically emerge, which exhibit softening in later stages due to the presence of secondary (local or global) buckling, see also Fig. 5a and the discussion therein. A close-up on the local versus global buckling intersection is shown in Fig. 8a, where mild snap-backs for slenderness ratios H/W ∈ [6.33, 8.33] can be observed. The micromorphic computational homogenization is capable of reconstructing the stability behaviour ( Fig. 7b) with an accuracy that decreases with decreasing scale ratio (i.e. larger errors are observed for smaller scale ratios). In particular, for H/W = 1 we see a large discrepancy in the post-bifurcation nominal stiffness (i.e. in the slope of the P 22 versus u/H curve), which corresponds to 25% of error relative to the initial pre-bifurcation stiffness (which is practically constant for all considered H/W ratios). With increasing slenderness ratio, however, the error drops rapidly down to 0.1%. For better clarity, the bifurcation curves corresponding to the DNS and micromorphic results are compared in Fig. 8b. Here the shapes as well as slopes of the DNS bifurcation curves (shown in blue) are captured accurately by the micromorphic scheme (shown in red), although micromorphic homogenization systematically overestimates the DNS results. The maximum relative error in terms of the critical buckling stress of the first instability is of the order of 12%, but not lower than 7% even for large scale ratios.
The buckling strains expressed as a function of the slenderness ratio H/W and corresponding to the first bifurcation points of Fig. 7 are shown in Fig. 9 for several specimen widths W ∈ {6, 8, 10} . Bertoldi et al. (2008) reported that the buckling strain of a single RVE corresponds to approximately 3% (shown as the black dash-dot line in Fig. 9), and hence this value is expected to be the theoretical local buckling strain. Both the DNS as well as micromorphic results attain this limit for the range of slenderness ratios H/W ∈ [1, 7], although for very small slenderness ratios a mild increase in the critical strain is observed. This effect is explained by the growing influence of the stiff boundary layers constraining the evolution of the microstructural patterns. Note that short columns may still buckle globally, upon further increase of the external load. For a slenderness ratio of approximately H/W = 7, local and global buckling occur simultaneously, whereas higher slenderness ratios converge asymptotically towards the theoretical bound of the global Euler buckling strain given by Eq. (57) (shown as the black dashed line). Both the DNS and the micromorphic scheme approach this limit from below, although the micromorphic results overestimate systematically the critical buckling strain obtained by the DNS. With increasing width of the specimen W , the size effects present for small slenderness ratios slowly decrease. For short columns the local buckling strain converges towards the theoretical bound 3%, whereas the global buckling strain shows little change. The overall relative error of the micromorphic scheme in terms of the buckling strain does not exceed 3% for local and 6% for global buckling. With increasing scale ratio, this error drops down to 0.5% for local and 3% for global buckling. Note that the theoretical global Euler buckling strain given by Eq. (57) significantly overestimates the DNS results due to a substantial amount of shear and changes triggered in the microstructure upon buckling (cf. Fig. 6c), which become important especially for intermediate slenderness ratios H/W ∈ (7, 20].

Example 2: Multiple Local Modes
The second example considers an infinite microstructure composed of a hexagonal stacking of holes, shown in Fig. 10a. The periodic cell, considered later as RVE for the micromorphic scheme, thus comprises two holes in each of the three directions along the hole centres, see Fig. 10b. As reported by Ohno et al. (2002a) for the case of hexagonal honeycombs, three different patterns can emerge under biaxial compression, depending on the biaxiality ratio where ε ii = F ii − 1 are the nominal compressive normal strains in the e i , i = 1, 2, directions, and F ij are the components of the overall deformation gradient tensor F , with F 12 = F 21 = 0. The three distinct patterns, depicted in Figs. 11a-11c, correspond to the following cases: (i) Pattern I, uniaxial or shear pattern denoted π 1 , occurs when the compressive load on the vertical cell walls is higher than the load on the other cell walls, i.e. γ > 1. The multiplicity of the bifurcation point corresponds to one. The displacement field leads to the formation of horizontal layers of holes sheared alternatingly to the right and to the left, see Fig. 11a. (ii) Pattern II, also called biaxial or butterfly-like pattern and denoted π 2 , emerges when the inclined cell walls at θ = ±30 • are compressed more than the vertical cell walls, i.e. γ < 1. In this case, the multiplicity of the bifurcation point is two and the pattern exhibits horizontal layers of holes buckled along the horizontal and vertical directions, see Fig. 11b. (iii) Pattern III, also referred to as the equi-biaxial or flower-like pattern π 3 , is observed when all three cell walls are subjected to an equal compressive load, i.e. γ = 1. The multiplicity of the bifurcation point equals three, and the displacement field corresponds to a virtually undeformed central hole, surrounded by ellipses, see Fig. 11c. Because the individual patterns π i are not mutually orthogonal, it is convenient for further treatment to introduce the so-called modes ϕ i , i = 1, 2, 3, (see Ohno et al., 2002b;Okumura et al., 2002;Rokoš et al., 2020a), which satisfy orthogonality. Linear combinations of these modes result in the previously introduced patterns as follows: The individual modes ϕ i , i = 1, 2, 3, correspond to the shear pattern I (Figs. 11a and 11d) developing perpendicular to each of the cell wall directions, i.e. at θ = 90 • and ±30 • , which can be expressed in an analytical form. The first mode reads (see Rokoš et al., 2020a, Eq. (3)) were C 1 is a normalization constant ensuring that 1 |Q| Q ϕ 1 ( X) d X = 1 for the periodic cell Q of Fig. 10b, whereas modes II and III are obtained by rotating ϕ 1 by ∓60 • , see Figs. 11e and 11f.
The example analysed in this section represents an infinite specimen, made of a hexagonal cellular structure with a hole diameter d = 1.28 mm and a centre-to-centre spacing = 1.386 mm, subjected to biaxial compression. For the micromorphic computational homogenization it is modelled with a 10 × 10 mm 2 periodic square domain discretized with eight identical quadratic triangular elements of size h m = 3.6 with a three-point Gauss integration rule, see Fig. 10c. Again, the same macroscopic discretization is used for all three micromorphic fields v i , i = 1, 2, 3, as well as for the mean solution v 0 . The RVE, shown in Fig. 10b, discretized with isoparametric quadratic triangular elements of average size h m = /10 using a three-point Gauss integration rule, is assigned to each macroscopic (1) and (2), i.e. n = 3. A similar preliminary analysis of this case has been reported in Rokoš et al. (2020a), which was limited by the capabilities of the employed (quasi-Newton) solver. Here a more detailed study is presented because a full Newton solver and a bifurcation analysis are used instead. Because an infinite specimen is considered, the DNS solution directly corresponds to the behaviour of a single periodic cell (i.e. RVE) subjected to a compressive load of biaxiality ratio γ. Since the deformation state is periodic, the ensemble average reduces to volume averaging, easily obtained from a single microstructural translation. The question arises, however, whether the micromorphic computational homogenization is capable of reproducing such a behaviour, i.e. yielding an affine mean field v 0 and constant micromorphic fields v i , while providing patterns of Eqs. (59)-(61) as the outcome of the analysis. Fig. 12 collects the results for a uniaxial compressive strain ε 22 = −0.05 (γ = ∞). As expected, the mean displacement in the e 1 direction is zero, whereas in the e 2 direction it is linear and corresponds to the applied nominal strain. The normalized micromorphic fields v i = v i / max k=1,2,3 v k (t, X) ∞ are also spatially constant with v 1 = 1 and v 2 = v 3 = 0, resulting in an activation of mode I and, consequently, pattern I (recall Eq. (59)). Moreover, the deformed RVE shape in Fig. 12a matches the DNS solution in Fig. 11a, corroborating further the validity of the micromorphic results.
The evolution of the magnitudes corresponding to the individual micromorphic fields as a function of ε 22 is shown for the overall applied deformation gradient and three values for γ in Fig. 13. Prior to bifurcation, all micromorphic fields remain zero. Upon reaching the critical strain, activation of the micromorphic fields starts exactly at the bifurcation point where a negative or sufficiently small lowest eigenvalue is observed and the system is perturbed towards the corresponding eigenvector. The correct patterns are triggered, i.e. v 1 = 0 while v 2 = v 3 = 0 for pattern I (γ = ∞), v 2 = v 3 = 0 while v 1 = 0 for pattern II (γ = 3 10 ), and v 1 = v 2 = v 3 = 0 for pattern III (γ = 1), recall Eqs. (59)-(61). To verify that the observed patterns in all three cases correspond to the correct solutions (i.e. the one related to the lowest strain energy), the existence of multiple local minima is explored. To this end, all micromorphic fields are initialized as constant fields, with magnitudes spanning the entire cube [ v 1 , v 2 , v 3 ] ∈ [0, 1] × [0, 1] × [0, 1] for a fixed applied overall strain which corresponds to a buckled state, while assuming the exact mean fields v 0 . It is found that although other patterns may yield stable local minima, the global minima always correspond to the correct combinations of modes. Note that similarly to the DNS, the multiplicities of the bifurcation points associated with the second and third pattern occur also for the micromorphic formulation. In that case, the associated buckling modes have a zero mean v 0 = 0 and spatially constant micromorphic fields, spanning the same vector space as in the DNS case. The only reliable procedure to identify the proper solution is then to explore each equilibrium path separately, opting for the one requiring the least amount of elastic strain energy. Although it may seem at this point that a bifurcation analysis is not of much benefit for a hexagonal stacking of holes, it reduces the number of possible combinations that would otherwise have to be considered as initial guesses for a quasi-Newton solver. In the case of pattern I, the benefit is clearly substantial. For pattern II the dimensionality reduces from three to two, whereas for pattern III the entire space of dimensionality three should be considered. From numerical evidence, however, mode combinations approximately matching the three patterns are typically observed as eigenmodes corresponding to the three lowest eigenvalues obtained during simulations, thus reducing all possible options (spanning a vector space of dimensionality three) to only three options.
The entire quadrant [ε 11 , ε 22 ] ∈ [0, −0.025] × [0, −0.025] is further explored with the micromorphic computational homogenization to provide a phase diagram of the hexagonal microstructure. The obtained result is plotted in Fig. 14a, where the normalized micromorphic fields v i are shown. Four regions A-C are distinguished, as depicted in the corresponding contour plot in Fig. 14b. In the first region, A, no pattern is triggered, because the critical bifurcation strain has not been exceeded yet. Pattern I occurs in region B, whereas the second pattern is triggered in region C. In the fourth region, D, a mixture of both patterns is observed, with the special configuration for γ = 1 corresponding to pattern III (denoted by the black solid line). A similar behaviour has been observed for hexagonal honeycomb  (c), the corresponding contour plot is shown in (b) (normalization constant max k v k (t, X) ∞ = 3.3), whereas a circumferential section for ε ii = ε ii /0.025 in (c) (normalization constant max k v k (t, X) ∞ = 2.8). The circumferential section is taken along the black dashed curves in (a) and (b) in a clockwise direction.
structures in the work of Okumura et al. (2002); see Fig. 10 therein. Fig. 14c plots a circumferential section through the phase diagram of 14a taken along the dashed black curve highlighted in Fig. 14b, including its extension to the other three strain quadrants. The magnitudes of the individual normalized micromorphic fields v i are plotted as a function of angle ϑ ∈ [0, 2π], spanning the entire circle. The angle starts from the (ε 22 = −1, ε 11 = 0) direction and sweeps clockwise. Again, Fig. 14c confirms that equal magnitudes of all modes occur for ϑ = π/4. Furthermore, it is clearly visible which strain combinations yield which microstructural pattern. For instance, close to ϑ = 5π/8 and 15π/8 we notice that even though one of the applied strains ε ii is positive, a pattern transformation occurs due to a large negative magnitude of the other compressive strain.

Summary and Conclusions
This contribution has extended a recently developed micromorphic computational homogenization framework for mechanical metamaterials with a full Newton solver. The micromorphic framework decomposes the kinematic field by exploiting prior knowledge on the typical patterning modes, allowing to accurately capture non-local effects present in the microstructure. The derivation and implementation of a full Newton solver for this framework has been provided, including analytical expressions for the first and second variations of the total effective potential energy. Significant gains have been obtained compared to the existing quasi-Newton implementation, in particular with respect to the bifurcation analysis, which is essential for applications of elastomeric mechanical metamaterials relying on local and global buckling. Two examples have been tested to demonstrate the capabilities of the presented numerical scheme. In the first example, a metamaterial column consisting of a microstructure with a square stacking of holes has been analysed for various slenderness ratios, for which a competition between the local and global buckling exists. The second example elaborated a uniformly loaded infinite specimen with a hexagonal stacking of holes, which may buckle into different patterns depending on the loading direction.
The main conclusions of this paper can be summarized as follows: 1. The developed full Newton solver for the micromorphic computational homogenization framework is robust and efficient. 2. The micromorphic approach captures the behaviour of the reference Direct Numerical Simulation (DNS) accurately in terms of both local and global buckling as well as the pattern magnitudes. 3. The nominal stresses are reproduced by the micromorphic framework with a good accuracy, although the post-bifurcation results are in general systematically overestimated compared to DNS. The maximum error in terms of the critical buckling stress corresponding to the first instability point does not exceed 12%, and decreases with increasing scale ratio down to approximately 7%. 4. The buckling strain is captured with a higher accuracy compared to the nominal buckling stress. The relative error stays below 3% for local and 6% for global buckling, and decreases down to 0.5% for local and 3% for global buckling with increasing scale ratio. 5. The micromorphic scheme reproduces DNS results correctly even in the case of a hexagonal stacking of holes, for which multiple patterning modes occur. It predicts the correct patterns for the loading directions considered.
The full Newton solver presented here greatly reduces the dependency of the solution on the initial guess by perturbing the system towards the correct direction when a bifurcation point is encountered; therefore, it provides an indispensable numerical tool for modelling instability-based mechanical metamaterials.