A Newton Solver for Micromorphic Computational Homogenization Enabling Multiscale Buckling Analysis of Pattern-Transforming Metamaterials
S.E.H.M. van Bree, O. Rokoš, R.H.J. Peerlings, M. Doškář, M.G.D. Geers
AA Newton Solver for Micromorphic Computational HomogenizationEnabling Multiscale Buckling Analysis of Pattern-TransformingMetamaterials (cid:73)
S.E.H.M. van Bree a , O. Rokoˇs a,b, ∗ , R.H.J. Peerlings a , M. Doˇsk´aˇr b , M.G.D. Geers a a Mechanics of Materials, Department of Mechanical Engineering, Eindhoven University of Technology,P.O. Box 513, 5600 MB Eindhoven, The Netherlands b Department of Mechanics, Faculty of Civil Engineering, Czech Technical University in Prague,Th´akurova 7, 166 29 Prague 6, Czech Republic.
Abstract
Mechanical metamaterials feature engineered microstructures designed to exhibit exotic,and often counter-intuitive, effective behaviour such as negative Poisson’s ratio or negativecompressibility. Such a specific response is often achieved through instability-induced trans-formations of the underlying periodic microstructure into one or multiple patterning modes.Due to a strong kinematic coupling of individual repeating microstructural cells, non-localbehaviour and size effects emerge, which cannot easily be captured by classical homogeniza-tion schemes. In addition, the individual patterning modes can mutually interact in spaceas well as in time, while at the engineering scale the entire structure can buckle globally. Forefficient numerical predictions of macroscale engineering applications, a micromorphic com-putational homogenization scheme has recently been developed (Rokoˇs et al.,
J. Mech. Phys.Solids , 119–137, 2019). Although this framework is in principle capable of accounting forspatial and temporal interactions between individual patterning modes, its implementationrelied on a gradient-based quasi-Newton solution technique. This solver is suboptimal be-cause (i) it has sub-quadratic convergence, and (ii) the absence of Hessians does not allow forproper bifurcation analyses. Given that mechanical metamaterials often rely on controlledinstabilities, these limitations are serious. Addressing them will reduce the dependency ofthe solution on the initial guess by perturbing the system towards the correct deformationwhen a bifurcation point is encountered. Eventually, this enables more accurate and reliablemodelling and design of metamaterials. To achieve this goal, a full Newton method, entailingall derivations and definitions of the tangent operators, is provided in detail in this paper.The construction of the macroscopic tangent operator is not straightforward due to specific (cid:73)
The post-print version of this article is published in
Comput. Methods Appl. Mech. Engrg. ,10.1016/j.cma.2020.113333 ∗ Corresponding author.
Email addresses: [email protected] (S.E.H.M. van Bree),
[email protected] (O. Rokoˇs),
[email protected] (R.H.J. Peerlings), [email protected] (M. Doˇsk´aˇr),
[email protected] (M.G.D. Geers)
Preprint submitted to Comput. Methods Appl. Mech. Engrg. September 4, 2020 a r X i v : . [ c s . C E ] S e p odel assumptions on the decomposition of the underlying displacement field pertinent tothe micromorphic framework, involving orthogonality constraints. Analytical expressionsfor the first and second variation of the total potential energy are given, and the completealgorithm is listed. The developed methodology is demonstrated with two examples in whicha competition between local and global buckling exists and where multiple patterning modesemerge. The numerical results indicate that local to global buckling transition can be pre-dicted within a relative error of 6% in terms of the applied strains. The expected patterncombinations are triggered even for the case of multiple patterns. Keywords:
Mechanical metamaterials, computational homogenization, micromorphiccontinuum, Newton method, bifurcation analysis
1. Introduction
Acting like a carefully engineered structure, rather than a standard bulk material, is acommon characteristic of mechanical metamaterials. Recent advances in 3D printing andadditive manufacturing enable the production of such structures on a relatively small scale,allowing to treat them as a homogeneous medium. Metamaterials are typically designed toexhibit an exotic behaviour which cannot be found in nature, such as a negative compress-ibility (Nicolaou and Motter, 2012), negative Poisson’s ratio (Kolken and Zadpoor, 2017), ora high stiffness with an ultra low density (Zheng et al., 2014). In this contribution, we focuson elastomeric mechanical metamaterials, which under compression exhibit microstructuralbuckling resulting in a pattern transformation. Such a transformation induces an abruptchange in effective properties including Young’s modulus and Poisson’s ratio, with envi-sioned applications in, e.g., soft robotics (see Yang et al., 2015; Mark et al., 2016; Mirzaaliet al., 2018).Because elastomeric mechanical metamaterials rely mostly on local instabilities in theirmicrostructural morphology, large deformations, rotations, and strains occur. In particular,the microstructure undergoes a pattern transformation, due to coordinated buckling of theunderlying microstructure, resulting in a strongly non-local behaviour. If the specimen isrestricted, e.g. by applied essential boundary conditions, the expected pattern cannot fullydevelop, 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 thebulk of the (meta)material, which may significantly influence the overall response even atthe engineering scale. Such a configuration is depicted in Fig. 1a, in which an example ofan elastomeric metamaterial beam subjected to compressive load is shown. The buckledpattern vanishes close to the two vertical boundaries, resulting in a stiffening effect. Theextent to which the boundary layers influence the effective mechanical behaviour dependson the ratio of their thickness and the overall size of the specimen, or more generally on thescale ratio defined as the ratio between the overall size of the specimen H relative to thetypical size of the microstructural features (cid:96) , i.e. H/(cid:96) .For predictive modelling of engineering-scale applications, it is important to accuratelyyet efficiently predict the overall mechanical response of the structure. To this end, homog-2 a) local buckling (b) global buckling (cid:107) (cid:126)u (cid:107) /u D Figure 1: An elastomeric mechanical metamaterial containing circular holes in a square packing exhibitinga pattern transformation under a compressive load induced by clamping the right-hand side of the specimenwhile prescribing a horizontal displacement u D at the left-hand side. Restricted evolution of the patterningnear 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) islocalized only in the compressive regions of the specimen. The colour indicates the pointwise magnitude ofthe displacement field (cid:126)u relative to u D . enization techniques are employed, replacing the complex microstructural behaviour withan equivalent continuum model. However, due to the non-locality, patterning, and buck-ling of the microstructure, homogenization of mechanical metamaterials presents a difficultchallenge. First-order computational homogenization, outlined e.g. by Kouznetsova et al.(2001), can significantly reduce computing time compared to full scale Direct NumericalSimulations (DNS). However, its inherent assumptions on locality and scale separation pre-vent it from predicting any size effects in the microstructure. Ameen et al. (2018) showedthat as a consequence relative errors up to 40% can be induced in terms of force quantities bythe first-order method in the post-bifurcation regime for small scale ratios. With an increas-ing scale ratio, the accuracy of the predicted overall behaviour typically improves, with anexact match for H/(cid:96) → ∞ . The second-order computational homogenization (Kouznetsovaet al., 2004), which incorporates the gradient of the macroscopic deformation gradient inthe micro-to-macro transition, permits to reflect non-locality and ensuing size-effects. Theeffective behaviour captured by this method coincides satisfactorily with DNS results evenfor low scale ratios, albeit at the cost of additional complexity stemming from a higher-ordercontinuum formulation at the macroscale level; see Sperling et al. (2020) for more details.Recently, Rokoˇs et al. (2019) proposed a micromorphic computational homogenization frame-work, specifically designed to predict the effective behaviour of mechanical metamaterialsby decomposing the displacement field into three components: (i) a smooth, mean displace-ment field, (ii) a spatially correlated microfluctuation field, and (iii) an uncorrelated, localmicrofluctuation field. This decomposition ensures an adequate performance and accuracyby introducing prior knowledge on the patterning fluctuation. See Sperling et al. (2020)and Rokoˇs et al. (2020a) for more details.A major practical limitation of the previously reported implementations of the micro-morphic computational homogenization framework is that a quasi-Newton solution methodhas been employed. This method is not as efficient as a full Newton scheme, it has shownto be quite sensitive to the initialization of the (equilibrium) iteration process, and in par-ticular to the perturbations applied to trigger buckling and, related to the latter point, it3oes 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 atransition in the effective mechanical properties, see Fig. 1a) or global buckling (indicatinga 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 microscopicbifurcation has been studied by Saiki et al. (2002), whereas Wadee and Farsi (2015) haveinvestigated geometrical effects on the buckling behaviour of cellular structures in whicha transition between local and global buckling is observed as a function of the specimenslenderness. Experimentally and numerically, Niknam and Akbarzadeh (2018) have com-pared in-plane and out-of-plane buckling of various types and sizes of architected cellularstructures. Specimens with a hexagonal honeycomb microstructure, which exhibit multiplebuckling patterns under different compressive biaxiality ratios, were studied by Ohno et al.(2002a). Rokoˇs et al. (2020a) demonstrated the same multi-pattern character for hexago-nally stacked cells with circular holes using the micromorphic homogenization framework.Obtaining correct patterning of the microstructure required, nevertheless, intervention of theuser based on insight in the mechanics of the system, precisely because no reliable solver wasavailable to tackle the buckling.The main goal of this paper is to derive the tangent operator for the micromorphiccomputational homogenization framework, enabling a more efficient and robust solutionprocedure using a full Newton algorithm and allowing for bifurcation analyses (Miehe andBayreuther, 2007). Unlike the derivation of the Hessians (i.e., the macroscopic tangents orstiffnesses) for first-order computational homogenization—see detailed explanation in (Mieheand Koch, 2002) or (Miehe, 2003)—, the Hessians for the micromorphic scheme require anon-trivial extension. Additional orthogonality constraints acting within each RepresentativeVolume Element (RVE) need to be enforced in order to guarantee uniqueness of the adoptedkinematic decomposition. Using variational calculus, the first variation of the averagedenergy resulting in the microscopic and macroscopic governing equations is derived, followedby the second variation from which the micro-, and coupling macro-Hessians can be obtained.Following Rokoˇs et al. (2020a), the formulation involves an arbitrary number of patterningmodes, and introduces a slight reformulation of the orthogonality constraints with respectto gradients of individual modes as compared to the original framework (Rokoˇs 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 andstiffnesses are constructed, from which the local microfluctuation fields are condensed out,yielding a macroscopic Newton algorithm. We illustrate the performance of the method withtwo examples, one focusing on local versus global buckling of a metamaterial column, andone on local patterning of a hexagonally-voided microstructure.The remainder of this paper is organized as follows. After recalling the kinematic de-composition pertinent to the micromorphic framework, Section 2 details the derivation ofthe first and second variation of the ensemble averaged energy, resulting in both macro-and microscopic governing equations accompanied by the relevant macro- and microscopic4essians. Employing standard FE procedures, Section 3 describes the discretization of thegoverning equations at both scales and addresses the bifurcation analysis. Section 4 illus-trates the developed methodology with two examples: (i) a compressed metamaterial columnwith a varying slenderness ratio in which a competition between local and global bucklingexists, and (ii) a microstructure with hexagonally-stacked holes subjected to biaxial com-pressive loading exhibiting three distinct pattern transformations. Finally, the summary andconclusions follow in Section 5.Throughout the paper, the following notational conventions are used - scalars a , - vectors (cid:126)a , - position vector in the reference config-uration (cid:126)X = X (cid:126)e + X (cid:126)e , - second-order tensors A = A ij (cid:126)e i (cid:126)e j , - third-order tensors A = A ijk (cid:126)e i (cid:126)e j (cid:126)e k , - fourth-order tensors A = A ijkl (cid:126)e i (cid:126)e j (cid:126)e k (cid:126)e l , - matrices A and column matrices a , - (cid:126)a · (cid:126)b = a i b i , - A · (cid:126)b = A ij b j (cid:126)e i , - A · B = A ik B kj (cid:126)e i (cid:126)e j , - A : B = A ij B ji , - transpose A T , A T ij = A ji , - right transpose A RT , A RT ijkl = A ijlk , - left transpose A LT , A LT ijkl = A jikl , - gradient operator (cid:126) ∇ (cid:126)a = ∂a j ∂X i (cid:126)e i (cid:126)e j , - divergence operator (cid:126) ∇ · (cid:126)a = ∂a i ∂X i , - integration (cid:104) f ( (cid:126)X m ) (cid:105) Ω • = (cid:82) Ω • f ( (cid:126)X m ) d (cid:126)X m , - derivatives of scalar functions with re-spect to second-order tensors δ Ψ( F ; δ F ) = dd h Ψ( F + hδ F ) (cid:12)(cid:12)(cid:12)(cid:12) h =0 = ∂ Ψ( F ) ∂ F : δ F ,where Einstein’s summation convention is adopted on repeated indices i , j , k , l , and (cid:126)e i , i = 1 , , denote the basis vectors of a two-dimensional Cartesian coordinate frame.
2. Reformulation of the Micromorphic Computational Homogenization Frame-work and Derivation of the Tangents
The micromorphic computational homogenization framework (Rokoˇs et al., 2019), de-picted schematically in Fig. 2, relies on the decomposition of the kinematic field (cid:126)u intothe mean effective displacement (cid:126)v , long range correlated fluctuation components v i (cid:126)ϕ i , i = 1 , . . . , n , and the remaining local microfluctuation field (cid:126)w , i.e. (cid:126)u ( (cid:126)X ) = (cid:126)v ( (cid:126)X ) + n (cid:88) i =1 v i ( (cid:126)X ) (cid:126)ϕ i ( (cid:126)X ) + (cid:126)w ( (cid:126)X ) . (1)The vector field (cid:126)ϕ i corresponds to the i -th patterning mode of the underlying microstruc-ture, whereas the scalar field v i regulates, spatially and in time, its magnitude. Unlike theoriginal formulation (Rokoˇs et al., 2019), which considered only one such mode, we consider5ere an arbitrary number of modes, n ; see also Rokoˇs et al. (2020a). Because in general itmay not be possible to control the positioning of the microstructure relative to the speci-men’s boundary, all possible microstructural translations should be taken into account viaensemble averaging (cf. Ameen et al., 2018). The micromorphic scheme avoids this costlyprocedure by approximating the mechanical state of a point in a translated microstructureby evaluating the mechanical state of a microstructurally equivalent point in the referencemicrostructure. In addition, a separation of scales into a macroscopic position vector (cid:126)X anda microscopic position vector (cid:126)X m is introduced, assuming the fields (cid:126)v and v i to vary slowlyover a close vicinity of each macroscopic point (cid:126)X spanned by a microscopic RepresentativeVolume Element (RVE) with a domain Ω m . Consequently, the microfluctuation field (cid:126)w iscomputed only locally over each RVE, and is independent for RVEs associated with distinctmacroscopic points (cid:126)X ; they communicate only by means of the macroscopic fields (cid:126)v and v i .Therefore, (cid:126)v and v i become functions of the macroscopic position only, whereas (cid:126)ϕ i and (cid:126)w are functions of the microscopic as well as the macroscopic position. Because the patterningmode (cid:126)ϕ i is the same for each macroscopic point, it eventually depends on the microscopicposition vector (cid:126)X m only. Using a first-order Taylor expansion and the above considerations,the decomposition of Eq. (1) can be approximated as (cid:126)u ( (cid:126)X, (cid:126)X m ) ≈ (cid:126)v ( (cid:126)X ) + (cid:126)X m · (cid:126) ∇ (cid:126)v ( (cid:126)X )+ n (cid:88) i =1 [ v i ( (cid:126)X ) + (cid:126)X m · (cid:126) ∇ v i ( (cid:126)X )] (cid:126)ϕ i ( (cid:126)X m )+ (cid:126)w ( (cid:126)X, (cid:126)X m ) , (cid:126)X ∈ Ω , (cid:126)X m ∈ Ω m . (2)For more details on the decomposition (2), the reader is referred to Rokoˇs et al. (2019)—albeit for a single mode (cid:126)ϕ . The fields (cid:126)v , v i , and (cid:126)w , are unknown and need to be solved for,while the individual patterning modes (cid:126)ϕ i are characteristic for the underlying microstructuralmorphology and are assumed to be known a priori, computed either from a Bloch-type analy-sis (Bertoldi et al., 2008), estimated analytically from full-scale numerical simulations (Rokoˇset al., 2019), or identified experimentally (Maraghechi et al., 2020). Because the patterningmodes (cid:126)ϕ i are defined with respect to the reference microscopic configuration (cid:126)X m , the effectof the macroscopic rotations needs to be factored out from the macroscopic deformationgradients, F M = I + ( (cid:126) ∇ (cid:126)v ) T , using the polar decomposition in Eq. (2), i.e. F M = R M · U M ,where R M is the macroscopic rotation tensor and U M is the macroscopic stretch tensor.This effectively means that the term (cid:126) ∇ (cid:126)v 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, how-ever, omitted hereafter to simplify all derivations, and can even be neglected in the limitof small rotations as is the case for both examples shown in Section 4. Note also that theadditional micromorphic fields v i in Eqs. (1) and (2) relate directly to the displacementfield (cid:126)u . Such a setup contrasts with standard micromorphic formulations, in which the mi-crodeformation χ typically relates to the gradient of (cid:126)u , see e.g. (Forest and Trinh, 2011,Eqs. (29)–(32)). The decomposition of Eq. (2) further suggests that when all micromorphic6 (cid:126)e (cid:126)e (cid:126)X Macroscale boundary value problem, Eq. (14) (cid:126) ∇ (cid:126)v v i , (cid:126) ∇ v i (cid:126)e (cid:126)e (cid:126)X m Ω m Θ , Π i , (cid:126) Λ i Stiffnesses Solve for (cid:126)w
Microscale boundary value problem, Eq. (13)Figure 2: A schematic representation of the micromorphic computational homogenization framework. Themacroscopic displacement gradient (cid:126) ∇ (cid:126)v , micromorphic fields v i , and their spatial gradients (cid:126) ∇ v i at a macro-scopic point (cid:126)X are sampled and sent to the microscale, where a boundary value problem for the microfluc-tuation field (cid:126)w is solved. Based on the solution (cid:126)w , the macroscopic stresses, Θ , Π i , (cid:126) Λ i , and stiffnesses arecomputed and passed back to the macroscale, where the macroscopic boundary value problem is assembledand solved. fields v i vanish, the ansatz of the standard first-order computational homogenization is re-covered, i.e. (cid:126)u ( (cid:126)X, (cid:126)X m ) = (cid:126)v ( (cid:126)X ) + (cid:126)X m · (cid:126) ∇ (cid:126)v ( (cid:126)X ) + (cid:126)w ( (cid:126)X, (cid:126)X m ), cf. e.g. (Geers et al., 2010,Eq. (1)).The uniqueness of the decomposition (2) within an RVE is guaranteed by introducingthe following additional orthogonality conditions: (cid:10) (cid:126)w ( (cid:126)X, (cid:126)X m ) (cid:11) Ω m = (cid:126) , (3) (cid:10) (cid:126)w ( (cid:126)X, (cid:126)X m ) · (cid:126)ϕ i ( (cid:126)X m ) (cid:11) Ω m = 0 , i = 1 , . . . , n, (4) (cid:10) (cid:126)w ( (cid:126)X, (cid:126)X m ) · [ (cid:126)ϕ i ( (cid:126)X m ) (cid:126)X m ] (cid:11) Ω m = (cid:126) , i = 1 , . . . , n. (5)Recall that throughout this contribution the angle brackets indicate the integration overthe domain specified by the subscript. The first condition (3) requires zero mean of (cid:126)w over Ω m , effectively eliminating rigid body translations. The second condition establishesthe uniqueness in terms of the patterning modes (cid:126)ϕ i themselves, because, upon assuming ahomogeneous state with (cid:126) ∇ v i = (cid:126)
0, for instance, the patterning can be equally well representedby the product of the micromorphic field and the patterning mode, i.e. v i (cid:126)ϕ i , or by themicrofluctuation field (cid:126)w while keeping v i = 0. The last orthogonality condition has not beenpreviously introduced in the original formulation of Rokoˇs et al. (2019). It follows, however,from the decomposition outlined in Eg. (2), and eliminates non-uniqueness issues related tolinearly varying magnitudes of the patterning fields (cid:126)ϕ i , which may cause spurious macroscopic7 Ω Rm ∂ Ω Tm ∂ Ω Bm ∂ Ω Lm Ω m (cid:96)d(cid:126)e (cid:126)e (cid:126)X m Figure 3: A schematic example of a pattern-transforming 2 × ∂ Ω +m = ∂ Ω Tm ∪ ∂ Ω Rm and ∂ Ω − m = ∂ Ω Bm ∪ ∂ Ω Lm for theimplementation of periodic boundary conditions. oscillations for microstructures with a hexagonal stacking of holes. This condition actsmainly as a stabilization condition and does not significantly affect the overall mechanicalresponse. Following the same reasoning as in first-order computational homogenization, theremaining orthogonality condition, which might arise from the non-uniqueness related to theansatz (2) (i.e. orthogonality with respect to (cid:126)X m · (cid:126) ∇ (cid:126)v ( (cid:126)X ), or (cid:10) (cid:126)w ( (cid:126)X, (cid:126)X m ) (cid:126)X m (cid:11) Ω m = foran arbitrary (cid:126) ∇ (cid:126)v ), is accounted for differently. As a well-accepted modelling choice, whichprovides accurate results in the first- as well as second-order computational homogenizationschemes, a periodicity constraint on (cid:126)w is adopted ensuring this orthogonality, i.e. (cid:74) (cid:126)w ( (cid:126)X, (cid:126)X m ) (cid:75) = (cid:126) , (cid:126)X m ∈ ∂ Ω +m , (6)where (cid:74) (cid:126)w ( (cid:126)X, (cid:126)X m ) (cid:75) = (cid:126)w ( (cid:126)X, ∂ Ω +m ) − (cid:126)w ( (cid:126)X, ∂ Ω − m ) denotes the jump of the field (cid:126)w ( (cid:126)X, (cid:126)X m ) onthe RVE boundary split into two parts: ∂ Ω m = ∂ Ω +m ∪ ∂ Ω − m . As an example, a 2 × ∂ Ω +m = ∂ Ω Tm ∪ ∂ Ω Rm and ∂ Ω − m = ∂ Ω Bm ∪ ∂ Ω Lm . The same approach is used for polygons with multiple edges—see e.g. ahead to Fig. 10b for the case of a hexagonal RVE. The periodicity constraint incombination with the condition of Eq. (3) also eliminates all rigid body modes of (cid:126)w . Al-though periodic boundary conditions for (cid:126)w (and thus antiperiodic RVE boundary tractions)are adopted, such a choice might not necessarily result in the optimal performance of theproposed micromorphic scheme. See e.g. Forest and Trinh (2011), where antiperiodic mi-crofluctuation fields (cid:126)w with periodic RVE boundary tractions have been observed. Yet, itwill be demonstrated in the results Section 4 that the periodicity constraints (6) provide anadequate accuracy here. Because we are restricting ourselves to hyperelastic materials, the unknown parts (cid:126)v ( (cid:126)X ), v i ( (cid:126)X ), and (cid:126)w ( (cid:126)X, (cid:126)X m ) of the solution (cid:126)u ( (cid:126)X, (cid:126)X m ) according to (2) can be found by minimizingthe total potential energy of the system while accounting for the constraints introduced8bove, i.e.( (cid:126)v ( (cid:126)X ) , v ( (cid:126)X ) , . . . , v n ( (cid:126)X ) , (cid:126)w ( (cid:126)X, (cid:126)X m )) ∈ arg min (cid:126)u ( (cid:126)X, (cid:126)X m ) max (cid:126)µ ( (cid:126)X ) , ν ( (cid:126)X ) η ( (cid:126)X ) , (cid:126)λ ( (cid:126)X, (cid:126)X m ) L ( (cid:126)u, (cid:126)µ, ν, η, (cid:126)λ ) , (7)where the displacement field (cid:126)u ( (cid:126)X, (cid:126)X m ) is, through Eq. (2), also a function of all unknownmacroscopic, (cid:126)v ( (cid:126)X ), v ( (cid:126)X ), . . . , v n ( (cid:126)X ), and microscopic, (cid:126)w ( (cid:126)X, (cid:126)X m ), quantities in addition tothe two spatial variables (cid:126)X and (cid:126)X m , and where the fields (cid:126)µ ( (cid:126)X ), ν ( (cid:126)X ) = [ ν ( (cid:126)X ) , . . . , ν n ( (cid:126)X )],and η ( (cid:126)X ) = [ (cid:126)η ( (cid:126)X ) , . . . , (cid:126)η n ( (cid:126)X )] collect the Lagrange multipliers pertinent to rigid bodymodes and orthogonality constraints with respect to individual patterning modes as theircomponents. While the bulk constraints of Eqs. (3)–(5) are considered for each RVE sep-arately, i.e. (cid:126)µ ( (cid:126)X ), ν ( (cid:126)X ), and η ( (cid:126)X ) are a function of the macroscopic position vector (cid:126)X only, the periodicity constraint (6) is a continuous function on the boundary of each RVEand therefore, it is also a function of the microscopic position vector (cid:126)X m , i.e. (cid:126)λ ( (cid:126)X, (cid:126)X m ). InEq. (7), the Lagrangian L = E + C consists of the total potential energy EE ( (cid:126)u ( (cid:126)X, (cid:126)X m )) = 1 | Ω m | (cid:68)(cid:10) Ψ( (cid:126)X m , F ) (cid:11) Ω m (cid:69) Ω (8)and the constraint term CC ( (cid:126)u ( (cid:126)X, (cid:126)X m ) , (cid:126)µ ( (cid:126)X ) , ν ( (cid:126)X ) , η ( (cid:126)X ) , (cid:126)λ ( (cid:126)X, (cid:126)X m )) = 1 | Ω m | (cid:68) (cid:126)µ ( (cid:126)X ) · (cid:10) (cid:126)w ( (cid:126)X, (cid:126)X m ) (cid:11) Ω m + n (cid:88) i =1 ν i ( (cid:126)X ) (cid:10) (cid:126)w ( (cid:126)X, (cid:126)X m ) · (cid:126)ϕ i ( (cid:126)X m ) (cid:11) Ω m + n (cid:88) i =1 (cid:126)η i ( (cid:126)X ) · (cid:10) (cid:126)w ( (cid:126)X, (cid:126)X m ) · [ (cid:126)ϕ i ( (cid:126)X m ) (cid:126)X m ] (cid:11) Ω m − (cid:10) (cid:126)λ ( (cid:126)X, (cid:126)X m ) · (cid:74) (cid:126)w ( (cid:126)X, (cid:126)X m ) (cid:75) (cid:11) ∂ Ω +m (cid:69) Ω , (9)where the minus sign in front of the last constraint term is adopted for consistency rea-sons (Miehe and Bayreuther, 2007, Section 2.2.2), and the remaining terms take a plussign 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ˇs et al. (2019), where theconstraint term C was incomplete. A hyperelastic constitutive law specified through theenergy density function Ψ( (cid:126)X m , F ) is adopted for the description of the material behaviour; F ( (cid:126)u ( (cid:126)X, (cid:126)X m )) = I + ( (cid:126) ∇ m (cid:126)u ( (cid:126)X, (cid:126)X m )) T is the deformation gradient and (cid:126) ∇ m = (cid:126)e i ∂/∂X m ,i is themicroscopic gradient operator defined in the reference configuration. Note that the explicitdependency of F on (cid:126)u has been dropped in Eq. (8), and will be omitted for brevity hereafteras well. A minimizer of the total potential energy can be found by taking the Gˆateaux derivative(i.e. the first variation) of the Lagrangian L and requiring it to vanish:0 = δ L ( (cid:126)u, (cid:126)µ, ν, η, (cid:126)λ ; δ(cid:126)u, δ(cid:126)µ, δν, δη, δ(cid:126)λ ) = 1 | Ω m | (cid:68)(cid:10) P ( (cid:126)X m , F ) : (cid:126) ∇ m δ(cid:126)u ( (cid:126)X, (cid:126)X m ) (cid:11) Ω m (cid:126)µ ( (cid:126)X ) · (cid:10) δ (cid:126)w ( (cid:126)X, (cid:126)X m ) (cid:11) Ω m + δ(cid:126)µ ( (cid:126)X ) · (cid:10) (cid:126)w ( (cid:126)X, (cid:126)X m ) (cid:11) Ω m + n (cid:88) i =1 ν i ( (cid:126)X ) (cid:10) δ (cid:126)w ( (cid:126)X, (cid:126)X m ) · (cid:126)ϕ i ( (cid:126)X m ) (cid:11) Ω m + n (cid:88) i =1 δν i ( (cid:126)X ) (cid:10) (cid:126)w ( (cid:126)X, (cid:126)X m ) · (cid:126)ϕ i ( (cid:126)X m ) (cid:11) Ω m (10)+ n (cid:88) i =1 (cid:126)η i ( (cid:126)X ) · (cid:10) δ (cid:126)w ( (cid:126)X, (cid:126)X m ) · [ (cid:126)ϕ i ( (cid:126)X m ) (cid:126)X m ] (cid:11) Ω m + n (cid:88) i =1 δ(cid:126)η i ( (cid:126)X ) · (cid:10) (cid:126)w ( (cid:126)X, (cid:126)X m ) · [ (cid:126)ϕ i ( (cid:126)X m ) (cid:126)X m ] (cid:11) Ω m − (cid:10) (cid:126)λ ( (cid:126)X, (cid:126)X m ) · (cid:74) δ (cid:126)w ( (cid:126)X, (cid:126)X m ) (cid:75) + δ(cid:126)λ ( (cid:126)X, (cid:126)X m ) · (cid:74) (cid:126)w ( (cid:126)X, (cid:126)X m ) (cid:75) (cid:11) ∂ Ω +m (cid:69) Ω , where the local first Piola–Kirchhoff stress tensor P ( (cid:126)X m , F ) is defined as P ( (cid:126)X m , F ( (cid:126)u ( (cid:126)X, (cid:126)X m ))) = ∂ Ψ( (cid:126)X m , F ( (cid:126)u ( (cid:126)X, (cid:126)X m ))) ∂ F T . (11)Making use of the decomposition introduced in Eq. (2), (cid:126) ∇ m δ(cid:126)u reads (cid:126) ∇ m δ(cid:126)u ( (cid:126)X, (cid:126)X m ) = (cid:126) ∇ δ(cid:126)v ( (cid:126)X ) + n (cid:88) i =1 (cid:126) ∇ δv i ( (cid:126)X ) (cid:126)ϕ i ( (cid:126)X m )+ n (cid:88) i =1 (cid:2) δv i ( (cid:126)X ) + (cid:126)X m · (cid:126) ∇ δv i ( (cid:126)X ) (cid:3) (cid:126) ∇ m (cid:126)ϕ i ( (cid:126)X m ) + (cid:126) ∇ m δ (cid:126)w ( (cid:126)X, (cid:126)X m ) . (12)Upon substituting this expansion, application of the divergence theorem, and rearrangementof individual terms, Eq. (10) leads to a set of microscopic and macroscopic balance equations.At the microscale, only the variations of the microfluctuation field (cid:126)w ( (cid:126)X, (cid:126)X m ) and its relatedLagrange multipliers (cid:126)µ ( (cid:126)X ), ν i ( (cid:126)X ), (cid:126)η i ( (cid:126)X ), and (cid:126)λ ( (cid:126)X, (cid:126)X m ) matter. Since Eq. (10) must holdfor arbitrary variations of these quantities, the following set of microscale balance equationsresults δ (cid:126)w : (cid:126) ∇ m · P T = (cid:126)µ + n (cid:88) i =1 ν i (cid:126)ϕ i + n (cid:88) i =1 (cid:126)η i · ( (cid:126)ϕ i (cid:126)X m ) , in Ω m , P · (cid:126)N m = ± (cid:126)λ, on ∂ Ω ± m ,δ(cid:126)λ : periodicity constraint for (cid:126)w , Eq. (6) ,δ(cid:126)µ, δν, δη : kinematic constraints for (cid:126)w , Eqs. (3)–(5) . (13)As a consequence of the constraint term introduced in Eq. (9), the right hand side of thefirst 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 orthogonalityconstraints. These terms were not present in the original formulation, where the constraintequations were enforced differently. Although the body forces (cid:126)µ , v i (cid:126)ϕ i , and (cid:126)η i · ( (cid:126)ϕ i (cid:126)X m ), aredirectly associated with the orthogonality constrains of Eqs. (3)–(5) in the form of Lagrangemultipliers, they can be introduced also directly at the level of governing equations, seee.g. (Yvonnet et al., 2020, Eqs. (8) and (9)). At the macroscale, only the slowly varying10elds (cid:126)v ( (cid:126)X ) and v i ( (cid:126)X ) are relevant, and their governing equations read δ(cid:126)v : (cid:40) (cid:126) ∇ · Θ T = (cid:126) , in Ω , Θ · (cid:126)N = (cid:126) , on Γ N , δv i : (cid:40) (cid:126) ∇ · (cid:126) Λ i − Π i = 0 , in Ω ,(cid:126) Λ i · (cid:126)N = 0 , on Γ N , i = 1 , . . . , n, (14)with the following definitions of macroscopic stress-like quantities: Θ ( (cid:126)X ) = 1 | Ω m | (cid:10) P ( (cid:126)X m , F ) (cid:11) Ω m , (15)Π i ( (cid:126)X ) = 1 | Ω m | (cid:10) P ( (cid:126)X m , F ) : (cid:126) ∇ m (cid:126)ϕ i ( (cid:126)X m ) (cid:11) Ω m , (16) (cid:126) Λ i ( (cid:126)X ) = 1 | Ω m | (cid:10) P T ( (cid:126)X m , F ) · (cid:126)ϕ i ( (cid:126)X m ) + (cid:126)X m [ P ( (cid:126)X m , F ) : (cid:126) ∇ m (cid:126)ϕ i ( (cid:126)X m )] (cid:11) Ω m . (17)The dependence of the homogenized stress quantities on the macroscopic position (cid:126)X orig-inates from the dependence of the deformation gradient F ( (cid:126)u ( (cid:126)X, (cid:126)X m )) on the macroscopicposition through all smooth fields (cid:126)v ( (cid:126)X ) and v i ( (cid:126)X ), cf. also Eq. (2). The second variation of the Lagrangian L reads δ L ( (cid:126)u, (cid:126)µ, ν, η, (cid:126)λ ; δ(cid:126)u, δ(cid:126)µ, δν, δη, δ(cid:126)λ ) = 1 | Ω m | (cid:68)(cid:10) (cid:126) ∇ m δ(cid:126)u : D : (cid:126) ∇ m δ(cid:126)u (cid:11) Ω m + 2 δ(cid:126)µ · (cid:10) δ (cid:126)w (cid:11) Ω m + 2 n (cid:88) i =1 δν i (cid:10) δ (cid:126)w · (cid:126)ϕ i (cid:11) Ω m + 2 n (cid:88) i =1 δ(cid:126)η i · (cid:10) δ (cid:126)w · [ (cid:126)ϕ i (cid:126)X m ] (cid:11) Ω m − (cid:10) δ(cid:126)λ · (cid:74) δ (cid:126)w (cid:75) (cid:11) ∂ Ω +m (cid:69) Ω , (18)where the spatial dependencies on (cid:126)X and (cid:126)X m have been dropped for brevity, and where thelocal stiffness tensor D m ( (cid:126)X m , F ) has been introduced as D ( (cid:126)X m , F ) = ∂ Ψ( (cid:126)X m , F ( (cid:126)u ( (cid:126)X, (cid:126)X m ))) ∂ F T ∂ F T . (19)Positive definiteness of the modified Lagrangian, from which all constraint variables—i.e. La-grange multipliers—have been condensed out, reflects stability of the combined multiscalesystem including micro- as well as macro-quantities. Upon condensing out also the microfluc-tuation field (cid:126)w , stability of the macroscopic system can be assessed. The microfluctuationfield as well as all Lagrange multipliers can be eliminated for each macroscopic point (cid:126)X bystatic condensation. This is done by means of the Schur complement after the FE numericaldiscretization, as described in Section 3 below, whereas stability of the macroscopic systemis detailed in Section 3.3.Substituting Eq. (12) into the first expression on the right hand side of Eq. (18), oneobtains a set of p ( p + 1) / p = n + 2, coupled specific stiffness quantities for the macro-scopic, (cid:126)v , v i , i = 1 , . . . , n , and microscopic, (cid:126)w , fields. In particular, the individual termscorresponding to the macroscopic quantities read as D = (cid:10) D (cid:11) Ω m , (20)11 i N = D i N = (cid:10) D : (cid:126) ∇ m (cid:126)ϕ i (cid:11) Ω m , (21) D i B = (cid:0) D LT i B (cid:1) RT = (cid:10) D RT · (cid:126)ϕ i + [( D : (cid:126) ∇ m (cid:126)ϕ i ) (cid:126)X m ] (cid:11) Ω m , (22) D i N j N = (cid:10) (cid:126) ∇ m (cid:126)ϕ i : D : (cid:126) ∇ m (cid:126)ϕ j (cid:11) Ω m , (23) (cid:126)D i N j B = (cid:126)D i B j N = (cid:10) (cid:126) ∇ m (cid:126)ϕ i : [ D RT · (cid:126)ϕ j + ( D : (cid:126) ∇ m (cid:126)ϕ j ) (cid:126)X m ] (cid:11) Ω m , (24) D i B j B = (cid:10) (cid:126)ϕ i · D RT · (cid:126)ϕ j + ( (cid:126)ϕ i · D : (cid:126) ∇ m (cid:126)ϕ j ) (cid:126)X m + (cid:126)X m ( (cid:126) ∇ m (cid:126)ϕ i : D RT · (cid:126)ϕ j ) + (cid:126)X m ( (cid:126) ∇ m (cid:126)ϕ i : D : (cid:126) ∇ m (cid:126)ϕ j ) (cid:126)X m (cid:11) Ω m , (25)where the subscript 0 relates to (cid:126) ∇ δ(cid:126)v , i B to (cid:126) ∇ δ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 assessmentof the macroscopic system.
3. Numerical Implementation
At every macroscopic integration point, the macroscopic quantities (cid:126) ∇ (cid:126)v , v i , and (cid:126) ∇ v i , aresampled and passed down to the microscale, cf. Fig. 2, where the modes (cid:126)ϕ i and microfluctu-ation field (cid:126)w are defined. Taking into account the constraints (3)–(6), the microfluctuationfield (cid:126)w is computed by solving the microscale boundary value problem defined in Eq. (13).Knowing (cid:126)w , the homogenized macroscopic stresses and stiffnesses are computed followingEqs. (15)–(18) and (20)–(25). All of these quantities are required for the solution of themacroscopic boundary value problem of Eq. (14) using the standard Newton algorithm, lead-ing to a fast quadratic convergence and allowing for a bifurcation analysis. Although multipleapproaches 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 solutionof the microscopic problem is presented in Section 3.1. The macroscopic problem is subse-quently 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 atGitLab’s repository mm4mm. Using standard FE procedures, the microfluctuation field (cid:126)w and its gradient (cid:126) ∇ (cid:126)w areexpressed over each microscopic element e within an RVE in terms of the shape functionsand nodal values as (cid:126)w ( (cid:126)X, (cid:126)X m ) ≈ N ew ( (cid:126)X m ) w e ( (cid:126)X ) , (cid:126) ∇ m (cid:126)w ( (cid:126)X, (cid:126)X m ) ≈ B ew ( (cid:126)X m ) w e ( (cid:126)X ) , (26)where w e is a column of element nodal values of the (cid:126)w field, collected for all n e elements ina column matrix w , N ew is a matrix of the corresponding shape functions, and B ew a matrixof the shape function derivatives. Corresponding variations δ (cid:126)w and (cid:126) ∇ δ (cid:126)w are discretized inthe same way. The patterning modes are expressed similarly as (cid:126)ϕ i ( (cid:126)X m ) ≈ N ew ( (cid:126)X m ) ϕ ei , (cid:126) ∇ (cid:126)ϕ i ( (cid:126)X m ) ≈ B ew ( (cid:126)X m ) ϕ ei , (27)12here ϕ ei is a column storing element nodal values of (cid:126)ϕ i , collected over all elements in acolumn matrix ϕ i . Although analytical expressions for the patterning modes (cid:126)ϕ i are pro-vided below in Eqs. (56) and (62), we opted here for a discretized version for convenienceand generality, since alternative definitions of patterning modes may be more appropriate,e.g. based on linearized buckling analysis or true deformed shapes obtained numerically. Thediscretized version of the orthogonality constraint (3) takes the form (cid:10) (cid:126)w ( (cid:126)X, (cid:126)X m ) (cid:11) Ω m ≈ (cid:40) T M w = 0 , T M w = 0 , (28)where M = A n e e =1 (cid:104) ( N ew ) T N ew (cid:105) Ω e m is the symmetric Gramian matrix of microscopic shape func-tions, 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 (cid:126)w . A denotes the standard FE assembly operator, and the integration over each microcopicelement volume Ω e m is performed using a standard Gauss integration rule, for brevity notexpressed explicitly as a sum. Analogously, the discrete forms of the scalar constraints (4)read (cid:10) (cid:126)w ( (cid:126)X, (cid:126)X m ) · (cid:126)ϕ i ( (cid:126)X m ) (cid:11) Ω m ≈ ϕ T i M w = 0 , i = 1 , . . . , n, (29)whereas the set of vector constraints (5) is discretized as (cid:10) (cid:126)w ( (cid:126)X, (cid:126)X m ) · [ (cid:126)ϕ i ( (cid:126)X m ) (cid:126)X m ] (cid:11) Ω m ≈ (cid:40) q T i, M w = 0 ,q T i, M w = 0 , i = 1 , . . . , n, (30)where q i,k stores ( (cid:126)ϕ i · (cid:126)e ) X k and ( (cid:126)ϕ i · (cid:126)e ) X k at the positions of the DOFs corresponding tothe first and second component. The periodicity constraint of Eq. (6) is expressed with thehelp of the link topology matrix C pbc , described e.g. in (Miehe and Bayreuther, 2007), as (cid:74) (cid:126)w ( (cid:126)X, (cid:126)X m ) (cid:75) ≈ C pbc w = 0 . (31)Finally, the microscopic governing equation (13) is solved iteratively using the standardNewton method (see Bonnans et al., 2006, Section 14) for the linear system D ww C T pbc M . . . M q n, C pbc . . . T M ... . . . ...... q T n, M . . . dwλµνη = − f w , (32)or in a compact form as D (cid:63)ww dw (cid:63) = − f (cid:63)w . (33)In Eqs. (32) and (33), D ww = A n e e =1 (cid:10) ( B ew ) T DB ew (cid:11) Ω e m , (34)13s the microscopic stiffness matrix, f (cid:63)w is the column of microscopic nodal internal forces f w = A n e e =1 (cid:104) ( B ew ) T P (cid:105) Ω e m and unbalanced equality constraints, where P is a column storing thecomponents of the stress tensor P , and dw (cid:63) is a column storing the iterative change of themicroscopic fluctuation field and current iterative values of Lagrange multipliers. Note that,in analogy to η and ν in Eq. (7), µ stores individual components of (cid:126)µ . An asterisk indicatesthat the matrix accounts for all the equality constraints and thus also has the correspondingextra 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 theentire RVE from the knowledge of the current state of the macroscopic quantities (cid:126) ∇ (cid:126)v , v i , (cid:126) ∇ v i , patterning modes (cid:126)ϕ i , i = 1 , . . . , n , and iterative state of the microfluctuation field w .Because the primary buckling modes are captured by the patterning fields (cid:126)ϕ i , a microscopicbifurcation analysis and stability control analogous to Section 3.3 below is not required. The macroscopic fields (cid:126)v and v i are discretized within each macroscopic element E as (cid:126)v ( (cid:126)X ) ≈ N E ( (cid:126)X ) v E , v i ( (cid:126)X ) ≈ N Ei ( (cid:126)X ) v Ei ,(cid:126) ∇ (cid:126)v ( (cid:126)X ) ≈ B E ( (cid:126)X ) v E , (cid:126) ∇ v i ( (cid:126)X ) ≈ B Ei ( (cid:126)X ) v Ei , (35)where N E • and B E • are macroscopic element shape functions and v E • the corresponding columnmatrices of DOFs, collected globally for all n E elements in column matrices v and v i ; thesame forms and expressions are used to discretize their variations. The internal elementforces f E M = f E f E ... f En (36)are then obtained as a sum over n g macroscopic Gauss quadrature points, expressed explicitlyas f E = n g (cid:88) i g =1 f E,i g = n g (cid:88) i g =1 w i g J i g ( B E ) T Θ i g , (37) f Ei = n g (cid:88) i g =1 f E,i g i = n g (cid:88) i g =1 w i g J i g (cid:8) ( B Ei ) T Λ i g i − ( N Ei ) T Π i g i (cid:9) , i = 1 , . . . , n, (38)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 i g are integration weights with the corresponding Jaco-bians J i g , and the column matrices Θ i g , Λ i g i , and scalars Π i g i , store the components of thehomogenized stress quantities, defined in Eqs. (15)–(17), evaluated at appropriate positionsof associated integration points. 14n order to condense out the effect of the constrained microfluctuation field w and allLagrange multipliers, the following monolithic incremental system of equations is assembledat each macroscopic quadrature point i g of each element E , including both macroscopicas well as microscopic quantities (superscripts E and i g are dropped in Eqs. (39)–(49) forbrevity): K K · · · K n K (cid:63) w K K · · · K n K (cid:63) w ... ... . . . ... ... K n K n · · · K nn K (cid:63)nw K (cid:63)w K (cid:63)w · · · K (cid:63)wn K (cid:63)ww v E v E ... v En w (cid:63) + f f ... f n f (cid:63)w = r r ... r n . (39)In Eq. (39), r i , i = 0 , . . . , n , are the residuals reflecting the fact that equilibrium is not satis-fied at the level of an integration point, but only for the entire assembly over all quadraturepoints of all elements. The last row is zero, since the microscopic system of Eq. (33) has beenequilibrated, implying also that f (cid:63)w = 0. The specific stiffness for the macroscopic quantitiescan be then obtained as a Schur complement via static condensation: K E,i g M = K K · · · K n K K · · · K n ... ... . . . ... K n K n · · · K nn − K (cid:63) w K (cid:63) w ... K (cid:63)nw (cid:2) K (cid:63)ww (cid:3) − (cid:2) K (cid:63)w K (cid:63)w · · · K (cid:63)wn (cid:3) . (40)Note that whereas for a given element E an n g microfluctuation fields w (cid:63),E,i g are com-puted and condensed out, only one set of macroscopic DOFs for the coarse fields v E and v Ei pertinent to that element are involved. The asterisk superscript to the specific stiffness com-ponents in Eqs. (39) and (40) again indicates that these matrices have extra zero entriescorresponding to the Lagrange multipliers. The individual sub-matrices are defined as K = 1 | Ω m | B T D B , (41) K i = K T i = 1 | Ω m | (cid:0) B T D i B B i + B T D i N N i (cid:1) , (42) K ij = 1 | Ω m | (cid:0) N T i D i N j N N j + B T i D i B j N N j + N T i D T i B j N B j + B T i D i B j B B j (cid:1) , (43) K w = K T w = 1 | Ω m | B T D w , (44) K iw = K T wi = 1 | Ω m | (cid:0) N T i D i N w + B T i D i B w (cid:1) , (45) K (cid:63)ww = 1 | Ω m | D (cid:63)ww , (46)where D , D i N = D T i N , D i B = D T i B , D i N j N , D i N j B , and D i B j B , are matrix representations ofthe specific stiffnesses defined in Eqs. (20)–(25), obtained upon RVE discretization simply15s volume integrals using a standard Gauss integration rule. For instance, D is a 4 × D , etc. In additionto the expressions (41)–(46), there are 3 n + 1 coupled specific stiffnesses related to themicroscopic variation (cid:126) ∇ δ (cid:126)w and one of the macroscopic variations (cid:126) ∇ δ(cid:126)v , δv i , or (cid:126) ∇ δv i , andone microscopic stiffness quantity that relates twice to (cid:126) ∇ δ (cid:126)w . Because (cid:126) ∇ δ (cid:126)w depends on themicroscopic position (cid:126)X m , over which the integral in Eq. (18) is carried out, these specificstiffnesses can be written only upon RVE discretization, yielding D w = D T w = A n e e =1 (cid:10) DB ew (cid:11) Ω e m , (47) D i N w = D T wi N = A n e e =1 (cid:10) ( B ew ϕ ei ) T DB ew (cid:11) Ω e m , (48) D i B w = D T wi B = A n e e =1 (cid:10)(cid:8) ( N ew ϕ ei ) T D + X m ( B ew ϕ ei ) T D (cid:9) B ew (cid:11) Ω e m , (49)where the integration over Ω e m is again carried out numerically.The stiffness matrix of the entire macroscopic element is obtained by summing the con-tributions from all quadrature points, K E M = n g (cid:88) i g =1 w i g J i g K E,i g M , (50)which are eventually assembled into a global stiffness matrix K M = A n E E =1 K E M . The resultingmacroscopic system is again solved using the standard Newton method with an incrementalsystem of linear equations K M dv M = f ext − f M , (51)where f ext denotes a column of externally applied forces (acting only on v ), and dv M = dv dv ... dv n (52)is an iterative increment of the global macroscopic quantities. Following Miehe and Koch (2002), an equilibrated configuration v M of a system is con-sidered to be stable if the energy of this state is lower than the energy associated with a stateobtained by adding a small kinematically admissible perturbation δv M to the equilibratedconfiguration. That is, if (cid:98) E ( v M + δv M ) − (cid:98) E ( v M ) ≈ δv T M K M ( v M ) δv M > , (53)where the second-order Taylor series expansion of the total energy has been used. Noticethat (cid:98) E corresponds to the total potential energy of the entire system (Eq. (8)) from which16icrofluctuation fields and Lagrange multipliers, w (cid:63) , associated with all macroscopic Gausspoints have been condensed out, and that the first-order term vanishes as a result of equi-librium. 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 α isnon-positive, the associated configuration is unstable and the corresponding eigenvector ψ determines the buckling mode. The equilibrated solution v M of the current increment is thenperturbed with the eigenvector ψ multiplied by a small perturbation factor τ > v M = v M + τ ψ , (54)and the system is equilibrated again. The factor τ is increased until a stable equilibriumis reached, i.e. until a possible energy barrier between the current unstable and a stablebuckled configuration is overcome, and at the same time until the lowest eigenvalue α ofthe updated macroscopic stiffness matrix K M does not become positive. If the system failsto find a stable equilibrium even for large τ , the previous increment is halved to decreasethe energy barrier, and the entire procedure is repeated.An outline of the overall micromorphic computational homogenization scheme for multi-ple modes, including stability control, is given in Algorithm 1.
4. Numerical Examples
In this section, predictions made using the micromorphic computational homogenizationscheme, introduced in Sections 2 and 3, are compared against Direct Numerical Simula-tions (DNS) for two examples. The first example represents a metamaterial column com-posed of a square stacking of holes subjected to compression, whereas the second exampleconsiders a specimen with a hexagonal stacking of holes subjected to a uniform compressiveloading with various biaxiality ratios, buckling locally into one of multiple possible patterns.The constitutive behaviour of the elastomer base material is modelled by a hyperelasticlaw with the following energy density ψ ( F ( (cid:126)X, (cid:126)X m )) = c ( I −
3) + c ( I − − c log J + 12 K ( J − , (55)where F = I + ( (cid:126) ∇ (cid:126)u ) T is the deformation gradient tensor, where the gradient operator (cid:126) ∇ is defined with respect to the reference configuration, J = det F , and I = tr C is thefirst invariant of the right Cauchy–Green deformation tensor C = F T · F . The valuesof the constitutive parameters employed, listed in Tab. 1, are based on the experimentalcharacterization of Bertoldi et al. (2008).The smallest RVE domains of the size 2 (cid:96) are adopted in both examples for the micromor-phic homogenization scheme, as shown in Figs. 4b and 10b below, which are large enoughto accommodate the longest microstructural patterning modes (see e.g. Bertoldi et al. 2008for the square and Ohno et al. 2002a for hexagonal stacking of holes). Although the cho-sen RVE size is sufficiently large to accommodate microstructural buckling, because of theperiodicity assumption on the microfluctuation fields (cid:126)w (recall Eq. (6) and the discussiontherein), choosing larger RVE domains may still slightly affect obtained results, especiallyfor a vanishing separation of scales. 17 lgorithm 1: Nested solution scheme for the micromorphic computational homogenization framework witha full Newton implementation and stability control. Initialization: (a): Initialize the macroscopic model, v ( t = 0) = 0, v i ( t = 0) = 0 for all i = 1 , . . . , n .(b): Assign an RVE to each Gauss integration point of the macro-model.2: for k = 1 , . . . , n T (loop over all time steps of associated parametrization time)(i): Apply macroscopic boundary conditions at time step k .(ii): while (cid:15) > tol (macroscopic solver, iteration l )(a): From v l and v li compute deformation gradient I +( (cid:126) ∇ (cid:126)v i g ) T , mode magnitude v i g i ,and its gradient (cid:126) ∇ v i g i for each macroscopic Gauss point i g .(b): Perform the RVE analysis for each macroscopic Gauss point i g : - Apply underlying deformation dictated by I + ( (cid:126) ∇ (cid:126)v i g ) T , v i g i , (cid:126) ∇ v i g i , and ϕ i . - Assemble and solve the nonlinear RVE problem, Eq. (13). For w enforce or-thogonality, periodicity, and rigid body motion constraints (3)–(6) over Ω m . - Average resulting microscopic quantities to obtain the homogenized macro-scopic stresses and stiffnesses.(c): Assemble the macroscopic gradient f l M and tangent K l M by condensing out thestiffness terms related to w and all Lagrange multipliers.(d): Update the macroscopic displacements v l +1M = v l M + dv l M , where K l M dv l M = f ext − f l M .(e): Update the iteration error (cid:15) = (cid:107) f l M (cid:107) + (cid:107) dv l M (cid:107) .(iii): end while (iv): If the lowest eigenvalue α of K l M is non-positive, perturb the system with corre-sponding eigenvector τ ψ , and equilibrate iteratively for an increasing perturbationfactor τ until the system becomes stable. Then continue to (i) for k = k + 1. Ifperturbation fails, halve the load increment and proceed to (i) with current k .3: end for Table 1: Constitutive parameters of the hyperelastic law specified in Eq. (55), used in both numericalexamples.
Parameter c c K [MPa] [MPa] [MPa]Value 0 .
55 0 . .1. Example 1: Local Versus Global Buckling The first example analyses a finite column of width W and height H with a microstructureconsisting of a square stacking of unit cells with edge size (cid:96) = 9 .
97 mm and circular holes ofdiameter d = 8 .
67 mm, cf. Fig. 4a. The bottom and top edges of the specimen domain aredisplaced by ± u/ (cid:126)e to induce 10% overall compressive strain, defined as u/H . Dependingon the slenderness ratio H/W , a competition between microstructural buckling (patterntransformation) and macrostructural buckling of the structure is expected. A similar examplehas been investigated numerically as well as experimentally by Coulais et al. (2015).For the DNS solutions, the entire domain is discretized using isoparametric quadratictriangular elements of typical size h m = (cid:96)/
10 with three Gauss integration points, as shownin Fig. 4b. For this case, a single microstructural realization adequately represents theensemble averaged DNS solution. This is shown in Fig. 5a, where nominal stress–straindiagrams corresponding to 100 microstructural translations (all possible combinations of 10translation steps in horizontal and 10 in vertical direction, covering together one period ofthe microstructure (cid:96) × (cid:96) ) are shown for a 6 (cid:96) × (cid:96) specimen. The overall response is initiallylinear until the first bifurcation point is reached, upon which a local patterning emergesand the specimen’s stiffness drops close to zero. Further increasing the compressive strainleads to the second bifurcation, corresponding to global buckling of the specimen, uponwhich the overall stiffness becomes negative. The maximum and minimum envelopes of allrealizations deviate less than 6% from the corresponding mean, suggesting that the referencemicrostructure 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ˇset al., 2019, Eq. (7)) (cid:126)ϕ ( (cid:126)X ) = 1 C (cid:104) − sin π(cid:96) ( X + X ) − sin π(cid:96) ( − X + X ) (cid:105) (cid:126)e + 1 C (cid:104) sin π(cid:96) ( X + X ) − sin π(cid:96) ( − X + X ) (cid:105) (cid:126)e , (56)were C is a normalization constant ensuring that | Q | (cid:82) Q (cid:107) (cid:126)ϕ ( (cid:126)X ) (cid:107) d (cid:126)X = 1, and where Q isthe 2 (cid:96) × (cid:96) periodic cell of Fig. 3. For the RVE discretization, the same type and density ofelements is used as for the DNS system (Fig. 4b). To provide sufficient kinematic freedom tothe macroscopic micromorphic system, a mesh convergence study is performed. A uniformmacroscopic mesh of quadratic isoparametric triangular elements with three Gauss points isconsidered with characteristic element sizes h M ∈ { , . . . , } (cid:96) . The same mesh is consideredfor both macroscopic fields (cid:126)v and v . An example of a particular mesh of an elementsize h M = 2 (cid:96) for a 4 (cid:96) × (cid:96) specimen (for which local buckling is expected to occur) is shownin Fig. 4c. The obtained results in terms of nominal stress–strain diagrams are plotted inFig. 5b. For element sizes h M ≤ (cid:96) , the behaviour is similar to the DNS result of Fig. 5a;moreover, the results of element sizes h M = 2 (cid:96) and h M = 1 (cid:96) are indistinguishable. Theelement size h M = 2 (cid:96) is thus adopted in what follows, although from the deformed shapeof Fig. 6d it may be clear that a locally refined mesh might be useful. For more details on19ppropriate choice of element types and associated integration rules see Rokoˇs et al. (2020b). H W u u Ω (cid:126)e (cid:126)e (cid:96)d (a) specimen geometry (b) micro-mesh (c) macro-meshFigure 4: (a) Specimen geometry for the local versus global buckling example, i.e. a specimen of a width W =4 (cid:96) and height H = 8 (cid:96) made of an elastomeric mechanical metamaterial, subjected to an overall verticalcompressive strain of 10%. (b) Typical discretization of the full Direct Numerical Simulation (DNS) andRVE model, element size h m = (cid:96)/
10. (c) Macroscopic discretization of a 4 (cid:96) × (cid:96) specimen with elementsize h M = 2 (cid:96) . In both cases, isoparametric quadratic triangular elements with three Gauss points are used. Depending on the slenderness ratio
H/W , two basic and mutually interacting deformationmechanisms occur, as shown in Fig. 6 by the deformed configurations for 6 (cid:96) × (cid:96) and 6 (cid:96) × (cid:96) specimens. The first mechanism is local patterning (Fig. 6a), emerging for low slendernessratios upon reaching a critical compressive strain of approximately 3%. The cells fold in atypical pattern of alternating ellipsoidal holes and an auxetic effect is observed along withboundary layers where the local buckling is restricted. The second mechanism, occurring forhigher slenderness ratios, is global buckling (Fig. 6c), which is triggered upon reaching thecritical Euler buckling stress. The corresponding buckling strain can be estimated as ε cr ( H/W ) = π H/W ) . (57)In Figs. 6b and 6d, the micromorphic field normalized by its maximum value consideredover space and a parametrization pseudo-time t , i.e. (cid:98) v = v / (cid:107) v ( t, (cid:126)X ) (cid:107) ∞ , is shown in colour.Comparing Fig. 6a with 6b, and Fig. 6c with 6d, we conclude that the micromorphic ho-mogenization scheme is capable of accurately reconstructing the overall kinematic response,correctly capturing the auxetic effect for lower slenderness ratios, and accurately indicatingregions of localized patterning reflected by the magnitude of the micromorphic field (cid:98) v forlarger slenderness ratios. The pattering regions localize in the compressive parts of the bentdomain, situated close to the specimen’s centre and near the supports at both ends. Theoverall deformed shape for the large slenderness ratio, i.e. the (cid:126)v field, is captured with goodaccuracy as well. 20 a) DNS shifts (b) MM mesh convergence studyFigure 5: (a) Stress–strain diagrams for 100 translated DNS solutions (10 translation steps in horizontaland 10 in vertical direction, covering one period of the microstructure (cid:96) × (cid:96) ), W = 6 (cid:96) and H = 30 (cid:96) .(b) Macroscopic mesh convergence study for the micromorphic (MM) homogenization scheme with a uniformtriangulation, cf. Fig. 4c. Considered element sizes h ∈ { , . . . , } (cid:96) for W = 6 (cid:96) and H = 30 (cid:96) (slendernessratio H/W = 5).
Nominal stress–strain diagrams for W = 6 (cid:96) and slenderness ratios H/W ∈ [1 ,
30] areshown in Fig. 7. Here, the two mechanisms and their mutual interactions are visible moreclearly. For the DNS results (Fig. 7a) and the applied strain range u/H ∈ [0 , .
1] the slen-derness ratios up to
H/W = 2 .
33 buckle only locally, ratios
H/W ∈ [2 . , .
67] buckle firstlocally and then globally, the ratio
H/W = 7 buckles locally and globally at the same time,whereas ratios
H/W ∈ [7 . ,
30] buckle first globally and then locally. Bilinear stress–strainresponses typically emerge, which exhibit softening in later stages due to the presence ofsecondary (local or global) buckling, see also Fig. 5a and the discussion therein. A close-upon the local versus global buckling intersection is shown in Fig. 8a, where mild snap-backsfor slenderness ratios
H/W ∈ [6 . , .
33] can be observed. The micromorphic computationalhomogenization is capable of reconstructing the stability behaviour (Fig. 7b) with an accu-racy that decreases with decreasing scale ratio (i.e. larger errors are observed for smaller scaleratios). In particular, for
H/W = 1 we see a large discrepancy in the post-bifurcation nomi-nal stiffness (i.e. in the slope of the P versus u/H curve), which corresponds to 25% of errorrelative to the initial pre-bifurcation stiffness (which is practically constant for all consid-ered H/W ratios). With increasing slenderness ratio, however, the error drops rapidly downto 0 . H/W and corre-sponding to the first bifurcation points of Fig. 7 are shown in Fig. 9 for several specimenwidths W ∈ { , , } (cid:96) . Bertoldi et al. (2008) reported that the buckling strain of a single21 a) DNS, 6 (cid:96) × (cid:96) (b) MM, 6 (cid:96) × (cid:96) (c) DNS, 6 (cid:96) × (cid:96) (d) MM, 6 (cid:96) × (cid:96) Figure 6: Comparison of deformed shapes obtained from Direct Numerical Simulation (DNS) and the micro-morphic (MM) computational homogenization scheme. (a) DNS and (b) MM results for a 6 (cid:96) × (cid:96) specimen(slenderness ratio H/W = 2 . (cid:107) v ( t, (cid:126)X ) (cid:107) ∞ = 51 . (cid:96) × (cid:96) specimen( H/W = 5 . (cid:107) v ( t, (cid:126)X ) (cid:107) ∞ = 69 . v field normalized by itsextreme value in time and space, i.e. (cid:98) v . The results are shown for an overall applied strain u/H = 0 . P and Θ as a function of the applied nominal strain u/H for columnsof various slenderness ratios H/W ∈ [1 , W = 6 (cid:96) , shown in colour, obtained via (a) Direct NumericalSimulation (DNS), and (b) micromorphic (MM) computational homogenization. Instability points corre-sponding to global buckling are connected by the black dash-dot lines ( ), whereas local buckling pointsare connected by the black dashed lines ( ); cf. also Fig. 8. a) DNS stress–strain curves, close-up (b) bifurcation curvesFigure 8: (a) A close-up on stress–strain curves corresponding to the DNS solutions of Fig. 7a in the vicinityof the local versus global buckling intersection. (b) Bifurcation curves for direct numerical solutions (blue,cf. Fig. 7a) and micromorphic homogenization (red, cf. Fig. 7b). Instability points corresponding to globalbuckling are connected by the dash-dot lines ( ), whereas local buckling points are connected by thedashed lines ( ). RVE corresponds to approximately 3% (shown as the black dash-dot line in Fig. 9), and hencethis value is expected to be the theoretical local buckling strain. Both the DNS as well as mi-cromorphic results attain this limit for the range of slenderness ratios
H/W ∈ [1 , H/W = 7,local and global buckling occur simultaneously, whereas higher slenderness ratios convergeasymptotically towards the theoretical bound of the global Euler buckling strain given byEq. (57) (shown as the black dashed line). Both the DNS and the micromorphic scheme ap-proach this limit from below, although the micromorphic results overestimate systematicallythe 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 lo-cal buckling strain converges towards the theoretical bound 3%, whereas the global bucklingstrain shows little change. The overall relative error of the micromorphic scheme in terms ofthe buckling strain does not exceed 3% for local and 6% for global buckling. With increasingscale ratio, this error drops down to 0 .
5% for local and 3% for global buckling. Note thatthe theoretical global Euler buckling strain given by Eq. (57) significantly overestimates theDNS results due to a substantial amount of shear and changes triggered in the microstructureupon buckling (cf. Fig. 6c), which become important especially for intermediate slendernessratios
H/W ∈ (7 , igure 9: Buckling strain u/H as a function of the slenderness ratio H/W ∈ [1 , W ∈ { , , } (cid:96) , obtainedfrom the DNS and micromorphic (MM) computational homogenization algorithms, and from the theoreticalestimates for local ( ) and global ( ) buckling. The second example considers an infinite microstructure composed of a hexagonal stack-ing of holes, shown in Fig. 10a. The periodic cell, considered later as RVE for the micromor-phic 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 biaxialityratio γ = | F − || F − | = | ε || ε | ∈ [0 , ∞ ] , (58)where ε ii = F ii − (cid:126)e i , i = 1 ,
2, directions,and F ij are the components of the overall deformation gradient tensor F , with F = F =0. The three distinct patterns, depicted in Figs. 11a–11c, correspond to the following cases:(i) Pattern I , uniaxial or shear pattern denoted (cid:126)π , occurs when the compressive load onthe vertical cell walls is higher than the load on the other cell walls, i.e. γ >
1. Themultiplicity of the bifurcation point corresponds to one. The displacement field leads to theformation of horizontal layers of holes sheared alternatingly to the right and to the left, seeFig. 11a. (ii)
Pattern II , also called biaxial or butterfly-like pattern and denoted (cid:126)π , emergeswhen the inclined cell walls at θ = ± ◦ are compressed more than the vertical cell walls,i.e. γ <
1. In this case, the multiplicity of the bifurcation point is two and the patternexhibits horizontal layers of holes buckled along the horizontal and vertical directions, seeFig. 11b. (iii)
Pattern III , also referred to as the equi-biaxial or flower-like pattern (cid:126)π , isobserved when all three cell walls are subjected to an equal compressive load, i.e. γ = 1. Themultiplicity of the bifurcation point equals three, and the displacement field corresponds toa virtually undeformed central hole, surrounded by ellipses, see Fig. 11c.24 a) hexagonal stacking P P P (cid:126)e θ(cid:96) d (cid:126)e (b) hexagonal RVE (c) macroscopic meshFigure 10: (a) A microstructure with a hexagonal stacking of holes in the reference configuration. (b) Asingle hexagonal RVE with an average mesh size h m = (cid:96)/
10. The boundaries indicated by the same colourare coupled via the periodicity constraint (6). (c) A macroscopic periodic mesh for the micromorphiccomputational homogenization scheme, h M = 3 . (cid:96) . Because the individual patterns (cid:126)π i are not mutually orthogonal, it is convenient for furthertreatment to introduce the so-called modes (cid:126)ϕ i , i = 1 , ,
3, (see Ohno et al., 2002b; Okumuraet al., 2002; Rokoˇs et al., 2020a), which satisfy orthogonality. Linear combinations of thesemodes result in the previously introduced patterns as follows: (cid:126)π = (cid:126)ϕ , (59) (cid:126)π = (cid:126)ϕ + (cid:126)ϕ , (60) (cid:126)π = (cid:126)ϕ + (cid:126)ϕ + (cid:126)ϕ . (61)The individual modes (cid:126)ϕ i , i = 1 , ,
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 ± ◦ , whichcan be expressed in an analytical form. The first mode reads (see Rokoˇs et al., 2020a, Eq. (3)) (cid:126)ϕ ( (cid:126)X ) = 1 C (cid:20) sin (cid:18) πX √ (cid:96) (cid:19) (cid:126)e + 1 √ (cid:18) πX (cid:96) (cid:19) (cid:126)e (cid:21) , (62)were C is a normalization constant ensuring that | Q | (cid:82) Q (cid:107) (cid:126)ϕ ( (cid:126)X ) (cid:107) d (cid:126)X = 1 for the periodiccell Q of Fig. 10b, whereas modes II and III are obtained by rotating (cid:126)ϕ by ∓ ◦ , see Figs. 11eand 11f.The example analysed in this section represents an infinite specimen, made of a hexag-onal cellular structure with a hole diameter d = 1 .
28 mm and a centre-to-centre spac-ing (cid:96) = 1 .
386 mm, subjected to biaxial compression. For the micromorphic computationalhomogenization it is modelled with a 10 ×
10 mm periodic square domain discretized witheight identical quadratic triangular elements of size h m = 3 . (cid:96) with a three-point Gaussintegration rule, see Fig. 10c. Again, the same macroscopic discretization is used for allthree micromorphic fields v i , i = 1 , ,
3, as well as for the mean solution (cid:126)v . The RVE,shown in Fig. 10b, discretized with isoparametric quadratic triangular elements of averagesize h m = (cid:96)/
10 using a three-point Gauss integration rule, is assigned to each macroscopic25 a) pattern I (b) pattern II (c) pattern III(d) mode I (e) mode II (f) mode IIIFigure 11: Three pattern transformations (cid:126)π i and orthogonal modes (cid:126)ϕ i . (a) Pattern I, (cid:126)π , the uniaxial orshear pattern, corresponds to γ >
1. (b) Pattern II, (cid:126)π , the biaxial or butterfly pattern, occurs for γ < (cid:126)π , the equi-biaxial pattern or flower-like pattern, emerges for γ = 1. (d) Mode I, (cid:126)ϕ ,developing perpendicular to θ = 90 ◦ , (e) mode II, (cid:126)ϕ , developing perpendicular to θ = 30 ◦ , and (f) mode III, (cid:126)ϕ , developing perpendicular to θ = − ◦ . integration point. The three orthogonal modes (cid:126)ϕ i from Figs. 11d–11f are considered in theansatz in Eqs. (1) and (2), i.e. n = 3. A similar preliminary analysis of this case has beenreported in Rokoˇs 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 solverand a bifurcation analysis are used instead.Because an infinite specimen is considered, the DNS solution directly corresponds to thebehaviour of a single periodic cell (i.e. RVE) subjected to a compressive load of biaxiality ra-tio γ . Since the deformation state is periodic, the ensemble average reduces to volume averag-ing, easily obtained from a single microstructural translation. The question arises, however,whether the micromorphic computational homogenization is capable of reproducing such abehaviour, i.e. yielding an affine mean field (cid:126)v and constant micromorphic fields (cid:98) v i , while pro-viding patterns of Eqs. (59)–(61) as the outcome of the analysis. Fig. 12 collects the resultsfor a uniaxial compressive strain ε = − .
05 ( γ = ∞ ). As expected, the mean displacementin the (cid:126)e direction is zero, whereas in the (cid:126)e direction it is linear and corresponds to theapplied nominal strain. The normalized micromorphic fields (cid:98) v i = v i / max k =1 , , (cid:107) v k ( t, (cid:126)X ) (cid:107) ∞ are also spatially constant with (cid:98) v = 1 and (cid:98) v = (cid:98) v = 0, resulting in an activation of mode Iand, consequently, pattern I (recall Eq. (59)). Moreover, the deformed RVE shape in Fig. 12a26atches the DNS solution in Fig. 11a, corroborating further the validity of the micromorphicresults.The evolution of the magnitudes corresponding to the individual micromorphic fields asa function of ε is shown for the overall applied deformation gradient F ( t ) = (1 + ε ) (cid:126)e (cid:126)e + (1 + ε ) (cid:126)e (cid:126)e (63)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 thebifurcation point where a negative or sufficiently small lowest eigenvalue is observed andthe system is perturbed towards the corresponding eigenvector. The correct patterns aretriggered, i.e. (cid:98) v (cid:54) = 0 while (cid:98) v = (cid:98) v = 0 for pattern I ( γ = ∞ ), (cid:98) v = (cid:98) v (cid:54) = 0 while (cid:98) v = 0for pattern II ( γ = ), and (cid:98) v = (cid:98) v = (cid:98) v (cid:54) = 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 minimais explored. To this end, all micromorphic fields are initialized as constant fields, withmagnitudes spanning the entire cube [ (cid:98) v , (cid:98) v , (cid:98) v ] ∈ [0 , × [0 , × [0 ,
1] for a fixed appliedoverall strain which corresponds to a buckled state, while assuming the exact mean fields (cid:126)v .It is found that although other patterns may yield stable local minima, the global minimaalways correspond to the correct combinations of modes. Note that similarly to the DNS, themultiplicities of the bifurcation points associated with the second and third pattern occuralso for the micromorphic formulation. In that case, the associated buckling modes havea zero mean (cid:126)v = (cid:126) ε , ε ] ∈ [0 , − . × [0 , − . (cid:98) v i are shown. Four regions A–C are distinguished, as depicted in the correspondingcontour plot in Fig. 14b. In the first region, A, no pattern is triggered, because the criticalbifurcation strain has not been exceeded yet. Pattern I occurs in region B, whereas thesecond pattern is triggered in region C. In the fourth region, D, a mixture of both patternsis observed, with the special configuration for γ = 1 corresponding to pattern III (denotedby the black solid line). A similar behaviour has been observed for hexagonal honeycomb27 a) deformed RVE (b) (cid:126)e component of (cid:126)v [mm] (c) (cid:126)e component of (cid:126)v [mm](d) (cid:98) v [-] (e) (cid:98) v [-] (f) (cid:98) v [-]Figure 12: (a) A single deformed RVE and the resulting macroscopic fields of a homogenized infinite specimensubjected to a uniaxial compressive strain of ε = 0 .
05 ( γ = ∞ ). The two components of the mean displace-ment field (cid:126)v are shown in (b) and (c), whereas the fields (cid:98) v i that indicate the relative activation of individualmodes (cid:126)ϕ i , i = 1 , ,
3, are shown in (d)–(f), where the normalization constant is max k (cid:107) v k ( t, (cid:126)X ) (cid:107) ∞ = 3 . γ = ∞ (b) γ = (c) γ = 1Figure 13: Magnitudes of all three spatially constant modes (cid:98) v i , i = 1 , ,
3, as a function of the overall appliedvertical strain ε . Three biaxiality ratios are considered: (a) γ = ∞ (pattern I, − ε ∈ [0 , . ε = 0,max k (cid:107) v k ( t, (cid:126)X ) (cid:107) ∞ = 3 . γ = (pattern II, − ε ∈ [0 , . ε = ε , max k (cid:107) v k ( t, (cid:126)X ) (cid:107) ∞ = 3 . γ = 1 (pattern III, − ε = − ε ∈ [0 , . k (cid:107) v k ( t, (cid:126)X ) (cid:107) ∞ = 3 . a) phase diagram (b) magnitude contour plot(c) circumferential sectionFigure 14: A phase diagram of individual normalized mode magnitudes (cid:98) v i plotted as a function of the overallapplied strains ε and ε . The entire surface plot is shown in (a) with the colour coding as used in (c),the corresponding contour plot is shown in (b) (normalization constant max k (cid:107) v k ( t, (cid:126)X ) (cid:107) ∞ = 3 . (cid:98) ε ii = ε ii / .
025 in (c) (normalization constant max k (cid:107) v k ( t, (cid:126)X ) (cid:107) ∞ = 2 . structures in the work of Okumura et al. (2002); see Fig. 10 therein. Fig. 14c plots a cir-cumferential section through the phase diagram of 14a taken along the dashed black curvehighlighted in Fig. 14b, including its extension to the other three strain quadrants. Themagnitudes of the individual normalized micromorphic fields (cid:98) v i are plotted as a function ofangle ϑ ∈ [0 , π ], spanning the entire circle. The angle starts from the ( ε = − , ε = 0)direction and sweeps clockwise. Again, Fig. 14c confirms that equal magnitudes of all modesoccur for ϑ = π/
4. Furthermore, it is clearly visible which strain combinations yield whichmicrostructural pattern. For instance, close to ϑ = 5 π/ π/ ε ii is positive, a pattern transformation occurs due to alarge negative magnitude of the other compressive strain.29 . Summary and Conclusions This contribution has extended a recently developed micromorphic computational ho-mogenization framework for mechanical metamaterials with a full Newton solver. The mi-cromorphic framework decomposes the kinematic field by exploiting prior knowledge on thetypical patterning modes, allowing to accurately capture non-local effects present in the mi-crostructure. The derivation and implementation of a full Newton solver for this frameworkhas been provided, including analytical expressions for the first and second variations of thetotal effective potential energy. Significant gains have been obtained compared to the ex-isting quasi-Newton implementation, in particular with respect to the bifurcation analysis,which is essential for applications of elastomeric mechanical metamaterials relying on localand global buckling. Two examples have been tested to demonstrate the capabilities of thepresented numerical scheme. In the first example, a metamaterial column consisting of amicrostructure with a square stacking of holes has been analysed for various slendernessratios, for which a competition between the local and global buckling exists. The secondexample 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 homogenizationframework is robust and efficient.2. The micromorphic approach captures the behaviour of the reference Direct NumericalSimulation (DNS) accurately in terms of both local and global buckling as well as thepattern magnitudes.3. The nominal stresses are reproduced by the micromorphic framework with a goodaccuracy, although the post-bifurcation results are in general systematically overesti-mated compared to DNS. The maximum error in terms of the critical buckling stresscorresponding to the first instability point does not exceed 12%, and decreases withincreasing scale ratio down to approximately 7%.4. The buckling strain is captured with a higher accuracy compared to the nominal buck-ling stress. The relative error stays below 3% for local and 6% for global buckling, anddecreases 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 ahexagonal stacking of holes, for which multiple patterning modes occur. It predictsthe correct patterns for the loading directions considered.The full Newton solver presented here greatly reduces the dependency of the solution onthe initial guess by perturbing the system towards the correct direction when a bifurcationpoint is encountered; therefore, it provides an indispensable numerical tool for modellinginstability-based mechanical metamaterials.
Acknowledgements
The research leading to these results has received funding from the European ResearchCouncil under the European Union’s Seventh Framework Programme (FP7/2007-2013)/ERC30rant agreement no. [339392] (O. Rokoˇs 09/2016–03/2019, R.H.J. Peerlings, and M.G.DGeers) and from the Czech Science Foundation (GA ˇCR) grant agreement no. [19-26143X](O. Rokoˇs 03/2019–12/2019, and M. Doˇsk´aˇr). The authors would like to also acknowledgeProf. Jan Zeman from the Czech Technical University in Prague for fruitful discussions andcritical comments on the manuscript, and Dr. Geralf H¨utter from TU Bergakademie Freibergfor valuable email communication.
References
Ameen, M. M., Rokoˇs, O., Peerlings, R. H. J., Geers, M. G. D., 2018. Size effects in nonlinearperiodic materials exhibiting reversible pattern transformations. Mechanics of Materials124, 55–70.URL
Bertoldi, K., Boyce, M. C., Deschanel, S., Prange, S. M., Mullin, T., 2008. Mechanicsof deformation-triggered pattern transformations and superelastic behavior in periodicelastomeric structures. Journal of the Mechanics and Physics of Solids 56 (8), 2642–2668.URL
Bonnans, J. F., Gilbert, J. C., Lemar´echal, C., Sagastiz´abal, C. A., 2006. Numerical Op-timization: Theoretical and Practical Aspects (Universitext). Springer-Verlag New York,Inc., Secaucus, NJ, USA.Coulais, C., Overvelde, J. T. B., Lubbers, L. A., Bertoldi, K., van Hecke, M., 2015. Dis-continuous Buckling of Wide Beams and Metabeams. Physical Review Letters 115 (4),044301.URL https://link.aps.org/doi/10.1103/PhysRevLett.115.044301
Forest, S., Trinh, D. K., 2011. Generalized continua and non-homogeneous boundary condi-tions in homogenisation methods. ZAMM - Journal of Applied Mathematics and Mechanics/ Zeitschrift f¨ur Angewandte Mathematik und Mechanik 91 (2), 90–109.URL http://doi.wiley.com/10.1002/zamm.201000109
Geers, M. G. D., Kouznetsova, V. G., Brekelmans, W. A. M., 2010. Multi-scale computa-tional homogenization: Trends and challenges. Journal of Computational and AppliedMathematics 234 (7), 2175–2182.URL http://dx.doi.org/10.1016/j.cam.2009.08.077https://linkinghub.elsevier.com/retrieve/pii/S0377042709005536
Kolken, H. A., Zadpoor, A. A., 2017. Auxetic mechanical metamaterials. RSC Adv. 7, 5111–5129.URL http://dx.doi.org/10.1039/C6RA27333E
Kouznetsova, V. G., Brekelmans, W. A. M., Baaijens, F. P. T., 2001. An approach to micro-macro modeling of heterogeneous materials. Computational Mechanics 27 (1), 37–48.URL https://doi.org/10.1007/s004660000212 https://doi.org/10.1016/j.cma.2003.12.073
Kunc, O., Fritzen, F., 2019. Finite Strain Homogenization Using a Reduced Basis and Effi-cient Sampling. Mathematical and Computational Applications 24 (2), 56.URL
Maraghechi, S., Rokoˇs, O., Hoefnagels, J. P. M., Peerlings, R. H. J., Geers, M. G. D., 2020.Harvesting micromorphic fields from experiments on patterning metamaterials. Submitted,1–23.Mark, A. G., Palagi, S., Qiu, T., Fischer, P., 2016. Auxetic metamaterial simplifies softrobot design. 2016 IEEE International Conference on Robotics and Automation (ICRA),4951–4956.Miehe, C., 2003. Computational micro-to-macro transitions for discretized micro-structuresof heterogeneous materials at finite strains based on the minimization of averaged incre-mental energy. Computer Methods in Applied Mechanics and Engineering 192 (5), 559 –591.URL
Miehe, C., Bayreuther, C. G., 2007. On multiscale FE analyses of heterogeneous structures:from homogenization to multigrid solvers. International Journal for Numerical Methodsin Engineering 71 (10), 1135–1180.URL https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.1972
Miehe, C., Koch, A., 2002. Computational micro-to-macro transitions of discretized mi-crostructures undergoing small strains. Archive of Applied Mechanics 72 (4-5), 300–317.URL https://doi.org/10.1007/s00419-002-0212-2
Mirzaali, M. J., Janbaz, S., Strano, M., Vergani, L., Zadpoor, A. A., 2018. Shape-matchingsoft mechanical metamaterials. Scientific Reports 8, 965.Nicolaou, Z. G., Motter, A. E., 2012. Mechanical metamaterials with negative compressibilitytransitions. Nature materials 11 (7), 608–613.Niknam, H., Akbarzadeh, A. H., 2018. In-plane and out-of-plane buckling of architectedcellular plates: Numerical and experimental study. Composite Structures 206, 739 – 749.URL
Ohno, N., Okumura, D., Noguchi, H., 2002a. Microscopic symmetric bifurcation conditionof cellular solids based on a homogenization theory of finite deformation. Journal of theMechanics and Physics of Solids 50 (5), 1125–1153.32hno, N., Okumura, D., Noguchi, H., 2002b. Microscopic symmetric bifurcation conditionof cellular solids based on a homogenization theory of finite deformation. Journal of theMechanics and Physics of Solids 50 (5), 1125–1153.Okada, J. I., Washio, T., Hisada, T., 2010. Study of efficient homogenization algorithms fornonlinear problems. Computational Mechanics 46 (2), 247–258.URL https://doi.org/10.1007/s00466-009-0432-1
Okumura, D., Ohno, N., Noguchi, H., 2002. Post-buckling analysis of elastic honeycombssubject to in-plane biaxial compression. International Journal of Solids and Structures39 (13-14), 3487–3503.Rokoˇs, O., Ameen, M. M., Peerlings, R. H. J., Geers, M. G. D., 2020a. Extended micro-morphic computational homogenization for mechanical metamaterials exhibiting multiplegeometric pattern transformations. Extreme Mechanics Letters 37, 100708.URL https://linkinghub.elsevier.com/retrieve/pii/S2352431620300699
Rokoˇs, O., Zeman, J., Doˇsk´aˇr, M., Krysl, P., 2020b. Reduced integration schemes inmicromorphic computational homogenization of elastomeric mechanical metamaterials.Advanced Modeling and Simulation in Engineering Sciences 7 (1), 19.URL https://amses-journal.springeropen.com/articles/10.1186/s40323-020-00152-7
Rokoˇs, O., Ameen, M. M., Peerlings, R. H. J., Geers, M. G. D., 2019. Micromorphic compu-tational homogenization for mechanical metamaterials with patterning fluctuation fields.Journal of the Mechanics and Physics of Solids 123, 119–137.Saiki, I., Terada, K., Ikeda, K., Hori, M., 2002. Appropriate number of unit cells in arepresentative volume element for micro-structural bifurcation encountered in a multi-scale modeling. Computer Methods in Applied Mechanics and Engineering 191 (23), 2561– 2585.URL
Sperling, S. O., Rokoˇs, O., Ameen, M. M., Peerlings, R. H. J., Kouznetsova, V. G., Geers,M. G. D., 2020. Comparison of enriched computational homogenization schemes appliedto pattern-transforming elastomeric mechanical metamaterials. In preparation.Wadee, M. A., Farsi, M., 2015. Imperfection sensitivity and geometric effects in stiffenedplates susceptible to cellular buckling. Structures 3, 172 – 186.URL
Yang, D., Mosadegh, B., Ainla, A., Lee, B., Khashai, F., Suo, Z., Bertoldi, K., Whitesides,G. M., 2015. Buckling of elastomeric beams enables actuation of soft machines. AdvancedMaterials 27 (41), 6323–6327.URL https://onlinelibrary.wiley.com/doi/abs/10.1002/adma.201503188 https://linkinghub.elsevier.com/retrieve/pii/S0020768320300123
Zheng, X., Lee, H., Weisgraber, T. H., Shusteff, M., DeOtte, J., Duoss, E. B., Kuntz, J. D.,Biener, M. M., Ge, Q., Jackson, J. A., Kucheyev, S. O., Fang, N. X., Spadaccini, C. M.,2014. Ultralight, ultrastiff mechanical metamaterials. Science 344 (6190), 1373–1377.URL https://science.sciencemag.org/content/344/6190/1373https://science.sciencemag.org/content/344/6190/1373