Towards shock absorbing hyperelastic metamaterial design (II): a prospective multiscale buckling-lattice computational model

Computational Mechanics https://doi.org/10.1007/s00466-024-02593-y ORIGINAL PAPER Towards shock absorbing hyperelastic metamaterial design (II): A prospective multiscale buckling-lattice computational model J. Cante1,2 ·A. Nuñez-Labielle1,2 ·A. E. Huespe4,5 ·J. Oliver1,3 Received: 24 October 2024 / Accepted: 19 December 2024 © The Author(s) 2025 Abstract As a continuation of a previous work of the authors, on Computational Design of Shock-absorbing Metamaterials (Part I) (NunezLabielle et al. in Comput Methods Appl Mech Eng 393:114732, 2022), this work explores the potential of computational multiscale methods, in combination with massive buckling-lattice structures at the metamaterial core (meso/micro scale), to render a suitable framework for designing such a shock-absorbing metamaterials focusing on industrial applications. In this context, a prospective computational setting is considered under the hypothesis that, for a sufficiently complex microlattice topology, some localized regions might buckle within the lattice-structure core and propagate through it, giving rise to different loading-unloading paths, in such a way that hysteretic-like structural behaviours would take place, thus arising dissipative behaviours, even if the base material at the buckling micro-lattice behaves in a hyperelastic (thus intrinsically nondissipative) manner. Using the standard Hill-Mandel homogenization principle, and assuming that the necessary separation of scales holds, the homogenized body, now living in a classical solid-mechanics setting, displays a homogenized non-convex behaviour which, in agreement with the conclusions of Part (I) of the work, exhibits extrinsic dissipation and, thus, could be potentially used (at reduced computational cost) for shock absorbing metamaterials analysis and design purposes. A tentative industrial application, to a sneaker’s insole design, has been then considered as a work’s target for evaluation of the room offered by the explored setting in the context of shock-absorbing metamaterial design. Finally, remarks on the scope and limitations of the work, and its significance for further advances in the field are emphasized. Keywords Shock absorbing metamaterials ·Computational metamaterial design ·Multiscale material modeling ·Buckling micro-lattice materials BJ. Oliver oli[email protected] J. Cante [email protected] A. Nuñez-Labielle [email protected] A. E. Huespe [email protected] 1Centre Internacional de Mètodes Numèrics en Enginyeria (CIMNE), C/ Gran Capitá, S/N, Edifici C1, 08034 Barcelona, Spain 2Escola Superior d’Enginyeries Industrial Aeroespacial i Audiovisuals de Terrassa (ESEIAAT), Universitat Politècnica de Catalunya - BarcelonaTech (UPC), Campus Terrassa, C/ Colom 11, 08222 Barcelona, Spain 3Escola Tècnica Superior d’Enginyers de Camins, Canals i Ports de Barcelona (ETSECCPB), Universitat Politècnica de 1 Motivation This paper corresponds to the second part of a research work devoted to exploring the possibilities and benefits of computational multiscale methods to be used in the analysis and design of shock-absorbing metamaterials for industrial applications. In Part (I) of the work [1] a first approach to the overall problem was made based on a single-scale modeling of a solid, endowed with a classic computational large-strain hyperelastic model, with devised perturbations in the voluCatalunya - BarcelonaTech (UPC), Campus Nord, C/ Jordi Girona 1, 08034 Barcelona, Spain 4CIMEC-UNL-CONICET, Predio Conicet, Ruta Nac. 168 s/n - Paraje El Pozo, 3000 Santa Fe, Argentina 5Programa de Engenharía Mecânica, COPPE, Universidade Federal do Rio de Janeiro, Cidade Universitária, Rio de Janeiro, RJ 21941-972, Brazil 123 Computational Mechanics metric part of the free energy such us to override the original convexity of the model (loss of convexity). There, it was displayed that this action: •Breaks the non-dissipative paradigm for poly-convex hyperelastic models, •Yields the formation of mechanical shock waves, propagating across the solid, and •Results in the production of extrinsic dissipation associated both to the propagation speed of the waves and the intensity of the local Eshelby stress tensor jump [2–4]. This fact potentially endows homogenized continuum hyperelastic models, with partially overridden convexity in the free energy, the ability to be used in shock-absorbing computational models. •In addition, those models exhibit the null strain recovery under null stresses property i.e.: the produced extrinsic dissipation, and the corresponding shock absorbing capability, is not rate dependent (as, for instance, in classic viscoelastic dissipators), the dissipation character of any shock absorbing device (shock absorber) built with that type of material, is potentially instantaneous and independent of the rate of the produced strains. In other words: after a first shock, the absorber gets ready to absorb a subsequent shock immediately, with a fully intact dissipation potential. This is of crucial relevance in industrial applications aiming at object protection from fast repeated impacts (i.e. package protection during transportation) as well as in other fields (sport-wearing material design etc.) intending to control the amount of dissipation and its evolution. However, it is well known that such hypothesized elastic dissipative materials do not appear naturally. Most of elastic materials1exhibit, at least at initial stages of deformation, a pure elastic strain-behaviour and they behave as non-dissipative. But, what if in the modern setting of metamaterials2one could talk of, and even manufacture, hyperelastic materials exhibiting some kind of dissipation?. The purpose of this second part of the work is to explore the domain of existing computational methods and settings for computational metamaterial design, and provide some new insights on the subject of computational design of shockabsorbing metamaterials. 1.1 The role of scales Metamaterials, artificially engineered materials, yet manufacturable, showcasing unusual (and extreme) physical prop1Historically found in nature or even manufactured along the last decades of the industrial era. 2Artificially engineered materials showcasing unusual (and extreme) physical properties. erties, are currently a focal point in material science. This is primarily because these peculiar properties, though displayed at the observable scale of the material (the macroscale), are achieved through its interaction with artificial structures in the lower scales (the meso/micro scales) in a way that defies conventional physical intuition. In spite that these structures may sometimes be very intricate, they can nowadays be manufactured using advanced techniques, like additive manufacturing. The selection of the intended upper-scale properties, the observed properties, and the corresponding low-scale structures as well as their cause, are the goals of the concept of Materials by Design [5]. Computational Mechanics methods offer increasing room for such a design, in combination with experimental, in-lab, methods which has coined the term Computational (Meta) material Design (CMD) [6–8]. In the context of shock-absorbing metamaterials design a commonly chosen low-scale feature, to be responsible for the shock-absorbing properties of the metamaterial, is the structural buckling [9–13]. 1.1.1 Buckling as a source of dissipation: the buckling microlattice Buckling (or structural instability in the most general meaning) is a well-known phenomenon in structural mechanics appearing in slender structures which, despite being constituted by a convex hyper-elastic material, exhibit, beyond certain load limits, an unstable structural behavior, normally associated with sudden and large geometrical changes in the structural members: i.e. structural instability, snapthrough/snap-back responses, which manifest as hysteretic action-response histories. This suggests that if a massive buckling lattice3,thebuckling microlattice constituting the core of a certain device, is observed as a solid mechanics body subjected to external loading actions, it might experience local buckling at certain regions of the microlattice causing sudden instabilities and large geometrical changes [14–18]. If the size of these buckling regions is sufficiently small concerning the typical size of the body, these buckling phenomena can propagate throughout the microlattice causing local buckling/unbuckling behaviors, which translates into hysteretic global mechanical responses [19–21]. In other words: from the Solid Mechanics point of view, at the observable (macroscale) some global mechanical dissipation would occur. In fact, this has been experimentally observed: in [22] it is reported that observation of in-lab mechanical loading/unloading experiments using laser microscopy (3D laser lithography), on experimental specimens endowed with a certain low-scale buckling lattice morphology (i.e. a certain 3In the remaining of this work the term massive lattice should be understood as a topologically dense lattice. 123 Computational Mechanics number of buckling cells with determined initial topology and shapes), display hysteretic action-response (forcedisplacement) in loading/unloading histories. However, upon complete removal of the actions the response became null, and the specimen morphology remains unchanged concerning the original one. In addition, repetition of the experiments displays very similar action-response paths. 1.1.2 Hierarchical multiscale modeling: linking the scales During the last decades, hierarchical multiscale modeling has gained increasing credit, in the Computational Mechanics community, as a powerful numerical tool to obtain new insights for modeling complex material behavior [23–25]. Hierarchical multiscale modeling is an idealization of a highly complex body (to be designed and manufactured at a low scale), which is replaced by a much simpler one: the homogenized4material representation at the high-scale fulfilling: •The Continuum Mechanics conservation and energy balance laws •A scale-bridging principle: typically an energetic equivalence principle linking the idealized (homogenized) material representation and the original material resulting in a homogenized continuum body model, retaining the essential behavior of the original material and providing useful insights for Computational Material Design at affordable computational costs. The homogenization procedure can be, then, summarized as follows: 1. The highly complex, mechanical behavior of a deformable body at the observable scale, is considered the result of the hierarchical interaction of simpler mechanical behaviors linked to a nested sequence of representative material scales, assumed separable in terms of their corresponding, decreasing, representative sizes. 2. The material behavior at every scale is governed by specific mechanical laws, in terms of the values of some, conveniently identified, state and internal variables, which are aimed at being obtained via the computational model equations. 3. A hierarchical-link principle is then postulated for every two consecutive scales (from top to bottom). That link involves some energy-equivalence laws, expressed as 4The concept of homogenization refers here to the determination of the stress-strain measures evolution, at one point of the upper scale, in terms of their energetic equivalence with a portion of the low scale (the Representative Volume Element, RVE) displaying a certain actionresponse according to the so-called Hill-Mandel Principle. Details are given in section 1.2. variational equations, for the evolution of the state variables at every material point at the upper scale and the corresponding average density evolution, at the lower scale, in a conveniently shaped and sized domain termed the Representative Volume Element (RVE). The adequacy and physical significance of that link, in terms of the representation of the material behavior, is, obviously, crucial for the accuracy of the results, which is assumed to improve asymptotically with the so-called scale separation, understood as the ratio of each scale size and the immediate lower-scale RVE size (the larger the better). The solution of the resulting variational equations provides •Determination of the upper-scale state variables in terms of the lower-scale ones (the homogenization procedure) and •A nested variational problem at the lower scale (the RVE problem), to be solved, for every set of kinematic entries downloaded from the upper scale, in terms of some low-scale kinematic variables the fluctuating displacements,livingattheRVE. Recursive (bottom to top) resolution of the problem through scales, allows stating the evolution problem at all scales, now including all the physical interactions among material scales, which, once solved, provide the (top to bottom) state variables evolution. The effectiveness of this idealized multiscale material representation depends on the appropriate value of the so-called scales-separation (the larger the better), which might be assessed by asymptotic comparisons with simpler Direct Numerical Simulations (DNS) at the lower scale. 1.1.3 Hierarchical multiscale modeling for a buckling-microlattice shock absorber In light of the above considerations, an approach is explored in this paper for the hierarchical multiscale computational material design of a buckling microlattice to be used in 2D representations of lightweight shock-absorbing devices. It is considered that: 1. At the buckling microlattice, the material points exhibit alinear material behavior (elastic strains and convexnon-dissipative constitutive behavior) in a geometrically non-linear environment (characterized by large displacements and small strains) which, for a sufficiently complex microlattice morphology and topology may arise local buckling-regions in the lattice. In these conditions the following assumptions are made: 123 Computational Mechanics •Mechanical instability (microlattice buckling) can occur as dynamic effects, producing buckling phenomena in regions of the lattice, which can buckle, unbuckle, or exhibit a stable behavior over time. In turn, these buckling phenomena can propagate through different regions of the lattice. •This buckling propagation may be irreversible, i.e.: for certain structural loading/unloading paths the buckling/unbuckling affected microlattice regions are not the same, thus resulting in a hysteretic actionresponse history, i.e.: the global action-response history is different for loading and unloading processes. 2. Regarding the overall shock-absorbing specimen (the shock absorber) it is assumed that: •The microlattice regions behave as homogenized particles of a macroscale continuum, thus ruled by the standard conservation laws of solid mechanics and the homogenized constitutive law, in terms of the homogenized measures (homogenized stresses and strains). •According to these conservation laws the shock absorber, endowed with the homogenized material, exhibits a global dissipative behavior, emerging from the hysteretic action-response (structural response) of the microlattice, and ruled by a non-convex homogenized constitutive model. In other words: homogenization renders the structural (convex elastic) behaviour of the particles at the microlattice members into a non-convex (solid mechanics-like) constitutive model at the corresponding points in the homogenized macroscale continuum. This translates into an overall dissipation in the (macroscale) shock-absorbing device problem. 1.2 Proposed two-scale approach for lightweight buckling microlattice shock absorbers In this work, the general multiscale setting described above is specified as a two-scale approach as follows: 1. Only two scales are considered in the shock-absorbing metamaterial i.e.: •The (observable) macroscale: it constitutes a deformable solid (shock absorber device) ruled by the classical large strain kinematics, in terms of (homogenized) strains and stresses (here termed the homogenized measures), and fulfilling the standard conservation laws in solid mechanics problems (i.e. linear and angular momentum conservation). Their approximate solution, via a finite element discretization problem at the macroscale (the macroscale problem), provides a prediction of the displacement field, and the homogenized strains, and stresses in the global system (impactor-absorber). This solution also yields the evaluation of the accumulated system dissipation, in the impactor-absorber system. This dissipation is anticipated to be a scalar measure, computed as the difference between the external energy supplied to the shock absorber and the sum of its whole amount of internal energy (free energy plus kinetic energy), and it should be (for thermodynamic reasons) always time-increasing and positive. It is also anticipated that in this work, determination, optimization, and control of the system dissipation are considered the fundamental issues in the design of the shock-absorbing metamaterial. •The lower-scale: it is constituted by the material particles at a massive buckling microlattice in the shock absorber, which is characterized by the morphology and topology of the aforementioned chosen region of the lattice, the RVE, the whole massive buckling lattice being considered amenable to be generated through suitable repetitions of the RVE along the lowscale dimensions. 2. The low-scale kinematics (displacement and strains definition) is obtained as a (linear) Taylor’s expansion, along the RVE, of the macroscale displacements and the corresponding homogenized strain fields, supplemented by an unknown fluctuating displacements field,livingatthe RVE and playing the role of a correction, at the low-scale, of that linearized displacement field thus providing additional accuracy to the homogenized solution concerning the, theoretically exact, DNS (Direct Numerical Simulation) solution. 3. Imposition of the mechanical energy conservation paradigm (Hill-Mandel Principle) across the scales, which is numerically solved via a variational statement. The resulting solution provides [24]: •The homogenized (RVE-averaged) stress and strain fields in the macroscale •The RVE problem: a variational problem solving the fluctuating displacements, in a finite element discretized problem on the RVE domain. These additional (fluctuating displacements) unknowns are assumed to correct the macroscopic problem solution in terms of the homogenized solutions, yielding a higher accuracy, towards the (exact) DNS solutions, as the fluctuations decrease to zero (i.e. as the scale-separation increases). In a way, they compensate for the error of computing a problem solution, in the reduced-size RVE of the microlattice, instead of solving the (non-affordable in terms of computational cost) DNS solution. 123 Computational Mechanics 4. The RVE problem is additionally subjected to some restrictions on the fluctuating displacements: the consistency conditions. They enforce the nullification, at the RVE domain, of both the averages of the fluctuating displacements and their spatial gradient values. In a way, they guarantee the convergence of the multiscale solutions in terms of the scale separation: the larger the scale separation the more negligible the fluctuating displacements and strains, and the higher the accuracy of the obtained homogenized quantities, concerning the DNS solution. 5. In the hierarchical multiscale approach context, the representative character of the chosen RVE is assumed5.In this work, this is deemed achieved by using, as RVE, a suitable portion of the actual buckling microlattice fulfilling the following conditions: •The actual microlattice topology and morphology is obtainable by spatial repetition of the RVE, endowed with periodic boundary conditions, along all the dimensions of the problem. •At the same time, and in combination with the energy equivalence principle, the complexity of the RVE is enough for efficiently capturing the physical phenomena intervening in the propagation of buckling across the microlattice and, thus, to reproduce the corresponding hysteretic behavior at the macroscale causing the dissipation. In addition, the RVE size and complexity keep in balance the accuracy of the obtained simulations and their computational cost6. 6. Due to the lightweight character of the considered buckling microlattice, it is stated that the homogenized (macroscale) inertial forces in the shock absorber are negligible in comparison with the ones that arise in the impacting bodies. Therefore, the solid mechanics problem in the shock-absorbing metamaterial is considered quasistatic. Such an approach is adopted in this work aiming at rendering a suitable setting to face the computational design of shock-absorbing metamaterials, balancing the accuracy of the model with a moderate computational cost when compared with the more realistic, but unaffordable high-cost, DNS models. 5This is to say: the RVE is sufficiently small (and topologically detailed), with respect to the overall shock absorbing specimen, to yield structural responses which are statistically representative of the considered mechanical deformation process. 6When compared with DNS simulations. 2 Buckling microlattice modeling. Dimensional reduction of the RVE to a set of buckling beams Based on the motivational arguments in Sect. 1,sofarthelow scale has been considered a buckling microlattice, made of slender components, which, consistently with the 2D character of the assumed Solid Mechanics problem, constitute geometrical 2D objects (see Fig. 1b). Therefore, in principle, it should be consistently modeled as a 2D continuum. However, for the prospective purposes of this work, the slender character of the buckling components is hypothesized to be such that they could be replaced, with little error, by 1D buckling structural elements (buckling-beams), which could be modeled in the context of elastic large-strain beam theory. The corresponding Representative Volume Element (RVE) would then become the set of structural buckling beams in Fig. 1c. This dimensional reduction of the RVE model, considered so forth, will translate into a very relevant reduction in the computational cost, but still provides substantial accuracy of the results7. In consequence, the plane beam model sketched in Fig. 2, taken from Felippa [26], will be used to capture the structural behavior of the RVE domain at the lower scale. For reasons of completeness, a brief discussion about this model will be presented in this section. The main features of the model are then: (1) the reduced (1D), dimensional character, with respect to the original one (2D), which is assumed to yield enough accuracy, (2) large displacements and rotations, but small strains, kinematics in a Total Lagrangian formulation, and (3) a convex linear elastic model (Saint-Venant’s model) considered for the beam constitutive material. 2.1 Kinematics Let’s consider a beam of length L, with cross-section A, symmetric in the loading plane, having a straight neutral axis in the reference configuration passing through the centroid of the cross-section and identified by Y0(see Fig. 2). The unit vector ˆ Nis parallel to the beam neutral axis. Considering a local coordinate system (s,η) associated to the orthogonal basis {ˆ N,ˆ T}, where ˆ Tis the orthogonal unit vector to ˆ N,see Fig. 2, the points of the beam neutral axis in the reference configuration, denoted Y0, might be written in terms of the local coordinate (arc-length parameter) s:Y0(s)=Y1+ 7Notice that modeling the large deflections of the buckling components of the RVE, assumed endowed of very small thickness, would imply very dense discretizations of the corresponding (low-aspect-ratio) 2D finite elements meshes to avoid numerical ill-conditioning. 123 Computational Mechanics Fig. 1 The two scales problem. aThe macroscale, consisting of a deformable body governed by the classical large strain kinematics, described in terms of homogenized strains and stresses, bRVE, composed of slender structural 2D structural elements, cthe simplified RVE composed of 1D structural buckling elements (beams) Fig. 2 Felippa’s beam theory. (a) Kinematic description; (b) Considered 1D linear finite element s(Y2−Y1)/L=Y1+sˆ N. Considering now a material point Yoccupying the material configuration at: Y=Y0+ηˆ T(1) and assuming a deformed beam displaying small axial stretches and shear distortions, the deformation map of this material point might be expressed as: y(Y,t)=Y0+u0(Y0,t)+ηK(Y0,t); k=−sin(θ) ˆ N+cos(θ) ˆ T;(2) where the displacement vector of points on the beam neutral axis is denoted u0, the unit vector kis parallel to the beam cross-section in the current configuration, while in the reference configuration is orthogonal to the neutral axis and parallel to ˆ T. Thus, θis the cross-section orientation angle in the current configuration with respect to its direction in the reference configuration, γ=θ−ψis the shear distortion angle, and ψis the angle of the beam neutral axis direction in the actual configuration and the direction of the same axis in the reference configuration. Then, the displacement vector at Yis defined as: u(Y,t)=y−Y=u0(Y0,t)+η(k−ˆ T)(3) and the deformation gradient F=∇ Yyis evaluated in consequence, resulting: F=(ˆ N+u0 ,s+ηk,s)⊗ˆ N+k⊗ˆ T(4) where notation (·),srefers to the derivative of entity (·)with respect to s. Also, the displacement gradient Jis computed as J=F−1, where 1represents the second-order identity tensor. Note that in this approach, and from expression (2), the kinematic descriptors of the model presented in [26]areu0 and θ. Furthermore, based on the definition of the distortion angle γ, the curvature term is, κ=θ,s, assuming small strains and small γ, and disregarding higher order terms. Felippa in [26] derives a consistent linearization, which is here adopted. Taking into account all of these ingredients, the resulting Green-Lagrange deformation tensor E, after linearizing the conventional expression E=1 2(FTF−1), can be written 123 Computational Mechanics as follows: E=e−ηκ γ/2 γ/20 (5) where the axial, e, and shear, γ, strain measures are given by: e=1+d(u·ˆ N) dS cos(θ) +d(u·ˆ T) dS sin(θ) −1 γ=−1+d(u·ˆ N) dS sin(θ) +d(u·ˆ T) dS cos(θ) (6) respectively. Additional details about the linearization used to derive the Green-Lagrange strain tensor (5) from the deformation gradient (4) can be seen in [26]. The triplet (e,γ,κ) constitutes the generalized strain measures of the beam model. Then, the generalized strain vector will be denoted as h=⎡ ⎣ e γ κ⎤ ⎦(7) 2.2 Constitutive model and stresses By assuming a linear elastic material, the first PiolaKirchhoff stress Pin the proposed beam theory is given by: P=PnN(ˆ n⊗ˆ N)+PnT (ˆ n⊗ˆ T) +PtN(ˆ t⊗ˆ N)+PtT(ˆ t⊗ˆ T)(8) Equivalently the Second Piola-Kirchhoff stress tensor S is: S=SNN(ˆ N⊗ˆ N)+SNT(ˆ N⊗sym ˆ T)(9) where the notation ⊗sym indicates the symmetric tensor (open) product of the two vectors. Consistently with the beam theory developed in [26], the component STT is equal to 0. This non-null component can be expressed in terms of the generalized strains using the following constitutive equation: SNN SNT =E(e−ηκ) Gγ(10) where the symbols Eand Grepresent the Young’s and shear moduli of the beam material, respectively. Finally, the stress tensor S(Y)can be integrated through the beam crosssection, providing the generalized stresses: N=A SNNdA=EAe V=A SNTdA=GAγ M=A −ηSNNdA=EIκ (11) where Nis the axial force, Vis the shear force, Mis the bending moment, Ais the area of the cross section and Iis the second moment of inertia. 2.3 Homogenization of the buckling micro-lattice. Multiscale formulation As described in Sect. 1a multiscale problem displaying two characteristic well-separated length scales is analyzed. At the macro or coarse scale, a continuum model with an exact non-linear kinematic description is considered. At the micro or fine scale, a beam network is modeled using the approach presented in the previous section. An idealized sketch of this problem is depicted in Figure 1. For simplicity, it is also assumed that, at the microscale, the sections of each beam are uniform and that there are no discrete forces, discrete moments, or body forces applied to the microstructure. Therefore, the internal generalized forces remain uniform along each microstructural frame element. Thus, the mechanical variables at the macroscale are the position of the reference point, X, the spatial position of the point, x(X), the displacement vector, uM=x−X, the gradient of deformation, FM, and the Piola stress, PM. The microscale mechanical state associated with the macroscale point Xis defined by the position of the reference point Yand the spatial position of the point y(Y). Additionally, following the beam theory of Sect. 2, the displacement vector of the beam neutral fibers is u0 μ=y(Y0)−Y0,the gradient of deformation is Fμ, and the Piola stress is Pμ. 2.3.1 Scale bridging equations The displacement vector uμof any point Yat the microscale, in agreement with (2), can be written as follows: uμ=uM(X)+JM·(Y−YC)+˜ u0 μ(Y0)+ημ(˜ kμ−ˆ T) (12) where, after imposing the condition that YCis the mass center at the RVE, results: μ JM·(Y0−YC)dμ=0(13) 123 Computational Mechanics Under the hypothesis of periodic boundary conditions, the following mathematical expressions are fulfilled (see Fig. 19): u0 μY+=u0 μY−+JM·Y+−Y−, ˜ u0 μY+=˜ u0 μY−, ˜ θμY+=˜ θμY− (14) 2.3.2 Hill-Mandel principle The Multiscale Virtual Power Principle proposed in [24]is adopted in this work. It corresponds to a variational statement of the classical Hill-Mandel principle [27,28], in which an energetic equivalency between the two scales is postulated. It establishes that the stress power at a given particle at the macroscale, X, equals the mean value of the stress power at the corresponding RVE domain, μ(X), i.e.: PM(X):δFM(δuM)=1 |μ|μ Pμ(Y):(δFM(δuM) +∇Yδ˜ uμ)dμ(15) for any admissible variations, δFM, belonging to the full set of second order tensors and for any virtual kinematically admissible fluctuation,δ˜ uμ(see Appendix Afor detailed information on the notation and definition of spaces). Equation (15) defines a fundamental link between the two scales. Direct consequences of it are: •The homogenized stress expression, PM=1 |μ|μ Pμuμ(Y)dμ(16) •The RVE problem at the microscale, i.e.: Find uμfulfilling: μ Pμ(uμ):∇ Yδ˜ uμdμ=0, ∀δ˜ uμ,(admissible f luctuation)(17) As a consequence of Eq. (17), the homogenized stress (16) can be expressed in terms of the traction vector tμ=Pμ·Nμ evaluated exclusively at the boundary of the RVE, i.e., PM=1 μμ tμ⊗Ydμ(18) where μ=∂ωs μ∩∂μ, (see Fig. 1b). 2.3.3 From continuum to buckling micro-lattice Let’s assume that the arrangement of slender deformable solids can be represented by a straight beam network, μ (Fig. 1c), modeled using the approach presented in Sect. 2.1. By following the mathematical manipulations detailed in Appendix A, the direct consequences of the Hill-Mandel principle can be equivalently expressed in terms of the beam terminology as follows: •The RVE problem at the microscale Find (u0 μ,θ μ)fulfilling μ z(u0 μ,θ μ)·δ˜ hμδ˜ u0 μ,δ˜ θμdμ=0, ∀δ˜ u0 μ∈˜ V0 μ,∀δ˜ θμ∈˜ μ(19) where zμ=(Nμ,Vμ,Mμ)Tis the generalized stress vector, δ˜ hμ=(δ ˜eμ,δ˜γμ,δ˜κμ)Tis the variation of the generalized fluctuating strain vector, ˜ V0 μ=δ˜ u0 μ:μ→R2|δ˜ u0+ μ=δ˜ u0− μ,(20) and ˜ μ=δ˜ θμ:μ→R|δ˜ θ+ μ=δ˜ θ− μ.(21) •Homogenized stress expression computed in terms of the reaction forces at the boundary of the RVE: PM=1 μ nb  i Rμi⊗Y0i(22) where Rμi=μi Pμ·Nμdμ(23) represents the reaction force vector at node i:1...nb.8 3 Finite element approach Let h μ= nbeams  α=1 α μ(24) denote the finite element discretization of the polygonal curve that represents the midline of all the beams comprising the 8nbrepresents the number of boundary nodes in the RVE. 123 Computational Mechanics RVE. Based on [26], a C0continuous interpolation with two nodes under each element, α μ, can be used to interpolate the micro-displacements (of the neutral axis) and rotations as linear functions of the nodal parameters u0 μ(S,t) θμ(S,t)α =Nα 1(S)u0 μ(t) θμ(t)1 +Nα 2(S)u0 μ(t) θμ(t)2 (25) where Nα 1and Nα 2are the linear shape functions while u0 μ(t)=[ u0 μ,Y1u0 μ,Y2]Tand θμ(t)9are the interpolation parameters associated to nodes 1 and 2, respectively. According to the strains definition (6) and the interpolation expression (25), the discrete form of the variation of the generalized strain vector (7), can be expressed as: δ˜ hα μ=B B Bα μδ˜ dα μ(26) where B B Bα μ, the strain-displacement matrix, contains the variations of the strains δ˜eμ,δ˜γμand δ˜κμwith respect to generalized displacement fluctuations δ˜ dα μ=δ˜ u1 μ,δ˜ θ1 μ,δ˜ u2 μ,δ˜ θ2 μ. For the definitions of the components of B B Bα μand a detailed explanation of their derivation, please refer to [26]. By applying the finite element discretization (24), (25) and the element Eq. (26), the variational problem (19) can be reformulated in its discrete form as follows: •Discrete variational problem: Find dμ,periodic, fulfilling δ˜ dμ Tfμdμ=0∀δ˜ dμperiodic (27) where fμ(dμ)= αh α B B BμαT(dα μ)zμα(dα μ)dα μ.(28) Here, zα μrepresents the stress-resultant vector at element level. 3.1 Exact linearization of the RVE problem As a result of the geometrically non-linear features inherent to the beam theory [26] used to capture the mechanical behavior of the RVE the problem in Eq. (27) becomes non-linear. Consequently, an iterative approach is crucial for effectively solving the system. To this end, the Newton-Raphson (NR) method is proposed, involving the exact linearization of the problem at each iteration as follows: δ˜ dμ TKTd(k) μ,n+1d(k+1) μ,n+1+fμd(k) μ,n+1=0 (29) 9Here, tis considered a pseudo-time, t∈[0,T], governing the evolution of the deformation process. where d(k+1) μ,n+1represents the increment of the generalized displacement vector at load step n+1 and during the internal Newton-Raphson iteration k. The tangent stiffness matrix, denoted by KTis expressed as: KTd(k) μ,n+1=Kmat d(k) μ,n+1+Kgeom d(k) μ,n+1,(30) where Kmat is the material stiffness matrix and Kgeom is the geometric stiffness matrix. For detailed computation procedures, refer to [26]). 3.2 Imposition of periodicity conditions Considering JM n+1=FM n+1−1the displacement gradient tensor at load step n+1. Let us express the fluctuation displacement vector in a more convenient form: ˜ dμ=˜ dint μ˜ d− μ˜ d+ μT (31) where ˜ dint μis the fluctuations displacement vector of the nodes in the interior of the RVE, while ˜ d− μand ˜ d+ μare the fluctuations of the nodes in − μand + μ, respectively (see Fig. 19). In this context, the periodicity conditions are expressed as ˜ d− μ=˜ d+ μ,(32) or, in component form, ˜ u0 μ ˜ θμi− =˜ u0 μ ˜ θμi+ (33) for i+:1...mand i−:1...m.10 This condition can be equivalently rewritten in terms of the generalized micro displacements: u0 μ θμi+ =u0 μ θμi− +JMT n+1Y0+−Y0− 0(34) or equivalently as d+ μ=d− μ+¯ dM(35) where ¯ dMi−=JMT n+1Y0+−Y0− 0(36) 10 mis the number of nodes in + μ, or equivalently, in − μ. 123 Computational Mechanics Fig. 11 Structural responses and evolution of the total extrinsic dissipation energy for the vertical loading and unloading process of a variable cross-section specimen, for each of the RVEs complexities: aone unit cell; btwo unit cells; c three unit cells; dfour unit cells Fig. 12 Comparison of structural responses and accumulated dissipation evolution in the variable cross-section specimen for different RVEs complexities 123 Computational Mechanics Fig. 13 Sneaker-insole as a shock absorbing metamaterial: a sneaker insole, bsneaker impact on the ground, ceffective runner’s mass, W, impacting on a specific cleat (taken from [34, 35]) 5.1 Micro-lattice-buckling as a source of dissipation in shock absorbing metamaterials The considered key aspects to be analyzed in this prospective exercise are: •The ability of the shock absorbing metamaterial to produce extrinsic dissipation and its translation into kinetic energy attenuation, •The amount and time-evolution of the dissipation along astep 27. This dissipation will translate into attenuation of the kinetic energy introduced by the runner’s effective weight on the cleat (represented by the upper impacting mass Win Fig. 13c), along a typical runner’s step, •The quick shape recovery property of the cleat after the impact (in order to get ready to absorb an impact very son after another one). •The amount of energy dissipation and its evolution over one runner’s step are considered the optimality criteria28. Figure 14a sketches the impact process of the effective mass,W, of the runner impacting on a single representative cleat, which is to be analyzed in terms of two possible (cylindrical/trunk-conical) macroscopic shapes. Figure 14b, instead, focuses on the low-scale, displaying a generic RVE consisting of four different cells/buckling-layers29. Every cell of the RVE is assumed constituted by (see Fig. 14b): 27 Assumed normalized to one-second duration. 28 To be determined in the design process through appropriated biomechanical criteria. 29 Every color in the figure represents a different cell. •Inclined buckling beams, with low stiffness, which are in charge of inducing vertical-buckling effects in the cell and •Vertical and horizontal stiffening bars, with much larger stiffness, which are in charge of precluding horizontal buckling in the cell. The elastic material properties of beams and bars are the same for all cells in the RVE, with the following values •Inclined buckling beams. Young’s modulus: Eb=32.5 GPa. •Vertical and horizontal stiffening beams. Young’s modulus: Es=325 GPa. The structural properties of the buckling bars change from layer to layer, since the thickness of the inclined beams is set to increase by a 12% from one layer to the next one (from top to bottom). This is done to induce a sequential buckling (also from top to bottom) of the cells in the RVE. The resulting RVE-homogenized constitutive model is obtained following the procedure in Sect. 4.2.1, and it is shown in Fig. 14c. Observation of Fig. 14 yields the following comments: •The homogenized constitutive model exhibits a nonconvex shape similar to the one in Fig. 5, thus displaying several bumps, each of them coinciding with the intended sequential buckling of the RVE layers in Fig. 6. •In turn, this type of, non-convex, constitutive response matches the ones discussed in Part (I) of this work [1], concerning single-scale models endowed with artificially imposed non-convex free energy, where differentiation 123 Computational Mechanics Fig. 14 Shock-absorber elements for design. aconstant cross-section I, vs. variable cross-section II for a cleat impacted by an effective weight Wproducing an initial kinetic energy of 386 N·mm,bRVE consisting of four different buckling layers (unit cells), color-coded for easy identification, cHomogenized constitutive response for the proposed RVE of the free energy yields similar bumping constitutive models, thus inducing propagating mechanical shocks and the corresponding extrinsic dissipation. •Therefore, in a multiscale micro-buckling lattice metamaterials context, this allows identifying the propagating buckling phenomena at the microscale as the source of the extrinsic mechanical dissipation. 5.2 Room for design at the macroscale (cleat-shape design) Let us now explore the room for design provided by the present prospective approach to the problem considered in Fig. 13a in terms of the cleat shape. Figure 15a, displays the dynamic problem of the effective runner’s mass, W, impacting the cleat for both configurations. In Fig. 15b, c, d an oscillatory evolution of the vertical displacement, δat the top of the cleat, along with the kinetic energy of the impacting mass, and the accumulated dissipation over time, are, respectively, depicted for the two cases under analysis (case I, displayed with green dashed lines, and case II, with blue solid lines). The results indicate that the shape of the shock absorber significantly influences the physical response. Its observation indicates that: •The dissipation evolves much faster in case II than in case I. The amount of accumulated dissipation is similar in both cases (though a little higher for case II. •At the end of the dissipation process, some (small) oscillations persist, their magnitude being significantly smaller in case II than in case I. This could have been expected, since the hysteretic behavior in the action/response curves (see for instance Fig. 12a) ends before the static equilibrium point, (F,δ) ≡(0,0), and the model exhibits a (very small) non-dissipative oscillatory response around this point. •Finally, Fig. 15c confirms that the shape of the cleat significantly influences both the evolution of the accumulated dissipation and the total dissipated energy. At the same mass, a trunk-conical shape of the cleat produces faster and larger dissipation than a cylindrical one. In summary: the explored framework provides a substantial room for design in terms of the cleat’s shape. The cleat’s shape allows for control of both the evolution of the extrinsic dissipation and the total final amount of dissipation. 5.3 Room for design at the microscale buckling lattice (RVE morphology design) Figure 16 displays different RVE-morphologies to be considered, each one containing a different number of bucklinglayers/cells: 2, 4, and 8, respectively (to enhance visualization, every cell is shaded in a distinct color). The thickness of the buckling beams at each cell/buckling-layer is then adjusted to ensure that the average density of the three RVEs preserves the same total mass for all RVEs (same mass for the overall shock-absorber), and that, under the action of the compressive load, the different cells of the RVE buckle in a sequential (top-to-bottom) manner (Fig. 17). Figure 18 presents the results for the evolution of the displacements on the top face of the shock absorber, the kinetic energy, and the accumulated dissipation for the three proposed microscale configurations. The results demonstrate a 123 Computational Mechanics Fig. 15 Design variables in the process of designing a single cleat acylindrical-shaped cleat I, vs., trunk-conical shaped cleat II; bdisplacement attenuation curve; caccumulated dissipation; dhistory of the kinetic energy clear influence of the microscale morphology on the dissipation properties featured by the shock absorber. Figure 18c confirms that the dissipation performance, assessed through both the total accumulated dissipation and the time taken to achieve this value, strongly affected the RVE morphology. In these terms, case III appears the most efficient, while case I appears the least efficient one. This conclusions are supported by the evolution of displacements (Fig. 18a) and the corresponding changes in kinetic energy (18b). In case III, kinetic energy reaches its minimum rapidly, and the amplitude of subsequent oscillations is minimized. Finally, when considering the number of levels in each RVE as a reference, the proportion by which these levels increase in each configuration does not correspond to the increase in dissipated energy. This emphasizes the high degree of non-linearity and complexity of the phenomenon. 6 Final considerations and significance for future developments 6.1 Considerations on the explored computational setting for shock-absorbing metamaterials design In previous sections, a tentative computational setting for modeling shock-absorbing metamaterials has been explored. It relies on two main ingredients: 1. Consideration of shock-absorbing devices constituted by a massive buckling-microlattice, in which buckling can locally arise and propagate (Sect. 1.1.1). 2. The use of hierarchical multiscale modeling and the HillMandel principle for building an energetically equivalent homogenized solid to provide meaningful results, at moderate computational cost (Sect. 1.2). 123 Computational Mechanics Fig. 16 Sensitivitytothe microscale morphology. Three different RVE morphologies. for a constant cross-section cleat. I: 2 cells RVE, II: 4 cells RVE, and III: 8 cells RVE (every cell is depicted in a different color). The thickness of the buckling beams at each cell is adjusted to ensure that: (1) the average density of the three RVEs preserves the same total mass, and (2) under the action of the compressive vertical load, the RVE cells buckle in a sequential (from top to bottom) manner Fig. 17 Sensitivitytothe microscale morphology: a macroscale impact problem, b resulting homogenized constitutive models corresponding to the three RVEs in Fig. 16 In such a framework a shock absorbing metamaterial, constituted by a massive buckling microlattice, has been devised30 (Sect. 2), and virtually tested in a potential industrial application: a sneaker insole design. Then, the numerical simulations display encouraging results for the purposes of computational metamaterial design, i.e.: a) Significant mechanical (extrinsic) dissipation takes place in loading-unloading impact cycles (Sect. 4). b) Dissipation is increasingly sustained under quickly repeated impacts (Sect. 5)31. 30 Assuming that sufficient scale separation exists to provide reliable results. 31 A specific beneficial feature for some applications (e.g. packaging). c) There is clear room for design at both scales, i.e.: the resulting mechanical dissipation32 strongly depends on the high-scale shape of the impacted object but, also, of the considered buckling lattice morphology at the low scale. This allows the insertion of the explored framework into wide-purpose metamaterial design settings. 6.2 Remarks and significance for subsequent developments on the topic The authors are aware of the prospective character of this work on this specific branch of metamaterials design (multiscale-based computational shock-absorbing models). This is why, in the following paragraphs, they would like to 32 Considered as the goal function of an assumed optimization procedure 123 Computational Mechanics Fig. 18 Effects of microscale morphology on the shock-absorber dissipation. aHomogenized constitutive responses for the three analyzed cases, bdynamic histories of the vertical (downwards) displacement of the top face of the shock-absorber, caccumulated extrinsic dissipation histories devolution curves of the kinetic energy of the impactor transmit to the reader some additional remarks arisen during the research work and writing of this paper i.e.: •The statement, “the extrinsic dissipation arises from buckling propagation at the microscale”, made in this work (see Sect. 1.1.3), has gained much credibility at the authors’ eyes along the writing. However, they have not been able to provide graphical evidence of it. From one side, consequences of that propagation across the microlattice (production of extrinsic dissipation) have been widely displayed along the present manuscript by means of the representative examples. However, graphical representation of such propagation could not bee presented here, since it would require a highly complex, and computational costly, calculation process of the propagating buckling-bands at the microscale33. This is where Part (I) of this work [1] comes into play through the following reasoning: 1. Comparison of the presented examples in this Part (II) with the ones in Part (I), both producing extrinsic 33 May be after being de-codified,intopropagating shock-wave lines at the macroscale dissipation, suggests that the ultimate reason for that dissipation should be the same in both cases. 2. In Part (I) (single-scale model) the dissipation was identified as a direct consequence of the perturbation of the original (fully-convex) hyperelastic constitutive material34. Then, the propagation of the resulting shock-waves was naturally displayed and explained in the context of the classical single-scale shock waves propagation theory [2,4,36]. 3. In this Part (II) (multiscale model), it is clearly displayed the stress homogenization procedure translates into non fully-convex homogenized constitutive models as well (see, for instance, Fig. 14c). 4. Therefore, the manifestation of such lack of convexity should be the same, in both cases i.e.: shock-waves propagation. In other words, the buckling propagation appearing in the ideally homogenized solid in this Part (II) and, therefore, in the tackled buckling microlattice shock absorber device, could be explained, through a suitable de-codification method of, on a similar basis than in Part (I), in the context of the classical single-scale shock waves propagation theory, 34 Which overrides its original fully-convex hyperelastic character. 123 Computational Mechanics and not necessarily by resorting to hidden sources of intrinsic local dissipation35. 5. In contrast with the extrinsic dissipation (associated with the buckling propagation mechanisms accounted for here), intrinsic36 dissipation sources are not considered relevant in the presented idealized scenario. •In this context, issues such as the uniqueness and stability of the mechanical system might be more properly examined in the propagating shock waves setting instead of in the, more standard, equilibrium/stability settings. •From the Solid Mechanics point of view, the way in which the convex hyperelastic material (thus non-dissipative in terms of classical standard materials) endows the corresponding shock-absorbing device with energy dissipation capabilities, can be attributed to multiscale coupling. In the proposed computational setting, this coupling is taken into account by the adopted hierarchical multiscale framework37. The homogenization procedure results in the homogenized stress-strain measures in equations which are implicitly related to each other through the resulting homogenized constitutive model which turns out to be non-convex (see for instance Fig. 5). •A proper evaluation of the extrinsic dissipation can be expected as far as a sufficient scale separation exists, between the considered shock absorber and the buckling lattice sizes. In simpler words: as far as the typical micro buckling lattice size38 is sufficiently small concerning the device dimensions. Appendix: A Computational multiscale modeling: from continuum to lattice micro-structures As introduced in Sect. 1.2, a two-scale approach for lightweight buckling-microlattice shock absorbers is proposed. At the macroscale, a continuum model with an exact non-linear kinematic description is employed. The mechanical behavior of this model is captured using traditional 2D or 3D solid modeling techniques39. In contrast, the microscale 35 For instance, high-frequency micro-oscillations in the buckling lattice, as some times has been proposed in the literature on the subject (e.g., [37]). 36 Inelastic-strains-based dissipation or thermal-like dissipation i.e.: amenable to be quantified in terms of the buckling-lattice particles massdensity. 37 And the restriction of scale-separation associated to it. 38 A measure of the buckling-cell size 39 For the sake of simplicity in this work the problem will be restricted to the 2D case is constituted by a massive buckling microlattice composed of slender deformable solids arranged in a lattice structure devised to develop a propagating buckling behavior within the microstructure. For computational saving reasons, such micro-structures will be modeled using degenerated beam/frame 1D kinematics in the RVE (see Fig.1). The use of the simplified 1D degenerated kinematics translates into large savings in the computational cost for the RVE solution, when compared with the actual 2D kinematics, but introduces an incongruity when developing the classical multiscale theory40 which is here solved by resorting to the limit case of an Ersatz material containing slender 1D solid elements approach in a 2D RVE (see Sect. Appendix: A.2). The model details are outlined below, starting with a brief description of the classical 2D multiscale model used as a foundation. Appendix: A.1 Classical computational multiscale modelling Appendix: A.1.1 Problem set-up Consider a macroscopic solid body representing the shock absorber material occupying a domain ⊂R2withasmooth boundary ∂, and material particles labeled by X.Let uM:→R2(A.1) be the macro displacement field, JM(X)=uM(X)⊗∇ X≡∇ XuM(X)(A.2) the corresponding displacement gradient tensor, and FM(X)=1+JM(X)(A.3) the deformation gradient tensor, where 1represents the second-order unit tensor. Here, each material point Xof the macroscale is associated with a representative volume element (RVE) of the microscale. Microscale variables will be denoted with the subscript μ, and macroscale with the superscript M. The RVE will be denoted by μ⊂R2, with a smooth boundary ∂μand its material coordinates by Y. For readability, the angle bracket •μwill denote the RVE volume average integral of the field (•), •μ≡1 μμ (•)dμ.(A.4) A kinematic connection between both scales will be established by considering the first-order Taylor’s expansion of 40 Where the spatial dimensions of the macro and micro scales are assumed the same 123 Computational Mechanics the kinematic variables associated with point Xin the macroscale. Thus, uμ(Y)=uM(X)+JM(X)·(Y−YC)+˜ uμ(Y), (A.5) ∇Yuμ(Y)=JM(X)+∇ Y˜ uμ(Y),(A.6) and Fμ(Y)=1+∇ Yuμ(Y)=FM(X)+∇ Y˜ uμ(Y),(A.7) where YCrepresents the coordinates of the centroid of the RVE and ˜ uμdenotes the micro-displacement fluctuation field. To ensure the consistency conditions, the following kinematic homogenization relations must be satisfied: 41 ˜ uμ(Y)μ=0(A.8) and ∇Y˜ uμ(Y)μ=0.(A.9) Or equivalently, (A.9) on the boundary, ˜ uμ(Y)⊗Nμμ=0.(A.10) By assuming conventional periodic conditions, the following expressions should be satisfied (see Fig. 19): NμY+=−NμY− ˜ uμY+=˜ uμY− uμ+=uμ−+FM−1·Y+−Y− (A.11) where the second term in the last expression remains constant for each pair of boundary points Y+and Y−. Appendix: A.1.2 Hill-Mandel principle The Multiscale Virtual Power Principle proposed in [23]is adopted in this work. It corresponds to a variational statement of the classical Hill-Mandel principle introduced in [27,28], in which an energetic equivalency between the two scales is established. It establishes that PM(X):δFM=1 μμ Pμ(Y):(δFM+∇ Yδ˜ uμ)dμ (A.12) 41 Equations (A.8)and(A.9) are often referred to in the literature as minimal kinematic restrictions. Here, they will be satisfied by assuming that the RVE fulfills the periodical conditions. Fig. 19 Scheme of a square periodic RVE cell ∀δFMand δ˜ uμ∈˜ Vμ, where ˜ Vμ:= δ˜ uμ∈H1(μ)|δ˜ uμμ=0;∇Yδ˜ uμμ=0 (A.13) denotes the virtual kinematically admissible fluctuation space.InEq.(A.12), PM(X)and Pμ(Y)refer, respectively, to the First Piola-Kirchhoff stress tensor at the macroscale point Xand at the corresponding RVE point Y∈μ(X). Equation (15) defines a fundamental link between the two scales. Direct consequences of it are: •The homogenized stress expression, PM=Pμuμ(Y)μ(A.14) •The RVE problem at the microscale, Find uμfulfilling Pμ(uμ):∇ Yδ˜ uμμ=0,∀δ˜ uμ∈˜ Vμ.(A.15) Finally, applying the general tensor relation described in [38] and the obtained strong form in Eq. (17), yields an alternative expression for the homogenized stress based solely on the boundary information, i.e.: PM=1 μ∂μ tμ⊗Yd∂μ(A.16) where tμ=Pμ·Nμis the Piola stress vector on the RVE boundary ∂μ. Appendix: A.2 From continuum to lattice micro-structures The microscale consists of a massive set of slender deformable solids arranged in a lattice structure, as illustrated in Fig. (20). As introduced above, the use of simplified 1D 123 Computational Mechanics degenerated kinematics leads to significant savings in computational costs but introduces an incongruity when developing the classical multiscale theory described in the previous section, in which the spatial dimensions of the macro and micro scales are assumed to be the same. To address this incongruity, the Ersatz material approach (widely used in the field of topology optimization, [39,40]) is temporarily applied. Under this approach, it is assumed that the RVE is composed of two material phases: a solid material phase, denoted as ωs μ(see Fig. 20), which represents the slender deformable solid of the microscale, and a substitute material phase, denoted as ωv, representing a fictitious low stiffness material that completely fills the initially empty regions of the RVE. Subsequently, the classical multiscale model summarized in the previous section can naturally be applied, since the macro and micro spatial dimension scales are the same. In this context, the kinematic constraint (A.9) can be rephrased as μ ∇Y˜ uμdμ=ωs μ ∇Y˜ uμdωs μ+ωv ∇Y˜ uμdωv μ=0, (A.17) or equivalently on the boundary as ∂ωs μ ˜ uμ⊗Nm μd∂ωs μ+∂ωv ˜ uμ⊗Nv μd∂ωv μ=0, (A.18) where Nm μand Nv μrefer to the outbound unit normals to the solid material phase and the fictitious phase, respectively. Taking into account that Nm μ=−Nv μ, the last condition (A.18) reduces to μ ˜ uμ(Y)⊗Nμdμ=0(A.19) where μ=∂ωs μ∩∂μ. Thus, the modified kinematically admissible fluctuation space is given by ˜ Vω μ=˜ uμ∈H1(μ)˜ uμωs μ =0; μ ˜ uμ(Y)⊗Nμdμ=0 .(A.20) Assume, without loss of generality, that ψμand ψv μdenote the strain energy densities of the solid and substitute phases of the RVE, respectively. Under the ersatz assumption, ψv μcan be expressed in terms of ψμas ψv μ=ψμ, where represents a very small positive scalar value. Taking into account that Pμcan be defined through the corresponding strain energy function, and making tend to zero, the equilibrium condition (17) can be reformulated as: Find uμfulfilling 1 μωs μ Pμ:∇ Yδ˜ uμdωs μ=0∀δ˜ uμ∈˜ Vω μ.(A.21) Unlike problem (17), the equilibrium problem (A.21) needs to be solved only within the material domain ωs μ. Following a similar procedure, the calculation of the homogenized stress tensor PMreduces to the following expression, PM=1 μμ tμ⊗Yd∂μ(A.22) Remark 2 For convenience, problem (A.21) can be formulated equivalently in terms of the second Piola-Kirchhoff stress tensor Sμand the Green-Lagrange deformations δEμ as follows: 1 μωs μ Sμ:δEμδ˜ uμdωs μ=0∀δ˜ uμ∈˜ Vω μ. (A.23) Appendix: A.2.1 RVE representation via a microlattice of slender buckling beams For computational savings reasons, the slender deformable solids, ωs μ, constituting the RVE will be modeled using 1D degenerate beam/frame kinematics. Assume that ωs μcan be represented by a straight beam network, μ(Figure 20b), modeled with the approach presented in Sect. 2. For simplicity, it is also assumed that the sections of each beam, Aμ, are uniform and that there are neither discrete forces nor discrete moments nor body forces applied on the microstructure. Therefore, the internal generalized forces remain uniform along each microstructural frame element. Thus, the mechanical variables at the macroscale are the position of the reference point, X, the spatial position of the point, x(X), the displacement vector, uM=x−X,the gradient of deformation, FM, and the Piola stress, PM. The microscale mechanical state, related to the macroscale point X, is defined by the position of the reference point, Y, the spatial position of the point, y(Y). Also, considering the beam theory of Sect. 2, the displacement vector of the beam neutral fibers is u0 μ=y(Y0)−Y0, the gradient of deformation, Fμ, and the Piola stress, Pμ. 123 Computational Mechanics Fig. 20 Representative RVE constituted equivalently by (1) a set of slender deformable solids and (2) a set of slender buckling beams Scale bridging equations The displacement vector uμof any point Yat the microscale, in agreement with 3, can be written as follows: uμ=uM(X)+JM·(Y0−YC)+˜ u0 μ(Y0)+ημ(˜ kμ−ˆ T) (A.24) where the condition that YCis the RVE mass center position results in: !JM·(Y0−YC)"=0(A.25) Under the hypothesis of periodic boundary conditions, the following mathematical expressions are fulfilled: u0 μY+=u0 μY−+JM·Y+−Y− ˜ u0 μY+=˜ u0 μY− ˜ θμY+=˜ θμY− (A.26) RVE problem at the beam microlattice Under this beam network representation, the variational problem stated in Eq. (A.23) can be rewritten as ωs μ Sμ:δEμdωs μ=μAμ Sμ:δEμdA μdμ(A.27) An equivalent formulation of (A.27) using generalized stresses and strains is presented below (refer to Appendix B for details) μAμ Sμ:δEμdA μdμ=μ zμ·δhμdμ(A.28) where zμ=(Nμ,Vμ,Mμ)Tis the generalized stress vector and δhμ=(δeμ,δγ μ,δκ μ)Tis the variation of the generalized strains vector defined in Eq. (7). From (A.28)the equivalent RVE problem can be stated as follows: Find (u0 μ,θ μ)fulfilling μ zμ(u0 μ,θ μ)·δhμδ˜ u0 μ,δ˜ θμdμ=0,∀δ˜ u0 μ∈˜ V0 μ, ∀δ˜ θμ∈˜ μ(A.29) where ˜ V0 μ=δ˜ u0 μ:μ→R2|δ˜ u0+ μ=δ˜ u0− μ,(A.30) and ˜ μ=δ˜ θμ:μ→R|δ˜ θ+ μ=δ˜ θ− μ(A.31) Homogenization Taking into account the contribution of the buckling beams of the lattice, the homogenized stress tensor in Eq. (A.22) can be expressed as: PM=1 μ iμiPμ·Nμ⊗Yidμ(A.32) where each integral represents the contribution of the i-th beam that intersects the boundary of the RVE. Without loss of generality, it is assumed that Nμcoincides with the normal vector of the cross-section, ˆ Ni. Consequently, by (1), the material points of the cross-section μican be expressed as 123