Continuum representation of nonlinear three-dimensional periodic truss networks by on-the-fly homogenization

Thanks to scale-bridging fabrication techniques, truss-based metamaterials have gained both popularity and complexity, ultimately resulting in structural networks whose description based on classical discrete numerical calculations becomes intractable. We here present a framework for the efficient and accurate simulation of large periodic three-dimensional (3D) truss networks undergoing nonlinear deformation (accounting for linear elastic beams undergoing finite rotations). Although the focus is on elastic beams, the method is sufficiently general to extend to inelastic material behavior. Our approach is based on a continuum representation of the truss (and its numerical implementation via finite elements) whose constitutive behavior is obtained from on-the-fly periodic homogenization at the microstructural unit cell level. We pursue a semi-analytical strategy (previously reported only in two dimensions) which admits the analytical calculation of consistent tangents for convergent implicit solution schemes; the extension to 3D – through the addition of torsional deformation modes and the handling of 3D rotations – results in a powerful tool for the prediction of the complex mechanical response of large structural networks. We validate the small-strain response by comparison to analytical solutions, followed by finite-strain benchmarks that compare simulation results to those of fully-resolved discrete calculations. The homogenization of beam unit cells results in a regularized macroscale model with an intrinsic length scale, which manifests especially when modeling bifurcations or localization. We finally apply our approach to macroscopic boundary value problems involving complex-shaped truss metamaterials (with truss unit cells near the body’s boundary mapped onto a conformal surface), which reveal only an insignificant effect of boundary layers on the overall mechanical response, again supporting the applicability of our homogenization approach. 2020 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
To describe the response of large periodic structures such as those exploited in mechanical metamaterials (Greer and Deshpande, 2019), a number of analytical, numerical and semianalytical techniques have been developed  with the aim of replacing expensive discrete calculations by an efficient continuum approximation of the underlying finescale structural deformation mechanisms in macroscopic boundary value problems. The accuracy and efficiency of such approaches vary tremendous, and they may generally be classified into techniques based on pre-computed or on-the-fly computed effective constitutive laws of the homogenized truss.
Pre-computed approaches employ, e.g., analytical continuum models based on affine homogenization of a unit cell (Triantafyllidis and Bardenhagen, 1993;Wang and Cuitiño, 2000), also available for beam networks with rotational degrees of freedom (Chen et al., 1998;Kumar and McDowell, 2004), as well as linear elastic (Deshpande et al., 2001) and higher-gradient elastic continua (Abdoul-Anziz et al., 2019). For more complex unit cell architectures that do not admit a simple micro-to-macro transition based on, e.g., affine unit cell deformations as well as for nonlinear problems, analytical approaches reach their limits, so that the unit cell problem is instead solved on-the-fly as in FE 2 -type simulations involving the numerical solution of a microscale (unit cell) problem (Kouznetsova et al., Jan 2001;Feyel, 2003). A popular approach is the application of periodic boundary conditions to the truss unit cell (Vigliotti et al., 2014;Pal et al., 2016), which is equivalent to a multi-lattice Cauchy-Born description (Glaesener et al., 2019). Alternatives are hybrid techniques such as the quasicontinuum https://doi.org/10.1016/j.ijsolstr.2020.08.013 0020-7683/Ó 2020 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). method (Phlipot and , which localizes discrete calculations at full truss resolution in regions of interest while introducing an efficient continuum approximation elsewhere. Unfortunately, such on-the-fly approaches tend to suffer from high computational expenses, which is why much of previous research has been limited to two-dimensional (2D) problems; see, e.g., Vigliotti et al. (2014), . Although studies of the effective response of three-dimensional (3D) trusses date back to, at least, Nayfeh and Hefzy (1978) and have gained significant insight into, e.g., the popular octet truss (Deshpande et al., 2001;Challapalli and Ju, 2014), those approaches were limited to small strains and did not account for nonlinear behavior. Various studies have used finite-strain models to predict the nonlinear behavior of truss networks undergoing uniform deformation (see, e.g., Mohr (2005) for predictions of the plastic response of trusses, or Goda et al. (Jan 2014 for the homogenized response of truss and fiber networks), yet those were too expensive for twoscale simulations and hence not applicable for heterogeneous macroscale deformation.  extended their FE 2 approach to 3D, yet computational expenses limited its application to relatively simple benchmarks without significant nonlinearity or localization. Especially the latter scenarios benefit from a model with an intrinsic length scale, which has been achieved, e.g., by higher-order homogenization schemes (Kouznetsova et al., 2002;Kouznetsova et al., 2004).
In Glaesener et al. (2019), we recently introduced a 2D continuum approach for the simulation of periodic trusses undergoing finite deformations (considering beam networks under large rotations), which captured nonlinear behavior, localization, and bifurcation upon loading. We here extend that formulation into a versatile 3D framework for the efficient simulation of macroscale boundary value problems. Moreover, we extend our concept to arbitrarily shaped bodies with an underlying periodic truss architecture, and we examine the influence of boundary layers in case of macroscopic bodies that do not conform with the underlying truss. Our two-scale description of beam networks (based on onthe-fly homogenization) naturally introduces an intrinsic length scale and hence regularizes localization modes for convergent, mesh-independent numerical solutions. We validate the model by comparing simulation results to those of discrete numerical calculations for a series of 3D benchmarks, and we demonstrate the benefits and performance of the proposed scheme for examples of homogeneous deformation, localization, and buckling. Of course, in the linearized limit of small deformations, the continuum model converges to that of a homogenized linear elastic solid, which admits verification by classical examples from the literature (Deshpande et al., 2001;Valdevit et al., 2013). Importantly, the fact that a consistent tangent is available despite the numerical solution of a nested microscale problem allows for efficient macroscale calculations based on implicit solvers for the quasistatic boundary value problems of interest. We show that our 3D on-the-fly homogenization framework performs well in comparison to exact discrete truss simulations, serving as a computationally tractable alternative for the simulation of periodic truss architectures and truss metamaterials.

Two-scale description of truss lattices
We base our description of slender beam networks on the corotational beam formulation of Crisfield (1990), which assumes linearized strains but accounts for finite rotations. For each beam segment, the deformation is uniquely characterized by the (generalized) degrees of freedomx ¼ x; h ð Þ T at each of its two end points, h 3 Þ T the rotation of each end point from its reference configuration (denoting Euler angles in 3D space). Although we focus in the following on linear elastic beams, the overall setup and homogenization approach is sufficiently general to apply to other beam formulations and, in particular, to inelastic behavior through variational constitutive updates (see, e.g., Lestringant et al. (2020) for suitable viscous and viscoelastic beam formulations). We aim to describe truss architectures made by periodically tessellating a unit cell X UC within a macroscopic body X & R 3 of arbitrary shape. To avoid the expenses of discrete numerical calculations which resolve each and every beam within the body, we resort to the two-scale description schematically shown in Fig. 1. We replace the truss network by a continuous body on the macroscale (discretized by finite elements in order to numerically solve macroscopic boundary value problems), while the local constitutive behavior of the macroscale body emerges from on-the-fly solving a periodic boundary value problem on the microscale, i.e., at the unit cell level. Unlike classical FE 2 schemes , we couple both translational and rotational degrees of freedom between the micro-and macroscale. Beams on the microscale naturally possess rotational degrees of freedom, while on the macroscale a new rotation field is introduced like in micropolar continua (Eringen, 1965;Kumar and McDowell, 2004;Goda et al., Jan 2014). This introduces an intrinsic length scale via the unit cell size, as will be discussed in the following.

Formulation on the macroscale
By extending the 2D framework of Glaesener et al. (2019), we describe the macroscopic deformation of a body X 3 made of a periodic truss by two continuous fields: the deformation mapping u : X ! R 3 (describing the deformed position of an undeformed material point X 2 X as x ¼ u X ð Þ) and a rotation field h : X ! 0; 2p ½ 3 (describing local rotations through the three angles h 1 ; h 2 ; h 3 f g ). For convenience, we combine the latter two into the generalized deformation mapping containing three translational and three rotational degrees of freedom. From the various available parametrizations of rotations in 3D we choose Euler angles, since those admit an additive composition of incremental updates in the numerical treatment. Introducing the deformation gradient F ¼ Gradu and the curvature tensor j ¼ Gradh (with gradients defined with respect to the undeformed basis) leads to the generalized deformation gradient With the above kinematic description, we choose a variational setting based on the total potential energy with a strain energy density W and generalized surface tractionsT (containing both tractions and moments) acting on the Neumann boundary @X N & @X. For generality, we assume the energy density W depends on both the generalized deformation mapping and its gradient (which will become beneficial in the microscale description detailed below). Unlike, e.g., in micropolar or Cosserat theories, we do not postulate a particular functional form of the effective energy density W nor do we assume any a-priori coupling between the kinematic fields of translations and rotations. Instead, Wũ;F will be computed on the fly from a representative unit cell on the microscale, whose deformation is imposed to equalũ andF on average. The exact microscale treatment will be detailed in Section 2.2.
Knowledge of W ¼ Wũ;F is sufficient to compute the conjugate stress measures, e.g., Static equilibrium of an elastic medium (in the absence of body forces) implies stationarity of the total potential energy functional (3). This is equivalent to the balance equation and boundary condition, respectively, with unit outward surface normal N (and angular momentum balance is satisfied automatically by material frame indifference of W, which is invariant under rigid-body rotations). The resulting boundary value problem on the macroscale is solved numerically by the finite element (FE) method, which discretizes both translational and rotational fields with the same hexahedral (bi-linear) or quadratic tetrahedral (linear-strain) elements and solves the resulting weak form by implicit iterative techniques. We refer to Appendix A for a derivation of the weak form and a summary of the numerical implementation as well as to Glaesener et al. (2019) for a detailed derivation considering higher-order gradients (but restricted to 2D).

Formulation on the microscale
Our goal on the microscale is the identification of the effective energy density Wũ;F required by the continuum model on the macroscale, as introduced above. To this end, we impose a periodic deformation onto a representative unit cell (RUC), such that the RUC-average deformation matches a given generalized deformation mappingũ and gradientF. In practice,ũ X ð Þ andF X ð Þ will be passed from the macroscale to the microscale for every macroscale point X 2 X, and the microscale calculation returns Wũ X ð Þ;F X ð Þ and its gradients, as needed.
The RUC on the microscale is composed of a periodic arrangement of corotational beams (Crisfield, 1990), whose deformation is defined through the translations and rotations of all nodes (i.e., the junctions and end points of all beam elements). We assume a separation of scales between variations of all fields in the macroscale body X and the size of the RUC X m . This allows us to formulate the deformation within X m by considering a Taylor expansion of the deformation in an infinitesimal neighborhood of a macroscale point X 2 X. (Here and in the following, the superscript m indicates a microscale quantity.) If the RUC contains a total of n nodes located initially at X m a 2 X m (a ¼ 1; . . . ; n), their deformed positions x m a and rotations h m a follow from Taylor expansions as with the generalized shiftsd i ¼ dx i ; dh i ð Þ T . We here ignore higherorder gradients based on the observation of Glaesener et al. (2019) that retaining second-order terms in the Taylor expansions in (6) has a vanishing impact in most practical scenarios.
Since an affine deformation of the RUC yields an unreasonably stiff response (especially when dealing with bending-dominated truss topologies), we introduced in (6) translational and rotational shifts dx L a ð Þ and dh L a ð Þ , respectively, which enforce periodic boundary conditions. If the n nodes inside X m are decomposed into m simple Bravais lattices, then each Bravais lattice is characterized by unique constant shifts dx i and dh i (i ¼ 1; . . . ; m). For example, a simple cubic lattice is a simple Bravais lattice so no shifts are required, whereas a body-centered lattice is composed of two simple Bravais lattices, so one set of shifts is required to describe the relative displacements and rotations between the two lattices.
For any lattice, L a ð Þ 2 1; . . . ; m ð Þdenotes the Bravais lattice that node a is associated with. The expansions in (6) ensure that the RUC averages match the macroscale kinematics, i.e., as long as the RUC is centrosymmetric and if We note that the definition of the RUC averages in (8) is not unique (and could, e.g., be replaced by weighted averages involving masses or volumes); see Glaesener et al. (2019) for a discussion.
Based on the RUC deformation defined by (6) and introducing the we define the volume-averaged RUC strain energy as Fig. 1. The periodic truss is replaced by a continuous body on the macroscale, whose mechanical response is simulated numerically based on a finite element discretization. The effective constitutive behavior of this body (used within every finite element integration point) is extracted from a discrete representative unit cell (RUC) on the microscale, consisting of struts modeled as corotational beams. While the macroscale passesũ andF to the microscale, the RUC averages of energy, stresses, and consistent tangent are returned to the macroscale.
where W m is the energy of the mth beam in the RUC with B m denoting the set of all beams within the RUC on the microscale, and x m m1 ;x m m2 È É are the generalized degrees of freedom (translations and rotations) of the two nodes defining the mth beam, computed via (6). Derivatives of (10) with respect toũ andF can be computed analytically (Glaesener et al., 2019) and involve the nodal reactions of beams within the RUC on the microscale (which describes each strut as a corotational beam following Crisfield (1990Crisfield ( , 1991); see Appendix A for details. In equilibrium, the shiftsd are minimizers of the total potential energy. Since the energy is local in the shifts, that minimization can be carried out at the RUC-level, so we define the effective strain energy density of the truss lattice as in which the constraints (9) are to be enforced when computing the minimizing shifts. Due to the corotational beam formulation, this minimization is generally a nonlinear problem that requires an iterative optimization scheme. We here use Newton-Raphson iteration to obtain the minimizing shifts. Note that, although our examples here are linear elastic, the above setup also applies to inelastic beams, as long as a variational formulation is available. The effective energy density (11) enters the macroscale problem in (3) and (5), which completes the two-scale description. Since the macroscale energy density W in (3) equals the average RUC energy on the microscale and all macroscale stress measures are derived from W directly without averaging any stress components on the RUC-level, this micro-to-macro transition satisfies Hill's macrohomogeneity condition exactly. The derivatives @W=@ũ and @W=@F (which are required in the macroscale problem) are trivial to calculate, since the shifts are minimizers of W m , so @W=@ũ ¼ @W m =@ũ and @W=@F ¼ @W m =@F, see Appendix A. Calculating second derivatives (required for the consistent tangents of iterative solvers on the macroscale) are significantly more involved since @D=@ũ-0 and @D=@F-0, and shifts must be computed numerically. Yet, the analytical consistent tangents derived in Glaesener et al. (2019) for 2D can be employed here, which is beneficial for efficiency (avoiding costly numerical derivatives that are oftentimes employed in FE 2 -type simulations); see Appendix A for details. In addition, the small number of nodes and hence of degrees of freedom per unit cell renders this two-scale framework sufficiently efficient to be applied to complex 3D problems -examples of which will be presented in the following.

Small-strain validation
To validate the above framework in the linear regime for which analytical solutions are available (Zhu et al., 1997;Deshpande et al., 2001), we compute the full effective stiffness tensor for three frequently used 3D truss unit cells: the octahedron, bitruncated octahedron, and octet (which will also be used in subsequent numerical simulations in Section 3). For all unit cells we choose a linear elastic base material with Young's modulus E and Poisson's ratio of m ¼ 0:3. Beams are assumed to have a constant circular cross-section with a slenderness ratio of k ¼ D=L ¼ 0:0522, where D denotes the beam diameter and L its length.
The components of the fourth-order stiffness tensor are derived from the effective energy density (11) as For a comparison with the linear elastic moduli of Zhu et al. (1997) and Deshpande et al. (2001) (computed in linearized kinematics), we calculate the effective linear elastic moduli as . Appendix B summarizes the elastic moduli obtained from this approach in comparison to those of Zhu et al. (1997) and Deshpande et al. (2001). Errors stem primarily from the numerical shift calculation but are less than 1% throughout, which we deem satisfactory.

Results
To demonstrate the applicability and accuracy of our two-scale simulation framework for truss networks, we present numerical examples of quasistatic 3D boundary value problems -for each of which we compare the mechanical response as obtained from the two-scale approach to expensive discrete numerical calculations of the underlying structural network. From the truss topologies tested, the bending-dominated bitruncated octahedron is most complex due to the importance of the microscale shift calculations, while the stretching-dominated octahedron and octet underline the applicability to large, nonlinear deformations. Beyond visual inspection, we examine the effective load-strain (or torque-twist) response of loaded samples to assess the accuracy.
All subsequent examples assume beams of circular crosssections with diameter D, area A ¼ pD 2 =4, area moment of inertia I ¼ pD 4 =64, and polar moment of inertia I p ¼ 2I. Based on the individual beam length L, we define the beam slenderness ratio k ¼ D=L. While an extension to more complex material behavior is, in principle straight-forward (as long as a variational formulation is available), we here restrict ourselves to an isotropic linear elastic base material having Young's modulus E, Poisson's ratio Þ . (We note that linearized strains do not rule out large deformations, as finite rotations are accounted for.) The reported relative density q of each truss is defined as its fill fraction, i.e., the volume of all struts within a RUC divided by the volume of that RUC.
The choice of the RUC for each example depends on the problem: while the smallest primitive unit cell of each lattice is sufficient for capturing homogeneous deformation, the emergence of bifurcation patterns at larger loads requires larger RUCs (composed of multiple cells). Specific RUCs will be reported below for each example, for which we know a priori (e.g., from Bloch analysis) the critical bifurcation mode(s) that are expected under loading. Such knowledge allows us to select a sufficiently large RUC to capture all relevant bifurcations. Of course, this excludes larger bifurcation patterns that cannot be represented by the chosen RUC, which is an inherent limitation of the chosen homogenization approach. (One could use on-the-fly Bloch analysis to adaptively adjust the RUC size, yet we here restrict our analysis to fixed-size RUCs.) We solve the macroscale FE problem using hexahedral and quadratic tetrahedral meshes along with the Toolkit for Advanced Optimization (TAO) from the Portable, Extensible Toolkit for Scientific Computation (PETSc) (Abhyankar et al., 2018;Balay et al., 1997;Balay et al., 2019;Balay et al., 2020).

Bending-dominated bitruncated octahedral truss under torsion
The bitruncated octahedron is a low-density, strongly bendingdominated truss topology, which implies that the expected deformation on the microscale is non-affine (thus highlighting the importance of the periodic Bravais shifts introduced within the RUC). Applying torsion to a bitruncated octahedron lattice results in a complex, inhomogeneous deformation which we utilize to highlight the performance of our multiscale framework. Fig. 2 shows a twisted cube containing 50 Â 50 Â 50 unit cells at a twist angle of 45 . Bottom and top surfaces are forced to deform in an affine manner to impose the twist while being clamped vertically. The chosen slenderness ratio of k ¼ 0:127 allows for good convergence of our two-scale framework. We solve the FE problem using a hexahedral mesh. Discrete simulations of the bitruncated octahedron cube are expensive (involving approximately 1; 530; 000 nodes and 3; 030; 000 beam elements) but within reach, thus providing a means for validation.
The torque-vs.-rotation curves in Fig. 2 show convincing agreement between the discrete simulation and the two-scale approximation, as long as we allow the microscale to deform in a periodic fashion through shifts. Also shown are results for affine and partially relaxed (only translational shifts) RUC deformation, which yield a significantly stiffer response. As can be expected, the best performance is achieved when both translational and rotational degrees of freedom are periodic at the RUC-level (the latter gaining importance with decreasing beam slenderness). We note that, in case of highly slender beams, the body behaves in an approximately incompressible manner, in which case the macroscale problem is better addressed by a suitable FE discretization based on augmented elements for modeling incompressible material behavior. Fig. 3 illustrates the convergence of our two-scale simulation results with h-refinement in the macroscale FE problem (going from 100 to 625 finite elements on the macroscale). Shown is the error of the measured net torque on the top and bottom faces of the cube at a twist angle of 45 . Rather than demonstrating volume-averaged error metrics, this serves as a particular worstcase scenario since boundary effects and associated localized deformation modes at the clamped top and bottom faces can considerably affect the torque, as highlighted in the mesh displayed in Fig. 3 showing differences in the deformation at the corners or the discrete and the FE mesh.  Bertoldi et al. (2010) and similar in spirit to the 3D designs of Shen et al. (2014Shen et al. ( , 2017Shen et al. ( , 2018, who studied connected arrays of solid cubes and prisms undergoing similar bifurcations. As a consequence of this bifurcation pattern, the truss is expected to display auxeticity (Evans et al., 1991;Reda et al., 2018) in all principle directions at large compressive strains, which we here demonstrate. As described in Overvelde et al. (2012), this concept can be extended to various geometrically more complex beam assemblies.

Instabilities in a reinforced octahedral truss under compression
The RUC in the two-scale approach is chosen to contain 2 Â 2 Â 2 unit cells to allow us to capture the bifurcation pattern . (As shown in Fig. 4, the octahedral units fold under compression to gradually collapse into the void spaces they enclose.) To improve numerical convergence, we introduce small initial imperfections in the form of perturbations of all undeformed nodal locations of magnitude 10 À4 L and in the shape of the expected bifurcation mode. We impose affine displacements across the top surface and fix the bottom surface to enforce a compressive strain, while leaving the rotations unconstrained and measuring the resultant forces (as before using a hexahedral mesh). Comparing the obtained effective force-strain curve to that of a discrete numerical calculation (for a beam slenderness of k ¼ 0:05) reveals good agreement across the full nonlinear deformation regime (including the bifurcation under compression), as shown in Fig. 4. The bifurcation appears (slightly) earlier for the fully resolved truss than for the homogenized model, which we attribute to surface effects (free boundaries are captured only in an average sense in the two-scale approach). Fig. 5 demonstrates convergence with hrefinement of the macroscale mesh in the two-scale simulation. Shown is the relative error of the total compressive force after bifurcation (at a compressive strain of 0:6%).
As a metric of practical interest, we quantify the Poisson effect by measuring the ratio of (observed average) lateral strain to compressive strain of the reinforced octahedral truss under compression, m ¼ DL z =DL x with L x and L z denoting the side lengths of the deformed cube (L z being imposed and L x ¼ L y measured as average across all boundary nodes in the respective directions). Fig. 6 plots m vs. the compressive strain for both the homogenized and discrete simulations, both demonstrating auxeticity beyond a compressive strain of about 0:4% (and differences in the two curves stemming from the shifted onset of instability already seen in Fig. 4). Fig. 2. Simulated torque vs. applied twist angle for a cube made of 50 Â 50 Â 50 bitruncated ocathedral unit cells, comparing the exact structural response to the two-scale approximation. For illustration purposes, we present results obtained by considering a) no unit cell shifts (affine RUC behavior), b) only translational (reduced) shifts, and c) both translational and rotational (complete) shifts. Shown on the right and above are visual comparisons of the resulting deformed truss and its FE counterpart (at a twist angle of 45 ). The torque is normalized by the base material's Young's modulus E, the side lengthL of the overall cube, and the relative density of the truss (q ¼ 1:35%).

R.N. Glaesener et al.
International Journal of Solids and Structures 206 (2020) 101-113 Fig. 3. Convergence of the relative error in the net torque at a fixed twist angle of 45 vs. the (inverse) finite element size of the macroscale boundary value problem (L is the side length of the cube; the exact solution is taken as the result of a discrete truss simulation; the error is defined as the relative deviation from the result of a discrete simulation). The magnification of one of the top corners shows good agreement but also points out deviations in such regions of strong deformation, thus over-predicting the stiffness in these areas. Fig. 4. Tension-compression of a cube made of 30 Â 30 Â 30 reinforced octahedral (body-centered cubic) truss unit cells, simulated by the two-scale homogenization approach (based on a 2 Â 2 Â 2 RUC) as well as by a discrete numerical calculation. Under compression, the structure shows an auxetic bifurcation pattern. The net compression force is normalized by the base material's Young's modulus E, the side lengthL of the cube, and the relative density q ¼ 1:13%. The discrete and two-scale FE cubes are visualized at a compressive strain of 1%. Insets show the RUC in its undeformed configuration and at a compressive strain of 5% (this large latter strain level, not contained in the plot, is chosen to emphasize the auxetic bifurcation pattern.

Fig. 5.
The relative error of the total compressive force after bifurcation (at a compressive strain of 0:6%) decreases with decreasing unit cell size in the discrete truss (i.e., keeping the overall truss sizeL constant while changing the number of unit cells within, n). The number of finite elements in the two-scale approach has no effect on the result, since the deformation is homogeneous. Included in the log-log plot is a linear regression fit. from the two-scale approach and an exact truss calculation (using a 30 Â 30 Â 30 truss). Both models reveal the transition to auxeticity upon bifurcation. The small offset between the transition strains is attributed to boundary effects in the discrete structure, which decrease with increasing unit cell size (Fig. 5).

Instabilities in an octet truss under tension and compression
While the above example displayed a bifurcation mode based on the rotational collapse of close-to-rigid octahedral units, we here demonstrate buckling instabilities in an octet lattice which engage individual beams in buckling, as shown in Fig. 7. Being a combination of octahedron and tetrahedron, the octet produces a buckling pattern which resembles that of the triangular 2D lattice (Glaesener et al., 2019). To observe such patterns under vertical compression, we rotate the unit cell by 45 about two principle axes, resulting in the truss shown in Fig. 7. To capture the expected short-wavelength instabilities, we refine each strut into two beam elements (thus adding one additional node at the center of each strut), resulting in a total of 5; 658 beam elements in the discrete truss. Akin to the 2D triangular example, the stiff 3D tetrahedral truss softens considerably under both tension and compression upon bifurcation. The two-scale FE problem again employs a hexahedral mesh. Fig. 7 illustrates the two distinct bifurcation modes that engage in compression and tension (for a beam slenderness ratio k ¼ 0:1). To capture such behavior in our two-scale approach, we choose a RUC consisting of four tetrahedra and two octahedra. (The resulting RUC is identical to that of an octahedron rotated by 45 about two principle axes.) As in the example of Section 3.2, we impose vertical displacements on the top and bottom nodes, while letting the latter free to slide in-plane. The results summarized in Fig. 7 confirm a convincing agreement between the homogenized twoscale simulation results and the discrete truss calculation in both tension and compression.

Non-confirming macroscale truss shapes and surface effects
So far we have restricted all validation examples to simple macroscopic bodies whose shapes conform with the underlying truss unit cells. For practical relevance, we proceed to consider generally shaped macroscopic bodies that do not conform with the underlying truss architecture, so the periodic assembly of unit cells inside the 3D body inevitably leads to geometric approximations near the surface (and non-unique choices of how to adjust unit cells near the surface) (Wu et al., 2019) as well as to approximation errors due to the FE discretization in the two-scale approach. Of course, this error vanishes in the limit of a true separation of scales, yet real-world examples always deal with finite sizes of the macroscopic body and the underlying unit cell. We here aim to assess the applicability our two-scale approach for this scenario and to elucidate the influence of surface effects.
We fill a pre-defined (ideal) macroscale body X (with boundary @X) by tessellating a given unit cell X UC throughout the body, starting at a random point inside X and retaining each unit cell if its geometric center falls into X and deleting unit cells otherwise (near and outside the boundary). This procedures results in steps on the surface of X, as shown in Fig. 8(a). The solid (blue) line illustrates @X, the boundary of the macroscopic body used in the twoscale finite element simulation, whereas in the discrete truss some nodes and struts are located inside X (marked by solid green circles) and others outside (open red circles). Consequently, the truss is perfectly periodic but the boundary @X is not reproduced exactly, which we refer to as a serrated surface. Alternatively, we may adjust the boundary unit cells by breaking the periodicity in a boundary layer on the order of the unit cell size. With the goal of creating a conforming surface from the periodic truss (i.e., modifying the periodic truss unit cells such that the target outer surface of the macroscopic body is the union of truss cell faces), we use the two approaches sketched in Fig. 8.
The first method, sketched in Fig. 8(b), projects (pushes) nodes from outside X along existing beams onto @X (within a search region indicated by the dashed red line). The second method, see Fig. 8(c), additionally relocates nodes from the inside (within the same search distance from the boundary, shown as dotted green line) onto the boundary. This latter technique involves a number of design rules such as the merging of overlapping nodes, removal of beams shorter than L=4 (with L being the smallest beam length in the unit cell), and rearrangement of beams to improve triangulation aspect ratios. This approach generally leads to improved Fig. 7. Tetrahedral truss showing bifurcation under both tension and compression due to short-wavelength buckling modes: average normalized stress vs. strain as obtained from both the homogenized two-scale approach and a discrete truss calculation, along with illustrations of the buckling patterns that appear in the discrete structure (outer images) as well as in the RUC of the homogenized model (insets). The average compressive stress (F=L 2 ) is normalized by the base material's Young's modulus E and the relative density q ¼ 6:67%. Images show axial strains inside the struts (neglecting bending contributions) to visualize beams undergoing compression vs. tension.
R.N. Glaesener et al. International Journal of Solids and Structures 206 (2020) 101-113 truss connectivity near the boundary. Details on the periodic mesh generation inside an arbitrary bounding shape can be found in Appendix C. The methods of Fig. 8 allow us to assess the influence of boundary effects arising from macroscale bodies that are nonconforming with the underlying truss unit cell. The two-scale approach is independent of those methods and simulates the ideal body X -here discretized by quadratic tetrahedral elements (meshes are generated using Gmsh (Geuzaine and Remacle, 2009)). By contrast, the response of the discrete trusses is expected to differ near the surfaces based on the above choices.
As an example, we consider a cylindrical macroscopic body filled by an octet truss with 50 unit cells each across the diameter and height of the cylinder (we assume circular beam cross-sections with a slenderness ratio of k ¼ 0:1). The cylinder is twisted by affinely rotating the top and bottom surfaces against each other. Fig. 9 compares the net torsional moment vs. twist angle, as obtained from the two-scale approach and from discrete calculations -the latter using structures based on the surface treatments of Fig. 8. Results overall show good agreement, considering that the smooth radial surface of the cylinder in the two-scale simulation is a rough approximation of the discrete structure. The specific surface adaption method used appears to have only a small effect on the effec-tive response, underlining that, for a sufficiently small unit cell, changes to the surface unit cells only insignificantly affect the mechanical response. The torque computed by the two-scale simulation serves as an upper bound to all exact simulations (using the three (un)corrected tessellations of Fig. 8), as may be expected Fig. 8. Turning a periodic truss into a conforming macroscopic body through truss modifications near the surface X of the macroscopic body (the sought target boundary of the macroscopic body is shown as a solid blue line). Shown is the example of a cross-reinforced square unit cell (schematically shown in 2D but analogously applicable in 3D), which is used to fill the body. (a) Unmodified, serrated surface including only perfect unit cells, which do not conform with the boundary @X. (b,c) Unit cells near the surface are modified by manipulating node positions to conform with the target boundary @X. (Not shown is the subsequent removal of overlapping nodes and beams as well as the fusion of short beams in (c).). Fig. 9. Computed net torque vs. twist angle for a cylindrical truss made of octet unit cells, comparing the two-scale approximation to the exact, discrete response of a fully resolved lattice with a serrated surface and also with a surface that has been manipulated using the methods of Fig. 8. Torques are normalized by the base material's Young's modulus E, height L and diameterD of the cylinder, and the relative density q ¼ 6:66% (not accounting for surface corrections). R.N. Glaesener et al. International Journal of Solids and Structures 206 (2020) 101-113 from the stiffening kinematic assumption in the homogenization approach.
To illustrate the vanishing impact of surface modifications with decreasing unit cell size, Fig. 10 shows the deviation in the net torque between cylinder trusses with serrated surfaces vs. modified confirming surfaces (using the method of Fig. 8(c)). As expected, with decreasing unit cell size (while keeping the overall cylinder size fixed), the deviation between the two computed torques (at a twist angle s ¼ 45 ) decreases and boundary effects lose importance.
Twisting of a soft cylinder may lead to torsion instability, which produces periodic shear bands on the cylinder's surface. Twisting our truss-filled cylinder is analogous to the torsion of soft homoge-   R.N. Glaesener et al. International Journal of Solids and Structures 206 (2020) 101-113 neous cylinders, for which a torsion instability at critical twist angle was observed, e.g., in soft stubby hyperelastic cylinders (Ciarletta et al., 2014). This phenomenon is indeed displayed in Fig. 12, which shows the discrete and FE cylinders form Fig. 9, here plotting the distribution of the vertical (normalized) displacement at three twist angles. The low relative density of the truss is responsible for its relatively soft behavior, which is why we observe the torsion instability -captured by both the discrete truss and its two-scale approximation. The strain distribution in Fig. 12 (and its surface projection in Fig. 13) reveals a buckling mode with circumferential wave number m ¼ 4.
Localizations of the above type as well as shear bands (previously reported in 2D examples (Glaesener et al., 2019)) require some means of regularization, so that the numerical scheme on the macroscale yields mesh-independent solutions. Our framework establishes the size of the unit cell as a natural intrinsic length scale of the two-scale model, which provides regularization through the Euler-Bernoulli beam formulation at the microscale and its coupling of deflections, rotations, and curvature. This becomes essential when instability and/or localization occurs, such as in Fig. 13. Fig. 11 demonstrates the mesh independence of the computed nonlinear torque-twist relation for the twisted cylinder, comparing results obtained from five different resolutions of the macroscale mesh (whose deformations are illustrated in Fig. 14 for some of the meshes). Despite the emergence of complex localized patterns, the shown response is mesh-independent, thanks to the regularizing influence of the micropolar-like two-scale formulation, as already observed in Glaesener et al. (2019).

Conclusions
We have outlined a two-scale computational framework for the simulation of large-scale truss networks, which replaces the discrete truss architecture by a continuous body on the macroscale whose effective mechanical constitutive behavior is computed on-the-fly by periodic homogenization on the microscale based on a representative unit cell. We have shown that this setup, which formulates the macro-problem based on both a deformation mapping and a rotation field, accurately captures the nonlinear deformation behavior of a variety of stretching-and bendingdominated 3D truss networks, including bifurcation and localization as well as regular and irregularly-shaped macroscopic bodies (the latter case involving boundary layers of incomplete or imperfect unit cells). We reported good qualitative and quantitative agreement, thus providing an efficient alternative to full-scale structural simulations of use, e.g., in predicting the complex response of mechanical metamaterials. We point out that the chosen homogenization approach depends critically on a reasonable choice of the representative unit cell on the microscale. Therefore, it only applies as long as instabilities can either be captured at the local unit-cell level (which is why we chose the unit cell sufficiently large in our examples with a-priori knowledge about the expected bifurcation modes) or at the large, long-wavelength level (where the macroscale finite element model captures those). There is an intermediate range of wavelengths between those two regimes, where the homogenization approach cannot capture instabilities, unless, e.g., Bloch analysis is used to adaptively select an appropriate unit cell size on the fly.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
With the concise notation introduced in Section 2.1, the variational problem on the macroscale becomes ðA:1Þ with generalized tractionsT applied on the Neumann boundary @X N and essential boundary conditionsũ ¼û imposed on the Dirichlet boundary @X D . Note that the effective energy density, defined in (11) That is, W is the average RUC energy minimized with respect to the shiftsD, andD Ã denotes the minimizing shift vector. This RUClevel minimization with respect toD must be carried out numerically due to the nonlinear nature of the problem and is in practice solved by Newton-Raphson iteration at the RUC-level (Glaesener et al., 2019).
The weak form is obtained by taking the first variation of (A.1) in the direction of an infinitesimal admissible perturbation fieldg, which leads to

ðA:3Þ
We discretize the macroscale fields by introducing a finite element mesh with a total of n s nodes, and we define shape functions N a X ð Þ (a ¼ 1; . . . ; n s ) to interpolate the nodal degrees of freedom.
Condensing all nodal unknowns intoŨ h ¼ũ 1 ; . . . ;ũ ns È É and analogously E h ¼g 1 ; . . . ;g ns È É lets us define the approximate fields in Bubnov-Galerkin fashion as ðA:5Þ f a int;i andf a ext;i denote, respectively, the (generalized) internal and external force vector components associated with node a (including both forces and torques). We solve the nonlinear system (A.5) by Newton-Raphson iteration, for which the components of the associated consistent tangent matrix arẽ with the local stiffness tensors ðA:7Þ Note that, as discussed in Glaesener et al. (2019), the definition of the effective energy in (A.2) requires caution when differentiating. For example, the stresses can be computed by straightforward differentiation since, ðA:8Þ and the other stress measures follow analogously. By contrast, the analogous does not hold for second derivatives (the second term does not vanish), so that, e.g.,

ðA:9Þ
To compute the last term in (A.9) we may use ðA:10Þ which admits the analytical calculation of the tangent (A.9) despite the iterative numerical solution forD Ã . Practically, W m is computed at the RUC-level by calculating the volume-averaged energy of all structural members within the RUC based on the corotational beam formulation of Crisfield (1990). Starting from the average RUC energy (10), all stress and stiffness measures are computed by differentiation. For example, we have  less than t 8 from the continuous line, the respective voxel is set to 1 resulting in the voxel matrix in Fig. C.15(b). For volumetric voxeling, we utilize the mesh generator Gmsh (Geuzaine and Remacle, 2009), which generates a volume-filling tetrahedral mesh inside any tringulation (defined in an STL file). The voxel matrix contains all geometric information necessary for a discrete tessellation. Fig. C.16 sketches the procedure to translate the voxel matrix into a periodic truss, with special attention required for avoiding the duplication of struts shared across unit cells while maintaining full struts in unit cells on the outer boundary. (The unit cell in this 2D example is a reinforced square; yet, the method is applicable to any periodic unit cell geometry with orthogonal translation vectors). The number of unit cells in the final assembly matches the sum over all voxels.