Helen Le Clézio

Claire Lestringant

Dennis M Kochmann

A numerical two-scale approach for nonlinear hyperelastic beams and beam networks

Keywords: Structure, Dimension Reduction, Hyperelasticity, Metamaterial, Finite Element Method

We introduce a numerical framework for modeling hyperelastic slender trusses, oftentimes used as elementary building blocks in architected materials, which accounts for both the geometric nonlinearity inherent in thin structures and the nonlinear constitutive behavior of the base material. Akin to the F E 2 method in homogenization, our approach is based on a formal, two-scale expansion. We decompose the three-dimensional (3D) description of a slender structure into a macroscale problem solving for the deformation of the beam center-line (based on an effective strain energy density, which depends on the stretching, bending, and torsional strains of the beam's center-line) and a series of two-dimensional (2D) microscale boundary value problems (defined over the beam cross-section). We solve a series of 2D problems (covering a range of macroscopic strain combinations) for a given cross-sectional geometry and material distribution in an offline, pre-processing step. Using this pre-computed energy landscape, we solve the macroscopic boundary value problem within a geometrically exact, nonlinear discrete beam framework. We demonstrate the accuracy of this approach through a set of benchmark problems highlighting the nonlinear effects of the cross-sectional geometry and constitutive material, including material heterogeneity and pre-strains. We further illustrate how this technique is applied to truss-based architected materials consisting of networks of slender struts, which are typically made of polymeric base materials and thus allow for large nonlinear elastic deformation.

Introduction

Advances in additive manufacturing (including jet-printing and fusion deposition melting) have enabled the fabrication of a wide range of lightweight architected materials (or metamaterials), composed of a variety of base materials from metals to polymers and elastomers [START_REF] Leary | Inconel 625 lattice structures manufactured by selective laser melting (slm): Mechanical properties, deformation and failure modes[END_REF][START_REF] Glaesener | Viscoelastic truss metamaterials as time-dependent generalized continua[END_REF][START_REF] Jiang | Highly-stretchable 3D-architected Mechanical Metamaterials[END_REF] with feature sizes ranging from nano-to macroscopic scales [START_REF] Gao | The status, challenges, and future of additive manufacturing in engineering[END_REF][START_REF] Meza | Resilient 3d hierarchical architected metamaterials[END_REF]. One such example are 3D-printed truss-based architectures, made of beams or hollow tubes [START_REF] Intrigila | The compressive response of additively-manufactured hollow truss lattices: an experimental investigation[END_REF] as elementary building blocks, which have demonstrated-among others-exceptional stiffness-and strength-to-weight ratios [START_REF] Chu | Design for additive manufacturing of cellular structures[END_REF][START_REF] Jang | Fabrication and deformation of three-dimensional hollow ceramic nanostructures[END_REF][START_REF] Meza | Strong, lightweight, and recoverable three-dimensional ceramic nanolattices[END_REF][START_REF] Meza | Reexamining the mechanical property space of three-dimensional lattice architectures[END_REF][START_REF] Hahn | Rapid assembly of small materials building blocks (voxels) into large functional 3d metamaterials[END_REF]. Current applications go beyond stiff and strong structures and leverage soft base materials such as elastomers [START_REF] Jiang | Highly-stretchable 3D-architected Mechanical Metamaterials[END_REF] or hydrogels [START_REF] Li | Stretchable 3D lattice conductors[END_REF], which exhibit hyperelastic behavior under large, nonlinear deformation [START_REF] Jiang | Highly-stretchable 3D-architected Mechanical Metamaterials[END_REF][START_REF] Weeger | Digital design and nonlinear simulation for additive manufacturing of soft lattice structures[END_REF] and display a nearly fully reversible mechanical response [START_REF] Wu | Proceedings of the ASME (2020) international design engineering technical conferences and computers and information in engineering conference[END_REF]. Moreover, nonlinear constitutive behavior has been exploited through material heterogeneity (e.g., by multi-material polyjet printing [START_REF] Patpatiya | A review on polyjet 3d printing of polymers and multi-material structures[END_REF]) or heterogeneous pre-stresses, which can be triggered by external stimuli such as thermal effects [START_REF] Wagner | Large shape transforming 4D auxetic structures[END_REF][START_REF] Ding | 4D rods: 3D structures via programmable 1D composite rods[END_REF][START_REF] Kotikian | Untethered soft robotic matter with passive control of shape morphing and propulsion[END_REF] or swelling [START_REF] Jiang | Bifurcation-based embodied logic and autonomous actuation[END_REF] to induce dramatic changes in shape. This has allowed for the design of active materials and so-called 4D-printing [START_REF] Lumpe | Computational design of 4d printed shape morphing mulit-state lattice structures[END_REF][START_REF] Xia | Responsive materials architected in space and time[END_REF].

In parallel to experimental advances, efficient predictive computational tools have been developed for the design process of such structures (see e.g., the review of [START_REF] Kochmann | Multiscale modeling and optimization of the mechanics of hierarchical metamaterials[END_REF]. A particular challenge is to accurately capture the mechanical response at the beam level (combining geometric and material nonlinearities) with a moderate computational cost. This typically reduces to a competition between efficiency and accuracy.

Towards efficient simulations, homogenization approaches have been introduced to model periodic truss structures as effective continua; see, e.g. the FE 2 -type formulations of [START_REF] Vigliotti | Linear multiscale analysis and finite element validation of stretching and bending dominated lattice materials[END_REF]; [START_REF] Vigliotti | Non linear constitutive models for lattice materials[END_REF]; [START_REF] Glaesener | Continuum models for stretchingand bending-dominated periodic trusses undergoing finite deformations[END_REF][START_REF] Glaesener | Continuum representation of nonlinear three-dimensional periodic truss networks by on-the-fly homogenization[END_REF] and references therein. To gain efficiency, the beam formulation underlying those techniques is usually limited to small strains-analytical treatments at finite strains were limited to networks of springs without flexural or torsional contributions [START_REF] Triantafyllidis | On higher order gradient continuum theories in 1-d nonlinear elasticity. derivation from and comparison to the corresponding discrete models[END_REF][START_REF] Wang | Three-dimensional nonlinear open-cell foams with large deformations[END_REF][START_REF] Desmoulins | Local and nonlocal continuum modeling of inelastic periodic networks applied to stretching-dominated trusses[END_REF]. [START_REF] Glaesener | Viscoelastic truss metamaterials as time-dependent generalized continua[END_REF] recently showed that such a homogenized description of periodic truss metamaterials, using an FE 2 -type setting, also performs well for linear viscoelastic base materials when undergoing large rotations but small strains. Their comparison with experimental results showed excellent agreement up to moderate strains, while the model failed to provide accurate predictions at increasing strain levels, where both geometric and material nonlinearity become appreciable. Obviously, such approaches depend on an accurate representation at the individual strut level.

The most accurate numerical description requires fully-resolved finite element (FE) analysis, which resolves the 3D geometry of each strut. For example, [START_REF] Truszkiewicz | Mechanical behavior of 3D-printed polymeric metamaterials for lightweight applications[END_REF] and [START_REF] Wagner | Large shape transforming 4D auxetic structures[END_REF] demonstrated that such FE simulations accurately capture the mechanics of 3D-printed polymer trusses. [START_REF] Bluhm | Nonlinear compressive stability of hyperelastic 2D lattices at finite volume fractions[END_REF] characterized the effects of truss joints and general finite-strain effects in plane strain, using hyperelastic constitutive laws. Similarly, [START_REF] Mirzaali | Shape-matching soft mechanical metamaterials[END_REF] simulated the response of 2D periodic trusses in plane stress, using a neo-Hookean hyperelastic model. Full FE simulations have also been employed to model active materials such as the 3D pop-up structures of [START_REF] Yan | Mechanical assembly of complex, 3d mesostructures from releasable multilayers of advanced materials[END_REF]. While such approaches provide a high level of detail and therefore excellent accuracy, they all come with high computational costs.

Beam theories-including those of Euler, Bernoulli, and Timoshenko-provide a computationally efficient alternative to fully resolved simulations, where those are applicable. [START_REF] Ding | 4D rods: 3D structures via programmable 1D composite rods[END_REF] analyzed 4D-printed structures that are thermally activated, relying on a Cosserat beam formulation while considering linear thermoelasticity. [START_REF] Weeger | Isogeometric sizing and shape optimization of 3D beams and lattice structures at large deformations[END_REF] explored linear elastic beam theory with rigid but multi-material cross-sections. Elastoplastic truss lattices were studied by [START_REF] Tankasala | Tensile response of elastoplastic lattices at finite strain[END_REF] using Timoshenko beam theory. While there is clearly a wide range of existing beam theories and associated computational implementations, the state of the art is limited when modeling beams with complex cross-sections and base material behavior undergoing large, nonlinear deformation. Besides, the range of validity of each particular model is limited and the choice of an adequate model depends on the scenario to be simulated.

In this contribution, we present a hybrid approach, which leverages the efficiency of structural beam models while accounting for both geometric and material nonlinearity in complex cross-sections. Based on asymptotic dimension reduction, our framework may be viewed as a special type of FE 2 with a beam representation at the macroscopic scale and a fully-resolved FE description of the cross-section at the microscopic scale.

Asymptotic theory assumes a separation of scales between the beam's length and its cross-sectional dimensions. We thus assume slow variations along the beam's length in the stresses and strains, so that the structure's kinematics can be captured effectively using a formal two-scale expansion [START_REF] Berdichevskii | On the energy of an elastic rod[END_REF][START_REF] Bermudez | Une justification des équations de la thermoélasticité des poutres à section variable par des méthodes asymptotiques[END_REF][START_REF] Cimetiere | Asymptotic theory and analysis for displacements and stress distribution in nonlinear elastic straight slender rods[END_REF]. The 3D problem is thus split into (i) a series of 2D microscale problems defined on the cross-section and (ii) an effective, macroscopic representation of the beam's center-line based on an effective strain energy capturing the mechanics of the beam's center-line undergoing a combination of stretching, bending, and torsional strains.

Numerical schemes following analogous approaches have been introduced, e.g., for composite structures with intricate cross-sections [START_REF] Cesnik | VABS : A New Concept for Composite Rotor Blade Modeling[END_REF][START_REF] Yu | GEBT: A general-purpose nonlinear analysis tool for composite beams[END_REF][START_REF] Gupta | Beam theory for asymptotic analysis of aperiodic and inhomogeneous structures[END_REF], thus reducing computationally expensive 3D problems to equivalent 1D beam problems-in the small-strain regime. Similar reduction procedures were recently extended to finite elasticity (Jiang and Yu, 2016;Jiang et al., 2016), prestrains [START_REF] Moulton | Morphoelastic rods III: Differential growth and curvature generation in elastic filaments[END_REF] and higher-order (strain gradient) effects (Lestringant and Audoly, 2020;[START_REF] Audoly | A discrete, geometrically exact method for simulating nonlinear, elastic and inelastic beams[END_REF].

While closed-form, semi-analytical expressions for the effective nonlinear asymptotic energy are available for simple geometries [START_REF] Moulton | Morphoelastic rods III: Differential growth and curvature generation in elastic filaments[END_REF][START_REF] Audoly | A discrete, geometrically exact method for simulating nonlinear, elastic and inelastic beams[END_REF], there is no general analytical solution for nonlinear deformations and arbitrary cross-sections and hyperelastic constitutive laws. Jiang and Yu (2016); Jiang et al. (2016) proposed a numerical implementation of the microscopic problem on the cross-section at leading order and analyzed the nonlinear coupling between stretching and bending/twisting for varying cross-sections, yet without a macroscale representation suitable for simulating boundary value problems such as, e.g., truss networks.

Here, we introduce a general numerical framework capable of handling arbitrary crosssectional geometries, material and pre-strain distributions and able to solve complex non-linear boundary value problems on the macroscale, including large assemblies of slender beams, as typically found in truss-based architected materials.

Our implementation effectively decomposes the 3D problem of each strut into a set of 2D boundary value problems, which we solve a-priori and offline (i.e., before solving any macroscopic problem) on a reference cross-section. Tabulation of the cross-sectional solution results in an effective strain energy density, which in turn leads to a nonlinear 1D beam problem on the macroscale, which we solve using a discrete beam framework (Lestringant et al., 2020;Lestringant and Kochmann, 2020). We thus realize numerically the theoretical recipe introduced by Audoly and Lestringant (2021) at leading order. For validation, we compare predictions of our hyperelastic beam implementation with fullyresolved FE simulations for single beams under various load conditions and for 2D truss architectures. Further, we discuss the difference between a linear beam model and our hyperelastic model for 3D truss architectures. Our results highlight the importance of accounting for nonlinear effects stemming from both geometric and material nonlinearity.

The remainder of this contribution is structured as follows. In Section 2, we introduce the numerical dimension reduction procedure: we parameterize the 3D beam problem based on the beam center-line, detail the derivation of the 1D model by formal asymptotic dimension reduction, and discuss the implementation of our numerical scheme. The performance of this framework is demonstrated in Section 3 in two parts.

First, Section 3.1 illustrates the effective strain energy density for various cross-sections, including multi-material and pre-strained cross-sections. Second, we solve macroscopic boundary value problems involving assemblies of beams in Section 3.2. Finally, Section 4 concludes our study.

Numerical dimension reduction procedure

Parametrization of the 3D elasticity problem

Following the theoretical approach proposed by [START_REF] Audoly | A discrete, geometrically exact method for simulating nonlinear, elastic and inelastic beams[END_REF], we start from a generic slender 3D prismatic solid, a beam, as shown in Fig. 1. We introduce a 3D Cartesian basis {e 1 , e 2 , e 3 } with e 3 being aligned with the longitudinal direction of the beam in its reference configuration and passing through the centroids of its crosssections. We define the position of a point within the beam in the reference configurations as X = X 1 e 1 + X 2 e 2 + Se 3 , where S is the longitudinal arc-length coordinate, and (X 1 , X 2 ) capture the undeformed geometry of the cross-section, as shown in Fig. 1.

We denote by A the cross-sectional area in the reference configuration. Here and in the following, we restrict ourselves to initially straight beams (prismatic solids in 3D), although the framework can readily be extended to initially curved beams by adapting the center-line parametrization to account for a uniform or slowly varying natural curvature.

d 3 (S) d 1 (S) F = @x @X e 3 e 2 e 1 X 2 S X X 1 reference configuration y 1 y 2 y 3 current configuration x r(S) d 2 (S)
r(S)

Figure 1: Reference (undeformed) and current (deformed) configurations of a slender beam in 3D space.

The reference configuration is parameterized by coordinates (X 1 , X 2 , S), while the current configuration is described by the combination of the center-line kinematics (the center-line r(S) and the frame directors

d i (S)
) and the microscopic (cross-sectional) positions y(X 1 , X 2 , S).

In the current configuration, let the deformation mapping x = x(X 1 , X 2 , S) assign to each reference point X its deformed counterpart. We define the position of the center-line r(S) as the average of the current position x(X 1 , X 2 , S) of all points in a cross-section, i.e.,

r(S) = ⟨x(X 1 , X 2 , S)⟩ with ⟨•⟩ = 1 A A (•) dA, (1)
where dA = dX 1 dX 2 . r(S) captures the kinematics of an imaginary line, which passes through the centroids of every cross-section of the prismatic solid. At every S, we introduce a set of orthonormal material frame directors d i (S) (with i = 1, 2, 3), which define the local average orientation of the cross-section. Let d 3 (S) = r ′ (S)/|r ′ (S)| be the unit vector tangent to the center-line; here and in the following, a prime denotes a derivative with respect to S. We enforce orthonormality via

d i (S) • d j (S) = δ ij ∀S, (2)
where i, j = 1, 2, 3 and δ ij represents Kronecker's delta. The geometry of the deformed cross-section is captured by the microscopic 3D position vector y = (y 1 , y 2 , y 3 ) T , whose components indicate-for a reference point X-its position in the deformed cross-section relative to the centroid at r(S) and in the directors' frame:

y i (X 1 , X 2 , S) = [x(X 1 , X 2 , S) -r(S)] • d i (S). (3)
With this choice of parametrization, the position vector in the current configuration writes as follows:

x

(X 1 , X 2 , S) = r(S) + y i (X 1 , X 2 , S)d i (S), (4)
where Einstein's summation convention on repeated Latin indices is implied (here and in the following). The uniqueness of this parametrization is enforced by the following set of kinematic constraints:

q(y) = ⟨y 1 ⟩, ⟨y 2 ⟩, ⟨y 3 ⟩, ⟨X 1 y 2 -X 2 y 1 ⟩ T = 0. (5)
These constraints ensure that the center-line runs through the centroids of the crosssections, and it guarantees uniqueness of the set of orthonormal vectors d 1 and d 2 in the plane orthogonal to d 3 .

Through this parametrization, the deformation depends explicitly on the macroscopic (i.e., center-line) translations and rotations through r(S) and d i (S) on one hand, and on the microscopic (i.e., cross-sectional) positions y(X 1 , X 2 , S) on the other hand. This allows us to define the apparent macroscopic strains based on the standard analysis of framed curves ( see e.g. [START_REF] Stoker | Differential Geometry[END_REF]), which are defined by the center-line and material frame directors. We start by defining the axial strain ε(S) via

r ′ (S) = [1 + ε(S)] d 3 (S). ( 6
)
The macroscopic bending and torsional strains are captured by the vector

κ(S) = κ 1 (S)d 1 (S) + κ 2 (S)d 2 (S) + τ (S)d 3 (S), (7)
which satisfies

d ′ i (S) = κ(S) × d i (S). ( 8
)
Here, κ i (S) captures bending about d i (S) (for i = 1, 2), and τ (S) the twist about d 3 (S).

These macroscopic strains1 are collected in the vector

E(S) = ε(S), κ 1 (S), κ 2 (S), τ (S) T . (9)
The deformation gradient follows from our parametrization as

F = ∂x(X 1 , X 2 , S) ∂X α ⊗ e α + x ′ (X 1 , X 2 , S) ⊗ e 3 . (10)
Here and in the following, Greek indices are in the 2D plane of the cross-section, so the above summation is over α = 1, 2. The deformation gradient can be expanded as

F = ∂y i (X 1 , X 2 , S) ∂X α d i (S) ⊗ e α + 1 + ε(S) δ i3 + η ijk κ j (S)y k (X 1 , X 2 , S) + y ′ i (X 1 , X 2 , S) d i (S) ⊗ e 3 , (11)
where η ijk is the Levi-Civita symbol. The right Cauchy-Green tensor, defined in the reference configuration, is given by

C = F T F or C E(S), y(X 1 , X 2 , S), y ′ (X 1 , X 2 , S) = ∂y i (X 1 , X 2 , S) ∂X α ∂y i (X 1 , X 2 , S) ∂X β e α ⊗ e β + a i ∂y i (X 1 , X 2 , S) ∂X α [e 1 ⊗ e 3 + e 3 ⊗ e 1 ] + a 2 i e 3 ⊗ e 3 (12)
with

a i = 1 + ε(S) δ i3 + η ijk κ j (S)y k (X 1 , X 2 , S) + y ′ i (X 1 , X 2 , S). ( 13
)
Note that, aside from the macroscopic strains E(S), C depends on both the microscopic field y and its derivative y ′ along the longitudinal axis.

The above allows us to formulate the total 3D strain energy functional for a hyperelastic rod of undeformed length L as

Φ [E, y] = L 0 A W C E(S), y(X 1 , X 2 , S), y ′ (X 1 , X 2 , S) dA dS, ( 14
)
where W (C) defines a hyperelastic strain energy density, which is assumed to be polyconvex.

Although our examples in Section 3-for ease of implementation-use a compressible hyperelastic strain energy density, incompressibility can readily be enforced by introducing an additional term with a Lagrange multiplier in the energy. The incompressibility term is then treated in the same fashion as the constraints in (5). No restrictions are imposed on the cross-sectional geometry; the cross-section may be inhomogeneous in terms of the material distribution, and it may be pre-stressed (which will be discussed in detail in Section 2.4).

The total potential energy functional follows as

Ψ[E, y] = Φ [E, y] + L 0 V r(S), d i (S) dS, (15)
where V (r(S), d i (S)) is a density of external potential, which we assume to vary slowly along S and to represent external forces. At equilibrium the total potential energy (15) is stationary with respect to (E, y), subject to the constraint (5). With this, we have fully described the 3D (hyper-)elasticity problem based on the center-line parametrization.

Note that in the potential (15) external forces act only on the macroscopic displacements and rotations. For slender structures this is a legitimate restriction. As a consequence, the addition of the external potential will not affect the dimension reduction procedure described below, and we will therefore focus on the potential (14) in the next section.

Effective energy from a formal asymptotic dimension reduction

To introduce an effective beam model depending on the macroscopic strain fields only, we aim at condensing out the microscopic unknowns in (14). Our approach is variational and formulated as a stationarity problem for the energy functional (14) under the constraint (5), similarly to Jiang and Yu (2016) and [START_REF] Audoly | A discrete, geometrically exact method for simulating nonlinear, elastic and inelastic beams[END_REF].

We formulate the associated Lagrangian L as

L = Φ [E, y] - L 0 f (S) • q y) dS, ( 16
)
where the Lagrange multiplier f (S) = (f 1 (S), f 2 (S), f 3 (S), f 4 (S)) T is a four-dimensional vector conjugate to the constraints q(y) in (5). To find y ⋆ and f ⋆ , which render ( 16) stationary for a given distribution of macroscopic strains E(S), we require the first variation of ( 16) for admissible perturbations δy and δf of the microscopic position field y and the Lagrange multiplier f , respectively, to vanish. Thus, for a prescribed distribution of macroscopic strains, we find y ⋆ and f ⋆ from solving

     L 0 A ∂W (C(E,y)) ∂y
• δy dAf • q(δy) dS = 0 for any admissible perturbation δy,

q(y) = 0, (17)
where we used the following notation for the directional derivative of a function f depending on a vector field u:

∂f (u) ∂u • δu = lim t→0 f (u + tδu) -f (u) t . (18)
Replacing y ⋆ in the expression of the energy in ( 14) leaves us with the condensed macroscopic potential energy functional

Φ ⋆ [E] = Φ [E, y ⋆ ] , (19)
which only depends on the macroscopic strain field E(S) and no longer on the microscopic unknowns. Ideally, we aim to express the above as

Φ ⋆ [E] = L 0 Ŵ E(S) dS, ( 20
)
where Ŵ is an effective strain energy density that depends only on the macroscopic strains E (and, of course, on the cross-sectional geometry and base material's constitutive behavior). However, the integrals in ( 16) are nonlocal in y(X 1 , X 2 , S), so that the stationarity condition must be satisfied over the entire domain S ∈ [0, L].

To find a local solution of the stationarity conditions in (17), we approximate y ⋆ via an asymptotic expansion; see, e.g., [START_REF] Berdichevskii | On the energy of an elastic rod[END_REF]; [START_REF] Bermudez | Une justification des équations de la thermoélasticité des poutres à section variable par des méthodes asymptotiques[END_REF]; [START_REF] Cimetiere | Asymptotic theory and analysis for displacements and stress distribution in nonlinear elastic straight slender rods[END_REF]; [START_REF] Hodges | Nonlinear Composite Beam Theory[END_REF]; Lestringant and Audoly (2020). The key assumption behind this approximation is that the macroscopic strains vary on a length scale significantly larger than that of the cross-sectional dimensions, so that the axial gradients scale as

E ′ = O(ϵ), y ′ = O(ϵ), E ′′ = O(ϵ 2 ), y ′′ = O(ϵ 2 ), (21)
where ϵ ≪ 1 is a small parameter, typically on the order of the beam slenderness ratio (i.e., the ratio of the characteristic dimensions of the cross-section w to the beam length L). With this scaling assumption, we seek the solution y ⋆ to the problem in ( 17) of the form

y ⋆ (X 1 , X 2 , S) = ỹ⋆ (X 1 , X 2 , S) + y ⋆ 1 (X 1 , X 2 , S) + y ⋆ 2 (X 1 , X 2 , S) + • • • with ỹ⋆ (X 1 , X 2 , S) = O(1), y ⋆ i (X 1 , X 2 , S) = O(ϵ i ). ( 22
)
Inserting ( 22) and ( 21) into the stationarity condition ( 17) allows us to solve the latter order by order in ϵ. At leading order, all longitudinal gradients vanish, and we can express the deformation gradient (11) as

F = F + O(ϵ) with F E(S), ỹ(X 1 , X 2 , S) = ∂ ỹi (X 1 , X 2 , S) ∂X α d i (S) ⊗ e α + [(1 + ε(S))δ i3 + ϵ ijk κ j (S)ỹ k (X 1 , X 2 , S)] d i (S) ⊗ e 3 . (23)
The right Cauchy-Green strain tensor at leading order follows as

C = C + O(ϵ) with C = F T F , which-owing to (23)-is local in ỹ(X 1 , X 2 , S)
, unlike the exact Cauchy-Green tensor in ( 12). The leading-order total potential energy functional writes

Φ = Φ + O(ϵ) with Φ [E, ỹ] = L 0 A W C E(S), ỹ(X 1 , X 2 , S) dA dS. ( 24
)
The stationarity condition (17) at leading order follows as

     L 0 A ∂W ( C(E,ỹ)) ∂ ỹ • δ ỹ dA -f • q(δ ỹ) dS = 0 for any admissible perturbation δ ỹ, q(ỹ) = 0. (25)
Equations ( 23) and ( 24) confirm that the nonlinear problem (25) does not involve derivatives ỹ′ = ∂ ỹ/∂S: the microscopic positions ỹ at leading order can hence be obtained by solving a local nonlinear elasticity problem on a reference cross-section for a fixed set E of macroscopic strain values. This writes in variational form as

     A ∂W ( C(E,ỹ)) ∂ ỹ • δ ỹ dA -f • q(δ ỹ) = 0 for any admissible perturbation δ ỹ, q ỹ) = 0. (26)
Let us denote by (ỹ ⋆ E (X 1 , X 2 ), f ⋆ E ) the solution of the problem in (26) for a given E. We obtain the leading-order approximation for the effective strain energy density W eff in (20) by replacing the microscopic positions ỹ(S) for every coordinate S by the leadingorder solution ỹ⋆ E(S) (X 1 , X 2 ). The energy to leading order thus becomes

Φ⋆ [E] = L 0 W eff E(S) dS with W eff (E(S)) = A W C E(S), ỹ⋆ E(S) (X 1 , X 2 ) dA. ( 27
)
This reduces the total energy functional (15) to

Ψ⋆ [E] = L 0 W eff E(S) dS = Φ⋆ [E] + L 0 V r(S), d i (S) dS. ( 28
)
This dimension reduction procedure requires that, for every load/time step of the structural macroscopic problem defined by the energy functional (28), we solve a series of microscale, local cross-sectional problems ( 26), as summarized in Fig 2.

When performing the asymptotic expansion up to the second order in ϵ, the strain energy ( 14) can be expanded as (Audoly and Lestringant, 2021)

Φ ⋆ [E] = L 0 W eff E(S) dS + 1 2 L 0 E ′ (S)B E(S) E ′T (S) dS + O(ϵ 3 )L (29)
plus boundary terms that are omitted here for concision. In (29), W eff (E(S)) is the strain energy density at leading order, as defined in ( 27), and B (E) is a second-order tensor of gradient moduli, which depends nonlinearly on the macroscopic strains E.

Such higher-order terms are typically required to capture localization effects along the length of the beam, for which the second term in (29) provides regularization. Note that computing the components of B (E(S)) requires solving for the microscopic unknowns at order ϵ (i.e., solving for y ⋆ 1 (E, E ′ ) in ( 22)), which is computationally expensive. As localization plays only a negligible role in our applications, we expect the leading-order strain energy (27) (the first term in ( 29)) to be sufficient to accurately capture the mechanics of truss-based architected materials-this assumption will be confirmed by our numerical results in Section 3. Therefore, in the following we focus on the leadingorder energy contribution, which implies reduced computational costs compared to the full expression in (29). Indeed, computing the matrix B (E(S)) requires also solving for the microscopic unknowns at order ϵ (i.e., solving for y ⋆ 1 (E, E ′ ) in ( 22)).

) 27)

W e↵ (E) = Z A W ( C(E, ỹ? E )) dA (
? [E] = Z L 0 W e↵ (E(S)) dS + Z L 0 V (r(S), d i (S))dS

Numerical implementation of the two-scale approach

By the above dimension reduction method, solving a macroscopic (initial) boundary value problem ((I)BVP)-making use of the effective strain energy potential (27)involves solving a series of local 2D problems (26) on the cross-section for each local combination of strains ε(S), κ 1 (S), κ 2 (S), and τ (S) at each S. These 2D elasticity problems lend themselves naturally to be solved using conventional FE with a discretization based on 2D solid elements. This can be done either on the fly at every load/time step and for every macroscale coordinate S (in an FE 2 fashion), or alternatively computed a priori for a catalog of macroscopic strain combinations for a given cross-sectional geometry and hyperelastic constitutive law(s) of the base material(s) within the cross-section.

We here take the route of tabulation in a pre-processing step: at the cost of an initial (offline) effort, we save computational costs when solving the macroscopic (I)BVP (the solutions of the microscopic problems are retrieved from look-up tables). More refined techniques such as kriging or machine learning could be explored to improve the efficiency of the approach.

tabulate W e↵ on a discrete grid ("

1 ...." N ) ⇥ ( 1 1 .... N 1 ) ⇥ ( 1 2 .... N 2 ) ⇥ (⌧ 1 ....⌧ N ) interpolate W e↵ ( Ēijkl ) ) W e↵ (E) and di↵erentiate to compute ⌃ i (37) and K i (38) W e↵ ( Ēijkl ) = Z A W ( C( Ēijkl , ỹ? Ēijkl )) dA
pass the discrete strain state E i (31) to the material model:

for given macroscale Dirichlet BC's & external loads F external :

solve the 2D BVP ( 26) on the cross-section for ỹ?

E ijkl material model
discretize the cross-section (using Gmsh)

2D nonlinear FEA on the cross-section (FEniCS ( 2022)) effective model on the macroscale is discretized using a discrete rod framework (Lestringant et al., 2020;Lestringant and Kochmann, 2020) and uses the a-priori obtained effective strain energy density to compute the effective stress and stiffness matrices for particular combinations of the macroscopic strains

x i 1 x i x i+1 ! i 1 ! i d i 1 3 d i
E i in each element.
Fig. 3 summarizes our numerical two-scale procedure. Given a particular crosssectional geometry and hyperelastic strain energy density W (C) (or several strain energy densities for heterogeneous cross-sections), we solve the microscopic nonlinear elasticity problem (26) on a 2D cross-section for a range of macroscopic strain combinations E(S). To that end, we discretize the local problem: we mesh the 2D domain using Gmsh [START_REF] Geuzaine | Gmsh: A 3-d finite element mesh generator with built-in pre-and post-processing facilities[END_REF]) and use quadratic triangular elements with full integration within FEniCS (FEniCS ( 2022)). The macroscopic strains E(S) are passed to the microscopic problem via the deformation gradient F (see ( 23)) which enters into the definition of the energy in ( 24) via the right Cauchy-Green tensor C, while the kinematic constraints on the cross-section ( 5) are enforced via Lagrange multipliers f . No additional boundary conditions are applied. The nonlinear variational problem ( 26) is solved monolithically for the Lagrange multipliers f ⋆ E (and p ⋆ E , if incompressible materials are considered) and for the microscopic unknowns ỹ⋆ E , using a Newton-Raphson scheme. For any macroscopic strain combination E, we use the solution to compute the effective strain energy density by numerical integration over the cross-section using the definition 27). The python-based FEniCS code used for solving a microscopic problem on the cross-section is provided in the Supplementary Material.

W eff (E) = A W C E, ỹ⋆ E (X 1 , X 2 ) dA (
We repeat this procedure for a catalog of N ×N ×N ×N different macroscopic strains combinations Ēijkl = (ε i , κj 1 , κk 2 , τ l ) and store the resulting effective strain energy density values W ijkl eff = W eff (ε i , κj 1 , κk 2 , τ l ) with i, j, k, l ∈ {1, • • • , N } in a four-dimensional array. This tabulated effective energy density is used as the material model in the macroscale problem. As the macroscopic problem will require derivatives of the effective energy density, we transform the tabulated W ijkl eff -values into an approximate albeit continuous and differentiable energy density W eff (E) by a four-dimensional cubic Hermite interpolation of the tabulated values [START_REF] De Boor | High accuracy geometric hermite interpolation[END_REF]Wolfram Research, 2008).

On the macroscale we use a discrete beam formulation based on discrete elastic rods to solve the structural (initial) boundary value problem ((I)BVP), following the implementation of Lestringant et al. (2020); Lestringant and Kochmann (2020), and using the open-source C++ library UtoBeams (2020). In the following, we only provide a brief summary of the salient implementation features, and we refer the reader to the above references for further details.

The center-line parametrization of Section 2.1 is discretized into a set of vertices connected by segments. We identify an individual element as a subset of three neighboring vertices connected by two segments, as shown in Fig. 3. Following Section 2.1, a set of orthonormal frame directors d i is defined at each segment, with d 3 defined as the unit tangent to the segment's center-line and the two other directors spanning the cross-section (in an average sense as the cross-sections are not constrained to remain planar). The full 3D kinematics of the discrete framed curve centre-line is parameterized by a reduced set of degrees of freedom, viz. the current vertex positions x i and the spin angles ω i capturing the center-line rotation about d 3 for each segment (note that we use sub-and superscripts to differentiate vertex and segment variables, respectively). We denote by N the total number of degrees of freedom and by n the number of elements in the mesh (note that there is no general relationship between N and n, since this depends on the connectivity of the mesh). For a single beam, we define element i as the one connecting vertices i -1, i, and i + 1 and segments i -1 and i. Its eleven degrees of freedom are collected in the vector

u i = (x i-1 , x i , x i+1 , ω i-1 , ω i ) T . (30)
Note that with this definition beam elements overlap (i.e., each segment is part of two beam elements). We use an updated-Lagrangian implementation, in which the current configuration serves as the reference configuration for the next discrete load/time step, such as to accurately capture large rotations through the current positions x. The deformation of discrete beam element i is characterized by the discrete strain measures

E i = (ε i , κ 1,i , κ 2,i , τ i ) T , (31)
where (31), ε i is the average axial strain of the two segments, defined as

ε i = 1 2 |x i -x i-1 | L i,i-1 -1 + |x i+1 -x i | L i+1,i -1 , (32)
where L i,j denotes the undeformed length of the segment connecting vertices i and j.

To ensure convergence of the discrete to the continuous strain measures in the limit of zero segment length, discrete bending and twisting strains are rescaled by the vertexbased undeformed Vornoi length li =

Li,i-1+Li+1,i 2

. The discrete bending strains are defined at the element's central vertex i and correspond to the projection of the curvature binormal k i on the directors d 1 and d 2 , averaged over the two neighboring segments:

κ α,i = k i li • d i-1 α + d i α 2 for α ∈ {1, 2}, (33)
where

(d k 1 , d k 2 , d k
3 ) denote the discrete set of directors at segment k, and the vertex-based discrete curvature binormal k i is defined as

k i = 2 d i-1 3 × d i 3 1 + d i-1 3 • d i 3 . ( 34
)
Finally, the discrete twist at vertex i is given by

τ i = (τ i ) ⋆ + 1 li ω i -ω i-1 + γ i , (35)
where (τ i ) ⋆ is the twist in the reference configuration (the configuration at the beginning of the load/time step in our updated-Lagrangian scheme), and γ i is a quantity called the holonomy in differential geometry, which depends on the rotations of the two neighboring segments. It captures the coupling between the torsional and flexural deformation modes.

For a detailed derivation of the above and a method for computing the holonomy as the sum of the signed area of two spherical triangles, we refer to Lestringant et al. (2020).

Rigid beam junctions are realized by coupling the 3D rotations of the segments connected at a junction through the introduction of so-called ghost elements, see Lestringant and Kochmann (2020). Besides, our formulation decouples the kinematics from the constitutive behaviour, thus allowing to leverage the results from our microscopic calculations. While both quasistatic and dynamic problems are solved in UtoBeams, we focus on quasistatic problems in the following to demonstrate our numerical two-scale approach.

Rendering the total potential energy stationary is equivalent to solving the quasistatic equilibrium equation, which for the discrete beam problem can be written as

n i=1 Σ i ∂E i ∂u i L i -F external = 0, ( 36
)
where n is the total number of elements, assembly matrix L i ∈ R 11×N maps local to global degrees of freedom via u i = L i U , where U = {x 1 , . . . , x v , ω 1 , . . . ω s } is the vector of global degrees of freedom, where v and s are, respectively, the total number of vertices and segments in the beam mesh. F external ∈ R N is the global external force vector.

Furthermore,

Σ i = li Σ(E = E i ) with Σ = ∂W eff ∂E = ∂W eff ∂ε , ∂W eff ∂κ1 , ∂W eff ∂κ2 , ∂W eff ∂τ T ( 37
)
is the integrated stress vector conjugate to the four macroscopic strain measures (evaluated with the local strain values E i ). Finally, matrix ∂Ei ∂ui ∈ R 4×11 encodes the beam kinematics and depends only on the discrete beam's geometry; and is detailed in Lestringant et al. (2020).

The generally nonlinear problem ( 36) is solved via a Newton-Raphson scheme, which makes use of the discrete stress Σ i and its derivative, the symmetric consistent tangent matrix

K i = li K(E = E i ) with K = ∂ 2 W eff ∂E∂E =           ∂ 2 W eff ∂ε∂ε ∂ 2 W eff ∂ε∂κ1 ∂ 2 W eff ∂ε∂κ2 ∂ 2 W eff ∂ε∂τ ∂ 2 W eff ∂κ1∂κ1 ∂ 2 W eff ∂κ1∂κ2 ∂ 2 W eff ∂κ1∂τ ∂ 2 W eff ∂κ2∂κ2 ∂ 2 W eff ∂κ2∂τ
sym.

∂ 2 W eff ∂τ ∂τ           . (38)
Both Σ and K can be obtained from the effective strain energy density W eff (E), which here serves as the material model and which was obtained a priori from interpolating the tabulated energy values from the microscale cross-sectional problem.

Modeling pre-strain and growth

Pre-strain within beams is often used as an actuation mechanism (e.g., in soft robotics applications) as well as in naturally growing fiber networks. The effects of finite prestrain can be included in our two-scale formulation by adapting the definition of the microscale problem as follows (keeping the macroscale problem unchanged). We restrict ourselves to pre-strain distributions that are invariant or slowly varying along the length of a beam, so that the scale-separation hypotheses underlying the asymptotic dimension reduction in ( 21) and ( 22) are satisfied. Pre-strain is incorporated via a multiplicative decomposition of the deformation gradient F ; see, e.g., [START_REF] Lestringant | Elastic rods with incompatible strain: Macroscopic versus microscopic buckling[END_REF] in the context mechanical pre-strain and Rodriguez et al. structure and G is a prescribed growth tensor. In compliance with the leading-order expansion, we require G = G(X 1 , X 2 ) to be invariant along the length of the beam in subsequent examples, although a slowly varying G would be considered locally invariant at leading order and hence be admissible within the current framework. Note that this does not exclude G to vary within the cross-section. The strain energy density (24) now depends on the elastic contribution A = F G -1 through C = A T A:

Φ[E, ỹ] = L 0 A W C(A(E, ỹ)) det G dA dS ( 39
)
with dA = dX 1 dX 2 in the (pre-grown) reference configuration.

The 2D microscopic problem (26) can thus be reformulated to account for the distribution G(X 1 , X 2 ) by replacing the expression of the strain energy ( 24) by ( 39). The numerical procedure described in Section 2.3 can be readily applied without any changes to obtain a tabulated material model (dependent on the choice of G(X 1 , X 2 )), which integrates the effects of the growth distribution. This will be demonstrated in the example in Section 3.1.3.

Note that, when the macroscopic strains vanish (i.e., when A = 0), the beam is in general not stress-free for a non-zero growth tensor G(X 1 , X 2 ). For particular choices of the growth distribution, the beam will relax its energy in a curved configuration (with zero macroscopic axial and torsional strains), akin to a beam with natural curvature, as will be shown in an example in Section 3.1.3.

Numerical examples

We demonstrate the efficiency and accuracy of our approach by a suite of examples capturing nonlinear effects in slender hyperelastic structures. In Section 3.1, we present effective material models obtained by solving a series of microscopic problems for a range of cross-sections, including multi-material cross-sections and cross-sections experiencing growth. In Section 3.2, we show results from the full two-scale procedure, including solutions of macroscopic boundary value problems for individual beams and for beam assemblies, relevant for the analysis of hyperelastic truss-based architected materials. In all examples, the constitutive behavior of the base material(s) is described as an isotropic, compressible neo-Hookean solid with strain energy density

W (C) = µ 2 (tr C -3) -µ ln √ det C + λ 2 ln √ det C 2 (40)
with shear modulus µ = E/[2(1 + ν)] and Lamé parameter λ = Eν/[(1

+ ν)(1 -2ν)],
where E is the Young modulus and ν Poisson's ratio. In the following examples, we consider a Poisson's ratio of ν = 0.3 and normalize results by Young's modulus E. Note that the presented approach is general and applies to any (polyconvex) strain energy density, so the above choice of the hyperelastic constitutive model is arbitrary.

Effective material models

We first showcase results for the effective strain energy density W eff (E), obtained from tabulating solutions from the nonlinear microscale problem (26) on a variety of homogeneous and heterogeneous (multi-material) cross-sections, first without and later with growth.

Homogeneous (single-material) cross-sections

Fig. 4 illustrates six different cross-sectional geometries, whose effective energy density W eff (E) we compute. All six cross-sections fit into a square of (normalized) side length w = 1, and we use triangular finite elements of approximate size 0.025w in the 2D meshes.

Fig. 4 also shows the effective energy density for the isolated modes of stretching, bending, and twisting. In the limit of small strains, we expect the structural response of the six structures to be identical to that of a linear elastic Euler-Bernoulli beam, viz.

W EB (E) = EA 2 ε 2 + EI 1 2 κ 2 1 + EI 2 2 κ 2 2 + µJ 2 τ 2 , ( 41
)
where A is the cross-sectional area, I 1 and I 2 the area moments of inertia about the two shown principal axes, and J is the torsion constant-all dependent on the geometry and all computed in the undeformed configuration. For comparability, we rescale W eff in Fig. 4 by EA for stretching, by EI 1 and EI 2 for bending about d 1 and d 2 , respectively, and by µJ for torsion. The Euler-Bernoulli energy density is included in the plots of Fig. 4, showing excellent agreement in the small-strain limit in all cases.

Also shown are undeformed (dotted) and deformed cross-sections as well as a color map of the (normalized) microscopic displacements (showing the dominant displacement Emax | w component for each load case at the maximum applied strain). Pure axial stretching depends only on the cross-sectional area, so that the normalized energy vs. strain curves of all cross-sectional geometries coincide (Fig. 4a). Cross-sections exhibit a similar Poisson effect of in-plane shrinking with increasing macroscopic stretch. By contrast, for bending and torsion, the cross-sectional geometry has a fundamental influence on the effective nonlinear macroscopic response (Fig. 4b-d). The tubular cross-sections in bending undergo significant flattening, which is a well-known nonlinear phenomenon, the Brazier effect [START_REF] Brazier | On the flexure of thin cylindrical shells and other" thin" sections[END_REF][START_REF] Coman | Bifurcation instabilities in finite bending of circular cylindrical shells[END_REF]. This effect is captured here: the thin tube crosssection shows a more pronounced deformation, and flattening starts at an earlier stage compared to the thicker tube. Note that we abort simulations before self-contact within the cross-sections occurs, which is why some of the curves end prematurely.

W e↵ (E) µJ W e↵ (E) EI 2 W e↵ (E) EI 1 W e↵ (E) EA
[⌧ w] 2 [ 1 w] 2 [ 2 w] 2 E = (0,  1 , 0, 0) E = (0, 0,  2 , 0) (d) (c) " 2 (
All results agree well with the Euler-Bernoulli model at small strains. The solid crosssections by contrast (i.e., the solid square and circle) do not considerably deviate from the Euler-Bernoulli response over the strain range considered here. The Poisson effect due to bending is noticeable in all cases, and warping of the cross-section is observed for the noncircular cross-sections. Except for the two solid cross-sections, all cases show a significant deviation from the linear theory at moderate to large strains, which demonstrates the importance of considering both a nonlinear constitutive law for the base material and geometric nonlinearity at larger strains.

Heterogeneous (multi-material) cross-sections

As an example of a heterogeneous cross-section, we analyze a two-phase composite cross-section, which embeds an auxetic re-entrant honeycomb/butterfly inclusion into a square matrix of side length w, as shown in Fig. 5a. The stiff butterfly inclusion has a Young's modulus ten times that of the surrounding matrix. Fig. 5 illustrates the deformed cross-section due to bending with κ 1 w about d 1 , highlighting the distribution of displacements and stresses in (b) and (c), respectively. The average distribution of the von Mises stress σ vM = 3 2 σ dev : σ dev , where σ dev is the deviatoric part of the Cauchy stress tensor, is shown in Fig. 5c for the deformed state. The scaled energy density again agrees with the Euler-Bernoulli model at small strains but shows significant deviation with increasing strain (Fig. 5d).

W e↵ (E)

EI 1 0 0.5 1 1.5 2 2.5 3 3.5 4 0 0.5

1 1.5 2 [ 1 w] 2 d 3 d 2 d 1  2  1 ⌧ (a) (b) (c) (

Effect of a pre-strain distribution

Spatially varying pre-strain or growth is another form of heterogeneity within a crosssection, used for example in [START_REF] Ding | 4D rods: 3D structures via programmable 1D composite rods[END_REF] to manufacture an active bilayer structure made of a pre-strained elastomer attached to a glassy polymer base: when activated by heat, the glassy polymer surpasses its glass-transition temperature, softens, and thus allows the bilayer to bend and relieve the stresses generated by the incompatible prestrain. Our model captures such behavior, as we illustrate by the following two examples based on growth, inspired by [START_REF] Moulton | Morphoelastic rods III: Differential growth and curvature generation in elastic filaments[END_REF].

We first consider a square cross-section made out of a laminate of two equally sized, perfectly bonded homogeneous layers stacked on top of each other vertically (Fig. 6a).

The top layer experiences a constant axial growth of 20%, while the bottom layer stays in its natural configuration. Following Section 2.4, the deformation gradient in the laminate cross-section is F = AG with

G(X) =       1 0 0 0 1 0 0 0 1 + 0.2H(X 2 )       (42)
with the Heaviside function H(•) being non-zero in the grown top layer. The resulting effective strain energy density is shown in Fig. 6a as a function of the scaled bending strain κ 1 w. The layered cross-section reaches an energy minimum when bent by κ 1 w = -0.24 when the bending deformation compensates for the growth. Note, however, that it does not reach a state of zero energy. As shown in the insets, the stress relaxes inside the cross-section to an inhomogeneous distribution due to the incompatibility at the interface between the two layers.

As a second example, we consider a solid circular cross-section that experiences a linearly varying growth throughout the cross-section with

G(X) =       1 0 0 0 1 0 0 0 1 + 0.2 X 1       . (43)
As shown in the upper right inset in Fig. 6b, the straight configuration (i.e., when the macroscopic curvature κ 2 = 0) has a non-zero cross-sectional stress. At a macroscopic bending strain of κ 2 w = -0.2, a zero cross-sectional stress and a vanishing effective strain energy density are observed.

Nonlinear coupling between macroscopic modes of deformation

Since our effective material model derives from 3D (hyper-)elasticity, it naturally includes the nonlinear coupling between the different macroscopic deformation modes, which are not described by classical beam theory. We illustrate some of these coupling effects here for the example of a laminate of two equally sized, perfectly bonded homogeneous layers stacked on top of each other vertically (Young's modulus of the bottom layer isE bottom = 2E top ).

Fig. 7 shows the components ∂W eff (E) ∂ε and ∂W eff (E)

∂κ1
evaluated for E = (ε, 0, 0, 0) and E = (0, κ 1 , 0, 0). When pure stretching is considered, i.e., E = (ε, 0, 0, 0), the axial component ∂W eff ∂ε ≈ EA ε matches the Euler-Bernoulli model at small strains, before softening as the axial strain increases. Due to the heterogeneous moduli in the laminate we observe a heterogeneous stress distribution across the beam's cross-section, as show in Fig. 7a. This induces a curvature in the beam-similar to the effect under a heterogeneous strain distribution in a homogeneous beam, discussed in Section 3.1.3. As the bending strains are zero, a reaction moment is observed.

As may be expected for a slender beam, the moment under bending coincides with the Euler-Bernoulli prediction up to moderate strains. However, a major difference from Euler-Bernoulli theory is the observed coupling between the deformation modes. Fig. 7b shows that a macroscopic bending strain causes an axial contraction of the beam (i.e., negative values of ∂W eff (E) ∂ε when evaluated at zero axial strain ε = 0).

To further illustrate the subtle nonlinear coupling between axial extension and torsion, we consider a soft cross-section with stiff inclusions (E inclusion = 1000E matrix ). As the stiff inclusions are nearly in-extensible, their presence in the matrix produces a natural contraction when the structure is subjected to torsion, thus displaying the negative

@W e↵ (E) @" 1 EA @W e↵ (E) @ 1 1 EI 1 Cauchy stress 33 E top
Poynting effect [START_REF] Billington | The poynting effect[END_REF]Jiang and Yu, 2016) in Fig. 8a. As the components are evaluated at zero axial strain ε = 0, this results in ∂W eff (0,0,0,τ ) ∂ε > 0 (i.e., the structure in tension). For reference, an inclusion-free beam exhibiting the positive Poynting effect is shown in Fig. 8b.

Macroscopic boundary value problems

Having discussed the reduced, effective hyperelastic energy density and the competing deformation mechanisms within the cross-section, we now move on to solving macroscopic boundary value problems involving individual beams and beam assemblies. The effective hyperelastic energy density serves as the material model in the discrete beam framework, which we use to solve the boundary value problems in the following.

@W e↵ (0, 0, 0, ⌧ ) @"

1 AE matrix @W e↵ (0, 0, 0, ⌧ ) @⌧ 1 Jµ matrix 0 0.2 0.4 0.6 0.8 1 -0.

Individual beams

We first consider a beam of length L with a square cross-section of side-length w = L/50. For stretching and torsion the results are shown in Fig. 9a and Fig. 9b, respectively, where we compare the predictions of our reduced-order framework with those of fullyresolved 3D FE analysis conducted in FEniCS ( 2022). For such a comparison, caution is required with regards to the boundary conditions imposed at the two ends of the beam.

We apply homogeneous boundary conditions to prevent any boundary layer effects, which we do not expect to capture accurately with the reduced model without gradient effects.

In the case of pure stretch in Fig. 9a, the beam is displaced axially at one end while the other end is held fix along the longitudinal direction and left free to contract in-plane.

In the torsion example in Fig. 9b, a rotation (twisting angle) is applied to one end of the beam. At the other end, displacements in the plane of the cross-section are suppressed, while the longitudinal displacement is set to zero on average; i.e., warping is permitted but the average axial displacement is constrained to vanish. The resulting reaction force F and (b) the reaction moment M T in Fig. 9a and Fig. 9b, respectively, show excellent agreement with the fully resolved FE results. A comparison at the microscale level of the cross-section also yields convincing agreement: Fig. 9c compares the stresses within the cross-section as obtained from the microscale problem in the two-scale setup and a slice of the 3D FE results for the twist example. The difference in the maximum stress is about 3%. The uniform stress within the cross-section for the stretching example is identical for both the microscale and 3D FE analysis.

However, a clear difference is observed when comparing the computation times of the two approaches in Table 1 for the torsion example. We divide the computational cost into the a-priori step of the microscale simulations (which are performed once for a given cross-section) and the macroscale simulations. The time spent on the microscale simulations for the energy tabulation is given for three different mesh sizes and two different grid sizes (i.e., step sizes of the strain values for tabulation). As expected, times increase significantly with finer mesh and grid sizes, and can be quite substantial. Yet, these are offline computation times, i.e., once completed, the microscale results can be reused for arbitrarily many beams and for arbitrarily many simulations of beams having the same cross-section. With our current implementation, we find the coarser mesh ( with average element size 0.1w) and the finer grid (with grid size 10 -2 )to accurately capture the 3D behaviour, as shown in Fig. 9. Note that the computation time of this pre-processing step scales with the tabulation grid size to the fourth power. Therefore, techniques such as kriging, which allow for coarser grids, could considerably reduce the time spent in this step.

Although the pre-processing step is expensive, most importantly, the subsequent macroscale beam simulation is significantly more efficient than the fully-resolved 3D FE analysis. However, caution is required when comparing computation times, since different implementations, different programming languages and code architectures are being used (for a fair comparison, both simulations did not exploit parallelism). Note that in our current implementation calls to the interpolated energy have not been optimized (they consume in total about 7 s per solver time step; we do not count this additional time in Table 1).

Table 1: Computational time comparison for a beam with a square cross-section subjected to twisting.

The 2D BVP on the cross-section is solved on a grid with three different grid sizes for E = (ε, κ 1 , κ 2 , τ ) as indicated below, where each strain measure varies from 0 -0.1. Two different 2D mesh resolutions (average element size) are given below. The solver time for one step in the macroscopic torsion problem, where a cantilever beam's free end is rotated up to θ = 4 rad in increments of ∆θ = 0.1 rad is listed below. The macroscopic beam converges in three and the 3D FE model in five iterations. The simulation times are given for a simulation run on a single core. * Note that the beam simulation time excludes the interpolation time (ca. 7 s per load step ). To further illustrate the nonlinear effects captured by our approach, we next consider a hollow tube with a thin-walled circular cross-section of diameter w and wall thickness is 0.1w, subjected to bending. The length is L = 50w. Here, we do not apply homogeneous boundary conditions but fully clamp the 3D beam on one end, while a rotation angle θ = 2 rad is applied to the free end in increments of ∆θ = 0.1 rad. For this 3D FE reference analysis only, the simulation was conducted in [START_REF] Abaqus | Simulia, abaqus unified fea[END_REF], where we exploited the kinematic coupling of a face to a reference point to apply the rotation at 29 the beam's free end. Fig. 10a shows the resulting bending moment at the clamped end.

The deformed beam obtained from the 3D FE analysis compares excellently with the center-line obtained from our reduced-order model, as shown in Fig. 10b. To highlight the significantly nonlinear behavior due to the Brazier effect, the linear response resulting from classical Euler-Bernoulli 1D beam theory is included in Fig. 10a. Small deviations between the reduced-order model and the 3D FE analysis are due to boundary effects, as the beam is fully clamped. This can be mitigated by including higher-order effects in the effective 1D model, as proposed in Lestringant and Audoly (2020). Here, the leadingorder reduced model accurately captures the nonlinear effects, owing to the slenderness of the beam, even in the presence of boundary effects.

Beam assemblies

Finally, we illustrate the application of the hyperelastic beam formulation to 2D and 3D truss networks. To begin, we consider two 2D truss lattices made of (a) square and (b) hexagonal unit cells (see Fig. 11). In both cases, our reduced-order framework is compared to classical Euler-Bernoulli beam theory and 2D FE analysis under the plane stress assumption. As before, we use a square cross-section of size w (with a length-towidth ratio of the individual struts as 125 and 80, respectively) and the compressible neo-Hookean constitutive law. Fig. 11 shows the force vs. strain curves under uniaxial vertical straining with free lateral edges in 2D as well as the resulting deformation of the truss at different strain levels.

For the square truss in Fig. 11a, all three simulation approaches match perfectly in the small-strain regime. At larger strains, the hyperelastic response is softer than the Euler-Bernoulli beam, as can be expected. Across the full tested strain range, our reduced-order model and the FE simulation agree well in both the deformation and reaction force. The square mesh consists of 167,996 second-order triangular elements, whereas the beam model contains 104 beams with 10 segments per beam, thus significantly reducing the problem size and saving computational effort.

The hexagonal truss in Fig. 11b changes from a bending-to a stretching-dominated behavior at about 20% applied vertical strain: while low strain levels are primarily accommodated by the bending of individual struts, struts align increasingly with the vertical direction with increasing strain, which leads to the response being dominated by the stretching of individual struts. Such behavior was previously reported also for experiments on 2D polymeric hexagonal trusses [START_REF] Glaesener | Viscoelastic truss metamaterials as time-dependent generalized continua[END_REF]. As seen in the reaction force vs. strain curve in Fig. 11b, the Euler-Bernoulli and hyperelastic responses are identical during the bending-dominated phase, as expected from the analysis of the cross-sectional response in Section 3.1.1. As soon as the deformation mode switches to stretching-dominated, the responses of the linear theory and of the hyperelastic models differ, while the hyperelastic beam model and the FE solution agree excellently up to more than 50% applied vertical strain. Here, the 2D FE mesh consists of 85,470 secondorder triangular elements, whereas 108 beams subdivided into 20 segments per beam (i.e., 2160 beam segments in total) are used in the reduced-order model, which implies significantly reduced computational expenses.

As our last example, we consider a 3D truss consisting of 3×3×3 octet unit cells, whose beams have a solid circular cross-section made of the same neo-Hookean base material used before. The truss is subjected to tension and compression in the vertical direction, as shown in Fig. 12. The bottom face is fully clamped, while the top surface is clamped in plane and displaced vertically. Due to the high computational cost, no full 3D FE analysis was conducted here for reference, but the comparison with Euler-Bernoulli beam theory serves to confirm the small-strain response. Fig. 12 reveals that the axial strain per beam is not particularly high, which leads to a close agreement between the results from the Euler-Bernoulli beam model and our hyperelastic beam setting up to moderate strains. However, when the truss is compressed to large vertical strains, the local strains within struts are higher and the response, as seen in the force-strain curve in Fig. 12b, predicted by the hyperelastic beam model becomes significantly softer with increasing

Conclusion

We have presented a two-scale numerical framework for the efficient simulation of slender, hyperelastic beam networks, which we applied to 2D and 3D trusses in the finitestrain regime. Our approach is based on an asymptotic expansion of the deformation within slender, hyperelastic primsatic solids, which decomposes the overall deformation into slowly varying macroscopic strains on the structural level and microscopic displacements of the beam's cross-section. By pre-computing the effective strain energy density for a given cross-section and a selection of base material constitutive laws, the subsequent structural-level simulations gain efficiency. Our model works in a generic hyperelastic setting and can handle beams with inhomogeneous cross-sections with arbitrary geometry and hyperelastic (polyconvex) constitutive behavior as well as pre-strains. The effective response of a cross-section is obtained from the finite element method, whose tabulated results are used to derive an approximate effective energy density to be used in the macroscale, structural-level calculations, which exploit the discrete rod code UtoBeams (UtoBeams, 2020). We have shown that, in all examples considered here, truncating the asymptotic expansion at leading order (resulting in a local effective energy density) is sufficient, while problems involving significant size effects or localization along individual struts may necessitate higher-order expansions leading to effective strain-gradient structural models, which can be included in our numerical framework.

We have demonstrated the performance of the framework for 2D and 3D trusses made of a compressible neo-Hookean material, showing excellent agreement with full FE analysis (across the full range of investigated strain levels) as well as with Euler-Bernoulli theory (at small strains). Note that we neglected contact within hollow cross-sections, which can be incorporated in a future extension of the model. The same applies to inelastic material models, for which the general two-scale framework still applies, but internal state variables must be considered and stored in full or a reduced format across the cross-section. This is the topic of ongoing research.

  Figure 2: A summary of the asymptotic dimension reduction procedure at leading order. The full 3D elasticity problem is split into a series of microscale problems (26) set on a reference cross-section for a local value of the macroscopic strains E. The leading-order solutions ỹ⋆ E (tabulated a-priori or computed on-the-fly) enter the definition of the macroscale effective energy functional (28).

3

  Figure 3: Schematic showing the numerical two-scale procedure. A series of microscale problems is solved on a reference cross-section for a range of combinations of macroscopic strains Ēijkl ; the resulting effective energy densities are stored in a four-dimensional array and subsequently interpolated. The

  (1994);[START_REF] Himpel | Computational modelling of isotropic multiplicative growth[END_REF];[START_REF] Genet | Heterogeneous growth-induced prestrain in the heart[END_REF];[START_REF] Moulton | Morphoelastic rods III: Differential growth and curvature generation in elastic filaments[END_REF] in the context of biological growth. We follow the latter approach, where the undeformed, pre-grown configuration is considered as the reference configuration, and the calculations are performed on the intermediate, grown configuration. The total deformation gradient at leading order, F in (23), is thus decomposed as F = AG, where A represents the elastic deformation of the grown

  10.20.30.40.50.60.70.80.9

Figure 4 :

 4 Figure 4: Homogeneous cross-sections subjected to (a) axial stretching, (b) bending about d 1 , (c) bending about d 2 , and (c) twisting. All undeformed cross-sections fit into a unit square of side length w; the base material is a compressible neo-Hookean solid. The normalized effective strain energy density is plotted vs. the applied effective (squared) strain in comparison to the Euler-Bernoulli model. Images on the right compare the undeformed (dotted lines) and deformed (solid lines) cross-sections at maximum applied strain. The color-code illustrates the distribution of displacements u = ỹ⋆ Emax -ỹ⋆ E=0 .

Figure 5 :

 5 Figure 5: A stiff butterfly inclusion in a soft square matrix (E inclusion = 10 E matrix ). The geometry and the FE mesh are shown in (a). When deformed in bending with κ 1 w = 2.0, the deformed cross-section is shown in (b), and the von Mises stress distribution is shown in (c). The Euler-Bernoulli strain energy density W EB = EI 1 2 κ 2 1 and the effective energy density are compared in (d). The energies are normalized by EI 1 = n(EI 1 ) matrix + (1 -n)(EI 1 ) inclusion , where n is the mass fraction of the matrix.

Figure 6 :

 6 Figure 6: The effective energy landscape as a function of the rescaled bending strain κ 1 w for (a) a square cross-section whose top half has experienced a constant growth of 20%, and (b) a circular crosssection with a linearly varying growth applied throughout the cross-section. The color-code in the insets illustrates the axial Cauchy stress σ 33 and the von Mises stress distribution in the cross-section for the indicated levels of κ 1 and κ 2 .

Figure 7 :

 7 Figure 7: The normal force and bending moment derived from the effective strain energy density for a bi-layered square cross-section (Poisson's ratio ν, Young's modulus Etop in the top half and E bottom = 2Etop in the bottom half) vs. (a) axial strain ε, (b) bending strain κ 1 w, as compared to the Euler-Bernoulli equivalents EAε and EIκ 1 , respectively, where EA = (EA)top + (EA) bottom and EI 1 = (EI 1 )top + (EI 1 ) bottom .

Figure 8 :

 8 Figure 8: Normal force and torsional moment for pure torsion, E = (0, 0, 0, τ ), of (a) a circular beam cross-section of diameter w with stiff inclusions (E inclusion = 1000E matrix ), and (b) a homogeneous circular cross-section of diameter w and Young's modulus E matrix . Magnifications show the finite element meshes used in the examples, illustrating the refinement around the stiff inclusions.

Figure 9 :

 9 Figure 9: The resultant normal force and moment of a cantilever beam modeled as a compressible neo-Hookean solid (Young's modulus E, Poisson's ratio ν, square cross-section of size w = L/50, and length L), whose free end is (a) displaced axially and (b) twisted, respectively. Results from our reduced hyperelastic beam framework are compared to full 3D FE analysis. In (c) the von Mises stress distribution within the cross-section from the two-scale approach and the 3D FE analysis is compared for case (b) at a twist angle of 0.4 [rad], showing excellent agreement with a difference in the maximum stress of 3%. A coarser mesh was used in the microscale analysis to facilitate a direct comparison with a slice of the 3D FE simulation. For the case of stretching in (a) the uniform stress within the cross-section is identical in the two-scale and 3D FE analyses.

Figure 10 :

 10 Figure 10: The reaction moment of a cantilever tube (fully clamped on one end) modeled as a compressible neo-Hookean solid (Young's modulus E, Poisson's ratio ν, cross-section of diameter w, wall thickness 0.1w, and length L = 50w) with a rotation applied to the free end. The reaction moment at the fully clamped end of the beam is shown in (a) for our reduced-order hyperelastic beam model, full 3D FE analysis, and Euler-Bernoulli beam theory. In (b) the result from our hyperelastic beam model (blue line) is compared to the full 3D FE result (red 3D model) at a bending angle of 1.5 rad.

  Figure 11: (a) Square and (b) hexagonal lattices in 2D made of struts of square cross-sections (side length w, lengths 125w and 80w, respectively) are loaded vertically in tension. The net reaction force at the clamped edge is compared for the reduced-order hyperelastic model, the Euler-Bernoulli model, and 2D FE analysis. The deformed trusses on the right are at (A) 0 %, (B) 15 %, and (C) 40 % vertical strain. The color code highlights the axial strain per beam segment, while the close-ups compare the deformation between 2D FE and our beam model.

Figure 12 :

 12 Figure 12: 3D trusses, composed of octet unit cells, are fully clamped at the top and bottom faces and loaded uniaxially (a) in tension and (b) in compression. The reaction forces are compared as obtained from our hyperelastic beam model and from Euler-Bernoulli beam theory. Individual struts have a circular cross-section of diameter w and a length of 200w. The axial (stretching and compression) strains per beam segment are shown by the color code.

Note that not all four quantities in E have units of strain, as is apparent from their definitions and as will be discussed in the following. Yet we refer to those as the macroscopic "strains" for conciseness.

Acknowledgements

We gratefully acknowledge Manuel Weberndorfer's support with the integration of the material model in UtoBeams.