Viscoelastic truss metamaterials as time-dependent generalized continua

Mechanical metamaterials provide tailorable functionality based on a careful combination of base material and structural architecture. Truss-based metamaterials, e.g., exploit structural topology and beam geometry to achieve beneficial mechanical and physical properties from stiffness and wave dispersion to strength and toughness. While the focus to date has been primarily on static metamaterial properties or elastic wave motion, 3D-printed polymeric base materials naturally come with significant viscoelasticity, making the effective truss response time-and rate-dependent. Here, we report a theoretical-numerical-experimental study which (i) deploys a linear viscoelastic corotational beam description (capturing finite rotations at small strains), (ii) implements the latter in a finite element framework, (iii) calibrates a generalized Maxwell model based on viscoelastic experiments on 3D-printed polymer samples, (iv) validates the theory and implementation through experimental truss benchmark tests, and (v) introduces a generalized continuum formulation for the efficient simulation of viscoelastic truss metamaterials containing large numbers of structural members. We show that the viscoelastic beam approach, calibrated via tension tests on individual strut samples, performs well when applied to complex truss lattices undergoing time-dependent stress relaxation — as verified by the effective mechanical response and full-field deformation maps. The resulting variational generalized continuum framework uses on-the-fly periodic homogenization based on a representative unit cell and is extended to dynamics by including inertial effects. By comparison to discrete numerical simulations we demonstrate the accuracy of the continuum approach, which is promising for modeling and optimizing 3D-printed truss metamaterials for engineering applications from shock-absorbing structures to rate-dependent architected materials and soft robotics.


Introduction
Architectural design plus base material determine the effective mechanical properties of metamaterials based on, e.g., trusses (Wadley et al., 2008;Meza et al., 2015), plates and shells (Zheng et al., 2018;Bonatti and Mohr, 2019), or spinodal structures (Vidyasagar et al., 2018;Kumar et al., 2020).The available design space has enabled metamaterials with optimized, application-dependent static (Surjadi et al., 2019) and dynamic properties, the latter including shock mitigating by impact energy absorption (Schaedler et al., 2014;Frenzel et al., 2016) and wave guiding (Zelhofer and Kochmann, 2017).As additive manufacturing is the primary means of fabrication, versatile 3D printing materials have been used as base materials for metamaterials, which requires detailed understanding of their mechanical behavior to accurately design, model, and optimize the engineered metamaterial response.While for additively manufactured trusses made of metallic (Zadi-Maad et al., 2018) and ceramic base materials (Zocca et al., 2015) the assumption of linear elasticity is oftentimes sufficient, compliant polymers require rate-dependent models to describe the complex time-dependent behavior and history dependence (Gibson et al., 2010;Wang et al., 2017).Although being a technical complication for modeling, the viscoelastic relaxation, creep, frequency-dependent harmonic or generally time-dependent response of such architectures facilitates fascinating new features and applications (Gomez et al., 2019) ranging from rate-dependent buckling patterns and bistable beams (Dykstra et al., 2019;Janbaz et al., 2020) to advanced vibration damping (Wang and J., 2018) and frequency control and dispersion tuning in soft phononic crystals (Frazier and Hussein, 2015;Parnell and De Pascalis, 2019).Seizing the opportunities provided by viscoelastic metamaterials requires accurate and efficient models for the design space exploration and property optimization.
While the classical theory of viscoelasticity was laid out long time ago (see e.g.Schapery (1969), Taylor et al. (1970), Christensen (1982), Lakes (1999)), the modeling of viscoelastic structural members is a more recent topic, including, e.g., finite element (FE) models for viscoelastic plates (Marques and Creus, 1994;Yi and Hilton, 1994), beams (Hilton, 2009), and sandwich structures (Galucio et al., 2004).Rate-dependent trusses were studied, among others, by Ghayesh et al. (2016) who studied the viscoelastic response of a single Euler-Bernoulli beam, whereas Hamed (2012) focused on a 2D Timoshenko beam, and Bottoni et al. (2008) modeled linear viscoelastic thin-walled beams.Recently, Ananthapadmanabhan and Saravanan (2020) reported numerical techniques to calculate the response of nonlinear viscoelastic truss networks (yet ignoring bending strains).In addition, Lestringant et al. (2020) and Lestringant and Kochmann (2020) introduced a general formulation of slender, geometrically exact beams made of viscoelastic base materials, which has not yet been applied to trusses.
With the objective of providing a simple and accurate truss description and an efficient finite element treatment, we here extend the corotational beam element of Crisfield (1990) to linear viscoelastic beams, which accounts for rate dependence in the axial, flexural and torsional stress-strain relations (based on a general one-dimensional linear viscoelastic constitutive law) and accounts for finite rotations.Implementation of the model in a FE framework allows us to simulate the linear viscoelastic response of complex truss networks in two and three dimensions (2D and 3D, respectively).Validation is realized by comparison to relaxation experiments on 3D-printed 2D hexagonal trusses, whose generalized Maxwell model parameters are extracted from tensile tests.The comparison between simulations and experiments on the strained hexagon lattice shows excellent agreement, which confirms the applicability of our viscoelastic corotational beam formulation to polymeric trusses.
Even though we present an efficient truss description, the cost of solving boundary problems involving thousands to millions of struts inside a truss becomes prohibitive, so multiscale techniques become beneficial (Kochmann et al., 2019).Here, we incorporate the linear viscoelastic beam into a generalized continuum formulation, which was previously introduced for linear elastic beam networks (Glaesener et al., 2019(Glaesener et al., , 2020)).The discrete truss is replaced by a continuous body, whose effective mechanical constitutive behavior is obtained from on-the-fly numerical homogenization of a representative unit cell, and whose deformation is solved by a finite element calculation on the macroscale (comparable to FE 2 ).This two-scale model captures finite strains on the macroscale (accommodated by finite rotations of truss members on the microscale), and it applies to stretching-and bending-dominated truss topologies.Unlike Vigliotti et al. (2014) and Pal et al. (2016), who pursued similar approaches in 2D but with rotational degrees of freedom condensed out on the microscale, we introduce both translational and rotational degrees of freedom on the macroscale and pass those to the representative unit cell.We compare simulation results of discrete numerical calculations to those of the two-scale generalized continuum approximation for several 3D relaxation and vibration benchmarks, using a selection of bendingand stretching-dominated truss topologies, which demonstrate the homogenization scheme's accuracy and efficiency as well as its applicability.In addition, our results highlight the exciting opportunities for exploiting viscoelastic base materials for the design of time-dependent metamaterial properties.
The remainder is structured as follows.Section 2 lays out the linear viscoelastic beam theory, introduces its FE implementation through corotational beam elements, and summarizes the formulation of the two-scale generalized continuum representation of viscoelastic trusses.Section 3 validates the viscoelastic truss model by comparing simulated results to experimental relaxation tests on 2D hexagonal trusses.3D quasistatic stress relaxation and dynamic vibration examples in Section 4 demonstrate good agreement between fully resolved discrete numerical calculations and efficient FE simulations based on the two-scale representation, before Section 5 concludes our study.

Linear viscoelastic corotational beams
We consider slender beams made of a linear viscoelastic material in 3D, which we describe by an extension of the corotational beam formulation originally introduced by Crisfield (1990) for linear elastic Euler-Bernoulli beams.Owing to the linear constitutive relations and the resulting applicability of Boltzmann's superposition principle, the formulation and numerical implementation of our viscoelastic beams follows that of their linear elastic counterparts and, particularly, admits the decoupling of stretching, flexural and torsional stress and strain components.

Viscoelastic beam theory
We assume our corotational beam to be sufficiently slender, so it experiences a linear combination of axial strains due to stretching and bending, on the one hand, and shear strains due to torsion, on the other hand, while we neglect all further strain components.Let the -axis be locally aligned with the neutral axis of the beam (Fig. 1(a)), so that for every point  = (, , ) T on the beam at time .We describe the displacement of the beam's center-line by the displacement field ) is an axial strain in the beam that is constant across the cross-section, whereas the second and third terms in (1) represent axial strains due to bending about the -and -axis, respectively, with curvatures   (, ) =  ′′  (, ) and   (, ) = − ′′  (, ).Here and in the following, primes denote derivatives with respect to .The only non-zero shear strain is due to torsion: where (, ) is the twist angle of the beam.Owing to the slender-beam assumption,   ,   , and  are constant per cross-section and hence depend only on the position  along the beam's center-line.
The linear viscoelastic constitutive relations relate the above strains to the only two non-zero stress components at every point  in the beam's cross-section: where dots denote derivatives with respect to time.For convenience, we approximate the relaxation functions by Prony series, defining with the long-term Young and shear moduli  ∞ and  ∞ , respectively, as well as moduli   (  ) and associated relaxation times   (  ), assuming a total of  spring-dashpot elements for each relaxation function (Fig. 1(b)).For simplicity, we assume that the beam is homogeneous, so that all material properties are constant within the beam.The generalized Maxwell (or Wiechert) model admits recasting the above integral representation (3) into the following form (Simo and Hughes, 2000;Marques and Creus, 2012).For the axial strain component we introduce  internal variables   p , and analogously   p for the torsional strain ( = 1, … , ), so that the elastic axial and torsional strains inside each of the  Maxwell elements are   −   p and   −   p , respectively.Consequently, ] . (5) The internal variables evolve at each point in the beam according to the kinetic relations ] .
(10) Inserting ( 7) into (5), applying the finite difference discretization and applying (10) leads to an explicit solution for the stress at the new time  +1 : and analogously The normal force , the two bending moments   and   as well as the torsional moment   follow from integration over the cross-section as The above model can be cast into a variational framework by exploiting variational constitutive updates (Ortiz and Stainier, 1999;Lestringant et al., 2020).Integration of the effective energy density of Ortiz and Stainier (1999) over the cross-section leads to an incremental energy density per beam length, which depends on the internal variables from the previous time step only: with the area and area moments of inertia defined by This completes the description of our linear viscoelastic corotational beam.Without finite rotations, the above framework admits closed-form solutions for monotonic loading and various boundary and initial conditions.Under finite rotations, the framework still applies (Crisfield, 1990;Crisfield and Cole, 1990) but the -axis refers to the current, deformed center-line of the beam, so that the resulting nonlinear governing equations require a numerical formulation for solving initial boundary value problems.

Finite element setting
Following the corotational beam formulation of Crisfield (1990) and its numerical implementation by Phlipot and Kochmann (2019), we represent a beam as a discrete set of segments separated by nodes, whose deformed and undeformed locations in space are given by, respectively, position vectors  ∈ R  and  ∈ R  in  dimensions.The deformed and undeformed nodal orientations are defined, respectively, by rotation vectors  ∈ R 2−3 and  ∈ R 2−3 .Altogether, the degrees of freedom of a node  are abbreviated by the vector of nodal degrees of freedom, x = (  ,   ) T ∈ R 3(−1) .For a beam segment with nodes at positions  1 and  2 in the current configuration, the deformed orientation of the beam is defined by  = ( 2 −  1 )∕ with the deformed segment length | denotes the undeformed segment length).To account for finite rotations of the beam that do not contribute to the beam's deformation, we define by   the rotations of node  relative to the current beam orientation , defined separately for each segment (for details see, e.g., Phlipot and Kochmann (2019)).
Let    and    denote the angles of rotation at node  about the -and -axes, respectively, with respect to the current orientation  of a segment, and let    denote the twist angle at node .The FE interpolation is based on linear ansatz functions for the axial displacement   () and twist (), while third-order Hermite polynomial ansatz functions are used for the transverse displacements   () and   () of the beam's center-line (which includes the exact solution for a beam loaded at its end points).Following a Rayleigh-Ritz approach, inserting those ansatz functions (with the essential boundary conditions enforced at nodes) into ( 14) and integrating over the beam segment yields the nodal reactions by differentiating the resulting segment energy with respect to the nodal degrees of freedom.This includes the axial forces with moduli where we point out that the axial strain is constant along the beam segment, so one internal variable   ax.,p per segment is sufficient (and the same applies to the torsional internal variables   p ). Assuming that the coordinate axes align with the principal axes of the beam's cross-section (so   = 0), the nodal bending and torsional moments are obtained as with moduli For the linear viscoelastic bending contribution, one set of internal variables { , ,p ,  , ,p } is required per node  to capture the history dependence.All nodal reactions must be rotated from the deformed frame to the undeformed reference frame within a Lagrangian FE setup (see Crisfield (1990)).
The above linear viscoelastic corotational beam element is implemented within a massively-parallel in-house computational mechanics code, so that we may subsequently use it either for discrete numerical calculations or within the generalized continuum framework described in the following.For quasistatic simulations, a Newton-Raphson solver is used, whereas dynamic problems are solved by a Newmark- implicit time integration scheme (Newmark, 1959) with the lumped mass matrix (Iura and Atluri, 1995) for the linear beam element.A validation example is presented in Appendix A, which simulates a viscoelastic cantilever beam utilizing both the above linear viscoelastic corotational beam formulation and a fully-resolved FE description based on solid elements for comparison.

Modeling viscoelastic trusses as generalized continua
We deploy the above linear viscoelastic beam description in our recently introduced generalized continuum framework (Glaesener et al., 2019(Glaesener et al., , 2020)), previously used only for static, linear elastic beam networks.The objective is to considerably reduce the computational costs associated with fully resolved discrete calculations involving large numbers of beams.Our strategy is to replace the periodic truss by an effective continuum whose response is approximated by 3D solid finite elements with the effective constitutive behavior on the macroscale obtained from a nested FE calculation on a representative unit cell (RUC) on the microscale -assuming a separation of scales.Efficiency is gained by our FE 2 -type formulation (Glaesener et al., 2019), which uses periodic homogenization on the microscale, while allowing for consistent analytical tangents to be used in implicit macroscale FE simulations in 3D (Glaesener et al., 2020).Fig. 2 illustrates the general procedure, which we briefly summarize here.On the macroscale, the deformation of a body  ∈ R  in  dimensions is described by a deformation mapping (, ) ∶ ×R → R  and a rotation field (, ) ∶  × R → [0, 2] 2−3 in  ≥ 2 dimensions, such that the deformed position of point  ∈  at time  is  = (, ).The deformation gradient tensor is  = Grad , while the curvature tensor is  = Grad  (with gradients defined with respect to the undeformed configuration).The finite-strain formulation is required to capture finite rotations.Each required constitutive model evaluation on the macroscale at a quadrature point   passes the kinematic quantities (  ) = {(  ),  (  ), (  ), (  )} down to a RUC, which consists of the smallest periodic unit representing the effective truss response (either a single unit cell or, if buckling gains importance, an array of unit cells, sufficiently large to capture the relevant buckling modes; see Glaesener et al. (2019)).In a nested RUC calculation, (  ) is applied on average to the RUC while enforcing periodic boundary conditions.Owing to the variational structure of our model, the resulting RUC energy is computed as the volume-averaged RUC energy, whose derivatives provide all required stress (and couple stress) measures on the macroscale along with their consistent tangent matrices.The interested reader is referred to Glaesener et al. (2019Glaesener et al. ( , 2020) ) for further information and implementation details (for linear elastic trusses).The extension to linear viscoelasticity implements internal variables in each beam segment within the RUC and hence stores those variables independently for each quadrature point on the macroscale.Otherwise, the numerical setup for linear viscoelastic trusses is analogous to that of linear elastic trusses.As an explicit form for the updated internal variables is available and no numerical condensation is required, the numerical efficiency is marginally affected by viscoelasticity, while computational storage is increased due to the requirement for storing internal variables.
As a further extension of our generalized continuum approach, we consider dynamic effects, which requires the introduction of effective inertia within the macroscale FE problem.We introduce the relative density (fill fraction) of the RUC as where  , is the th beam's cross-sectional area,   its undeformed length, and  the density of the beam's base material. RUC is the RUC volume, and we ensure that beams lying on the RUC boundary are uniquely accounted for.The total mass of a finite element  of volume   follows as   =  *   , which in turn defines the translational inertia in the lumped mass matrix of the macroscale element as where   represents the number of nodes of the macroscale element.The rotational inertia of the effective lumped mass matrix of a macroscale finite element is established by matching the kinetic energy   of element  with that of the underlying RUC,  RUC , each normalized by its respective volume (  ∕  =  RUC ∕ RUC ), which yields where  , is the (diagonal) mass moment of inertia tensor of the th beam in the RUC and is the angular velocity vector associated with its two end points 1 and 2. Analogously,  rot. is the effective (diagonal) lumped mass matrix of the macroscale finite element associated with rotations only, with  = ( φ1 , … , φ  ) containing the angular velocity vectors of all   nodes of the element.In 3D, rotations are accounted for about all three axes, whereas in 2D only a single (in-plane) rotational degree of freedom must be considered.The thus-obtained lumped mass matrix of an individual beam and of the macroscale finite element (including both translational and rotational contributions) is provided in Appendix B along with an illustration of the principle of matching kinetic energies between macro-and microscales.
We acknowledge that this approximation of matching the effective inertia in the continuum representation is only suitable in the long-wavelength limit, i.e., for relatively slow dynamic processes that engage the inertia of the overall truss rather than exciting resonances of individual truss members (e.g., we do not claim to capture wave dispersion with this approach).Yet, it will prove effective in computing the dynamic and vibrational response of large truss structures at low frequencies.The generalized continuum model is implemented inside an in-house 3D FE framework and uses 8-node brick elements with tri-linear interpolation along with the nonlinear Newton-Raphson solver from the Portable, Extensible Toolkit for Scientific Computation (PETSc) (Abhyankar et al., 2018;Balay et al., 1997Balay et al., , 2019Balay et al., , 2020)).Dynamic simulations use implicit Newmark- time integration of the FE equations of motion.

Model calibration and validation by experiments
Before assessing the continuum approximation, we demonstrate the applicability and accuracy of the viscoelastic beam formulation of Sections 2.1 and 2.2 by simulating a discrete, 3D-printed polymeric truss.We choose the representative example of a 2D hexagonal truss lattice made of a thermoplastic elastomer by selective laser sintering (SLS) on a Sintratec S1 printer.The hexagonal topology combines stretching and bending deformation, making it a valuable benchmark.(We deliberately choose a 2D example that admits detailed extraction of the truss deformation, rather than using a 3D truss with the added complication of imaging a cellular 3D structure).Fig. 3 illustrates samples clamped at their top and bottom and stretched under displacement control to total vertical strains of 8.1% (experiment 1) and 14.3% (experiment 2), leading to viscoelastic stress relaxation over time.To calibrate the Prony series of the base material, we performed separate uniaxial extension experiments on dog-bone samples and fitted the simulated response force vs. time of the generalized Maxwell model to the experimental force-time history, leading to the Prony series detailed in Appendix C along with all material parameters.
The undeformed truss in Fig. 3 measured 94mm × 55.4mm with individual beams being  = 8mm long with a measured beam width of  = 0.66mm and hence a slenderness ratio of  =  ∕ = 0.0825.The truss had a constant out-of-plane thickness of  = 5.0mm, so that deformation is approximately 2D (and out-of-plane displacements can be neglected in simulations).The sample was mounted on a tensile testing machine (Instron Model 5943, 1kN force capacity) using custom clamps, while a 450N load cell (mounted to the bottom clamp) measured the total vertical force applied to the truss.A high-speed camera (Photron AX200) was used to image the truss during deformation at a frame rate of 50 frames per second.Images were processed to track the center-lines of individual beams and the positions of nodes, using a custom-built image processing code.We define the (average) vertical strain of the truss as the relative displacement of the top-and bottom-most nodes on the vertical center-line of the sample normalized by the initial height of the sample (see Fig. 3).Note that, before applying any tensile strain, the starting configuration was already pre-stretched due to gravity (and this pre-stretched truss is used as the initial configuration in simulations).
Two sets of experiments were performed with maximum tensile strains of  = 8.1% and  = 14.3%; strain histories are shown in Fig. 4.Each experiment starts with a constant uniaxial displacement rate of 100 mm∕min applied to the clamped top edge of the truss (section a ⃝b ⃝ in Fig. 4), followed by the strain being held constant ( b ⃝c ⃝), until the structure is brought back to its original configuration ( c ⃝d ⃝) at a rate of 10 mm∕min in experiment 1 and 100 mm∕min in experiment 2 (all circled labels refer to the points marked in Fig. 4).
A visual inspection of the deformed trusses at the respective maximum strains, shown in Fig. 5, indicates a remarkable match between experimental observations and simulation results.The total vertical forces measured during the strain paths of Fig. 4 are summarized in Fig. 6.After an initial nonlinear rise of the tensile force, typical viscoelastic relaxation behavior is observed ( b ⃝c ⃝), before unloading sets in.Simulations use the linear viscoelastic corotational beam element (with 10 elements per strut) and the Prony series of Table C.3, and they demonstrate an overall convincing agreement with experiments.The error between numerical  and experimental data, which is almost constant during the relaxation phase, is attributed to imperfections in the printed material and to nodal effects (Portela et al., 2018), which are not considered in simulations (nodes are rigid in our corotational setup).In addition, due to experimental constraints we lack knowledge of the Prony series for faster relaxation times (as explained in Appendix C).
We evaluate the accuracy of the simulations for experiments 1 and 2 by the normalized root-mean-square error (NRMSE) where F are experimentally measured forces (for all  time steps) and   are the ones obtained from the discrete model;  denotes the total number of data points.Normalization is performed over the total range of forces in the experiment (difference between maximum and minimum force values: F − F ).The NRMSE between a ⃝ and d ⃝ evaluates to ∼ 4.0% and ∼ 15.48% in experiments 1 and 2, respectively.The error for experiment 1 is remarkably low, considering the moderate strains and resulting nonlinear deformation as well as the uncertainty in printing-induced imperfections.As expected, the error grows with increasing applied strain (cf.experiment 1 vs. experiment 2), which we attribute to increasing nodal effects and the decreasing legitimacy of the linearized-kinematics beam model -even though the deformed trusses in Fig. 5 show convincing agreement even at 14.3% strain.When stretching the structure up to 14.3% in experiment 2, especially those beams connected to the clamps undergo large strains that are likely to violate the small-strain assumption.The identically shaped but stiffer numerical force response leads to the conclusion that the elastic modulus  ∞ is slightly lower in experiments, whereas the contributions of the  dashpots in the Maxwell model match measurements well.A detailed analysis of the local strain and curvature values within individual struts in both experiments is provided in Appendix D, which reveals that local strain measures are significantly higher in experiment 2, which invalidates the assumption of linearized kinematics, hence leading to the observed larger error for experiment 2.
We point out that, like in any FE discretization, simulated results converge with ℎ-refinement, as demonstrated in Fig. 7, presenting the total force vs. time for the hexagonal truss relaxation in experiment 1 for six different mesh resolutions (ranging from a single element per strut to 10 beam elements per strut).With increasing refinement, bending deformation is increasingly concentrated at nodes, which leads to larger rotations and hence higher bending moments at the nodes.As a consequence, we observe convergence from below.As we deem 10 elements per strut to show a sufficient level of convergence, we used this resolution in our comparison to experiments in Fig. 6.Overall, we see convincing agreement between simulations and experimental measurements -both in the full-field deformation maps of Fig. 5 and the force-time histories in Fig. 6, based on which we use the viscoelastic truss framework in the following to simulate more complex examples.

3D simulations and the generalized continuum approximation
Having verified the applicability of the viscoelastic beam model for 3D-printed trusses, we proceed to assess the validity of the generalized continuum approximation introduced in Section 2.3 as a reduced-order model for efficient truss simulations.In the following, we refer to the FE 2 -type framework resulting from viscoelastic corotational beams used in the homogenized RUC employed at every quadrature point of the macroscale FE problem as two-scale model.As benchmarks, we illustrate the quasistatic relaxation of a 3D truss cube after torsional twisting as well as the dynamic response of the same truss undergoing damped free vibrations.
All struts in the trusses are described as viscoelastic corotational beams with a circular cross-section of diameter , which gives the cross-sectional area   =  2 4 , area moment of inertia  =  4 64 , and polar moment of inertia   = 2.To cover a range of stretching-and bending dominated trusses, we analyze the response of four different truss topologies based on cubic, octet, octahedral, and bitruncated octahedral unit cells (Portela et al., 2018).For a fair comparison, we ensure a constant effective density, while individual struts are ensured to remain below a maximum slenderness ratio of  = ∕ ≤ 0.1, so that Euler-Bernoulli beam theory applies.The chosen sets of parameters that satisfy those constraints are presented in Table 1.
For each of the four topologies, we simulate the truss response both via discrete FE calculations (resolving all struts by corotational beam elements as in Section 3) and using the two-scale model within an efficient macroscale FE calculation based on a mesh of brick elements.The highly slender beams in the octahedral topology cause the effective homogenized body to behave in an approximately incompressible manner (Phlipot and Kochmann, 2019), which causes convergence challenges and calls for specialized element types on the macroscale.We previously showed (Glaesener et al., 2020) that the two-scale approximation indeed captures instabilities and localizations well when applied to linear elastic trusses.To avoid the complexity of dynamically snapping structures, we here limit our benchmarks to relatively small strains at the macroscale which, however, yields representative examples that verified the applicability of the dynamic, viscoelastic continuum representation.For all benchmarks, we continue to use the previously calibrated Prony series of Table C.3.

Torsional relaxation test
We consider a cube filled by 40 × 40 × 40 truss unit cells (each of side length    = 10mm), which is undergoing torsional deformation.While the bottom of the cube is clamped to a rigid substrate, the top is free to move vertically but forced to undergo a prescribed in-plane twisting motion.A visual comparison of the fully resolved discrete octet structure and the generalized continuum approximation used in the two-scale model is shown in Fig. 8.
After twisting the top at a constant rate of 72 • ∕min up to a maximum angle of 3 • (section a ⃝b ⃝), the top is held fixed ( b ⃝c ⃝) for 5s, before being rotated back to its initial state at a rate of 72 • ∕min ( c ⃝d ⃝).Fig. 9 compares the average applied torque vs. time curves for both a fully resolved FE simulation of the discrete truss and the continuum approximation based on a mesh of 10 × 10 × 10 brick elements and the smallest possible RUC    for all calculations, shown within the plots of Fig. 9.
Analogous to the results of the 2D relaxation experiment in Section 3, the viscoelastic base material leads to the relaxation of internal stresses and, as a result, to a decrease in the applied torque while keeping the twist angle constant.The torque-vs.-timecurves for the different topologies in Fig. 9 show convincing agreement between the fully-resolved discrete calculations and the continuum approximation.We note that a proper representation of the cubic topology (Fig. 9d) within the two-scale approach requires an enlarged RUC, assembled from 2 × 2 × 2 cubic unit cells.We attribute the too stiff response obtained form only a single cubic unit cell (1 × 1 × 1) to the suppression of bifurcation modes that arise in the soft, bending-dominated cubic topology at relatively low strains -comparable to what Glaesener et al. (2019) showed for hexagonal and triangular linear elastic structures in 2D.Interestingly, we observe an almost identical response of the stretching-dominated octet topology (Fig. 9c) and the octahedron (Fig. 9b), which is due to their identical effective stiffness matrices when both trusses have the same relative density (see Table E.4).
The error between the discrete truss calculation and the homogenized FE approximation strongly depends on the mesh size of the macroscale FE problem, as shown in Fig. 10.With ℎ-refinement of the macroscale brick mesh, we see convergence from above of the error in the computed net torque to the solution of the discrete calculation.We note that, unlike in linear elasticity, the linear viscoelastic model's dependence on the previous time step's internal variables results in an accumulation of errors over time, so that an inappropriate FE mesh resolution may lead to significant errors over time.

Damped vibrations
Let us reconsider the cube-twist example of Section 4.1, now with the additional consideration of inertia.Instead of slowly twisting the cube after relaxation back in a displacement-controlled fashion, the truss is now released to dynamically return to its initial configuration under damped vibrations -with significant damping stemming from the viscoelastic base material (we ensure no numerical damping from the Newmark- scheme).Specifically, we twist a 30 × 30 × 30 unit-cell truss at a rate of 89 • ∕min up to 3.7 • , then allow the structure to relax until  = 5s, and then remove all constraints at the top face instantaneously, leading to vibrations.The evolution of the average twist angle of the top face over time is shown in Fig. 11 for each of the four truss topologies.
As in the quasistatic examples, the results obtained from the homogenized continuum approximation show convincing agreement with those from the considerably more expensive fully-resolved truss calculations, confirming that inertial effects within the RUC at the microscale are successfully passed to the macroscopic FE problem.The cube topology again necessitates the use of a 2 × 2 × 2 unit cell for an accurate agreement, so we use the same RUCs shown in Fig. 9. Also analogous to Section 2.3, the octet and octahedron show identical oscillations due to the same stiffness matrix and almost identical homogenized lumped mass matrices (only those terms associated with torsional inertia differ, yet those play only a minor role compared to the other inertial terms, see Appendix E). Results are presented for four types of truss topologies (bitruncated octahedron, octahedron, octet, and cube), all using the same linear viscoelastic base material based on the Prony series of Table C.3.All trusses have the same relative density (cf.Table 1).
Note that the underlying Prony series was identified from experiments without consideration of relaxation times  < 0.3s, which plays a significant role in high-frequency oscillations, so a re-calibration of the Maxwell model may be required to accurately account for such behavior if needed.

Conclusions
We have presented a framework of linear viscoelastic slender beams along with their numerical implementation as 3D corotational beam elements (undergoing finite rotations while considering linearized axial, flexural and torsional strains).By formulating the viscoelastic corotational beam in terms of reduced internal variables (one per each of the aforementioned strain contributions), the model becomes computationally tractable.A generalized Maxwell model was calibrated experimentally for the thermoplastic base material of 3D-printed truss by identifying the parameters of the associated Prony series from quasistatic uniaxial relaxation experiments.We have validated the model by comparing the simulated total force-time history and full-field deformation maps of 2D hexagonal truss lattices under tension to those of 3D-printed trusses.The comparison showed convincing qualitative and quantitative agreement with a normalized root-mean-square error as small as ∼ 4.0% for experiments with moderate strains.octahedron, octet, cube) in comparison to the results of the two-scale model used in a 10 × 10 × 10 brick mesh.The cube is initially twisted up to 3.7 • at a rate of 89 • ∕min and held at that twist angle for 2.5 s to allow the material to relax, and then released to return to the initial state through damped vibrations.All trusses have the same linear viscoelastic base material model with the parameters of Table C.3 and the same effective density according to Table 1.For comparison, dashed lines illustrate the quasistatic relaxation response of each structure without considering inertial effects.
2009).To estimate the validity limits of the linearized kinematics formulation underlying the corotational beam formulation, Abaqus simulations were performed using both updated-Lagrangian (UL) and total-Lagrangian (TL) settings -the former providing higher accuracy in case of increased deformation levels.As a worst-case scenario, we use a single corotational beam element within the discrete calculations, expecting the error to further decrease with refinement.Note that the use of Abaqus requires a representation of the Prony series in terms of the instantaneous elastic shear modulus  0 , viz.   , and torsional moments 12(c) from both types of simulations for the three loading scenarios with particular focus on the relaxation behavior.In tension (without lateral constraints), the beam deforms homogeneously.In bending and torsion, we fully clamp the cross-sections at both ends for a more realistic deformation mode.As done in Fig. 6, we evaluate the NRSME over the entire load history.The normal force shows convincing agreement for small axial strains (NRMSE ∼ 4%) but loses accuracy when the strain level increases (NRMSE > 10%).This effect is not observed under bending and torsion, where the accuracy remains high even up to large strains and large twists.
We further study the influence of the slenderness ratio  on the accuracy of the linear viscoelastic corotational beam description.The deformation within each of the three load case scenarios is kept constant, while the beam diameter is increased and the NRSME between the beam response and the FE results is measured.As expected, Fig. A.13 shows a constant error in case of tension (since the applied deformation is homogeneous and the NRSME therefore independent of the diameter).In bending, the NRSME increases with slenderness ratio for both UL and TL settings.Lastly, the error for the twisted beam increases with slenderness ratio  on the accuracy of a twisted beam when using the updated-Lagrangian formulation (while being marginal and constant trend when using the total-Lagrangian setting).
The presented comparison of the discrete beam calculation to a fully resolved FEM computation also underlines the efficiency of the proposed model concerning the required degrees of freedom and computational costs.On average, the solid FEM calculation with 6501 linear tetrahedral elements takes ∼ 350 times longer than the same calculation using a single corotational beam with the linear viscoelastic description (one full run is subdivided into 100 load steps).Especially when simulating large truss networks with several thousands of beams, we can achieve almost the same accuracy with the proposed beam model as with an expensive solid FEM calculation.

Appendix B. Mass matrix of the corotational beam and effective mass matrix of the macroscale finite element
We adopt the lumped mass matrix for an individual corotational beam element from Zuo et al. (2014).The translational lumped inertia is derived by dividing the total mass of the beam   =    on both nodes.The definition of the inertial terms associated with the rotational degrees of freedom are not unique, which is why the beam's kinetic energy is approximated by rigid body motions about each axis.This procedure results in the lumped mass matrix for the 3D corotational beam, which reads for each node in 3D (so the full lumped mass matrix is obtained by assembling the below block matrix for each of the two nodes of a corotational beam element) where  and   are the beam's length and cross-sectional area, and  denotes the base material's mass density.For simplicity, this mass matrix is constant (i.e., we do not take into account changes to the above mass matrix due to finite rotations of the beam, which can be accounted for in an extension).
The mass matrix of the individual beam is used to compute the kinetic energy of the RUC, which in turn is matched with the kinetic energy of the macroscale finite element to provide the latter's (lumped) mass matrix.Fig. B.14 schematically illustrates the concept of obtaining the effective rotational inertia in 2D.An average angular velocity  is applied to every node of the macroscopic (constant-strain) finite element as well as to each node in the RUC, whose matching kinetic energies yields the effective mass matrix of the macroscale element.Uniform -values applied to all nodes of an element result in  also being active at each quadrature point (and ∇ = ), which in turn results in the same  being applied to every node in the RUC.In 3D the kinetic energy matching procedure is repeated for each component of the angular velocity vector, so that three independent components of the lumped mass matrix are obtained.(In all our examples in 3D, the symmetry of the RUC leads to identical values about each of the three axes.)The translational components of the macroscale finite element simply follow from mass lumping onto the nodes of the macroscale element.Overall, the lumped mass matrix of the macroscale finite element for each node is then obtained as (to be duplicated for 12. Net normal force , transverse force , and torsional moment  vs. time of for a cantilever beam of  = 10mm and diameter  = 1mm.All plots compare the response of the presented corotational beam element to a detailed FEM simulation using 6501 linear tetrahedral solid elements in Abaqus (Smith, 2009), based on updated-Lagrangian (UL) and total-Lagrangian (TL) schemes.The deformation is ramped up linearly over 2.5s, then kept constant for 5.0s, until the beam is returned to its initial configuration (analogous to the hexagonal truss study).Tables indicate the NRSME between results from the corotational beam description and the 3D FE simulation.(c) includes a visualization of the fully resolved beam modeled using Abaqus.

Table E.4
Effective elastic moduli of the octahedron and octet lattices with the same relative density of 0.833% (all moduli are normalized by the base material's Young's modulus  ∞ ), computed by using the homogenized continuum framework.axial strains by a factor of ∼ 10 and in the maximum curvature by a factor of ∼ 3.5 in experiment 2 (even though the applied average vertical strain increases only a factor of ∼ 1.7).This strong increase in beam elongation and bending highlights the emergence of nonlinearity, which is why the Euler-Bernoulli-type beam model loses accuracy when simulating experiment 2. It is for this reason that we observe a strong increase in the NMRSE from experiment 1 to experiment 2.More accurate results can be obtained by (i) introducing finite kinematics in the beam formulation and (ii) extending the presented theory in Section 2.1 to a nonlinear viscoelastic constitutive response.(Both go beyond the scope of this study.)

Appendix E. Effective incremental elastic moduli and homogenized mass matrix of the octahedron and octet lattices
As discussed in Section 4, the effective response of the octahedron and octet lattices show only marginal differences, when considering slender struts and low relative densities as those chosen here.To underline the similarities in their mechanical behavior, Table E.4 summarizes their effective elastic moduli, computed from the generalized continuum model.
The homogenized lumped mass matrix of the two topologies also shows noticeable similarities, when the trusses have the same relative densities -with the exception of those components associated with torsional inertia.The latter for the octet is about half that of the octahedron.However, since torsional inertia of individual beams plays only a minor role in all examples reported here where axial and bending modes dominate, its effect on the damped vibrations in Fig. 11 is negligible.For comparison, the two diagonal lumped mass matrices (using the geometric and material parameters of all truss simulations, cf.

Fig. 1 .
Fig. 1.The strain distribution within a slender 3D corotational beam decomposes into axial, flexural, and torsional strains, whose resultant stresses are modeled by a 1D generalized Maxwell constitutive model (represented by a Prony series with  Maxwell elements).

Fig. 2 .
Fig. 2. Schematic of the multiscale formulation: the macroscale finite element extracts the effective constitutive behavior of the continuum approximation from a nested structural problem solved on the RUC-level with periodic boundary conditions (each beam modeled as a viscoelastic corotational beam).

Fig. 3 .
Fig. 3. Snapshots of experiments 1 (left) and experiment 2 (right), each comparing the respective undeformed to the deformed truss at the maximum applied average strain.

Fig. 4 .
Fig.4.Applied average strain vs. time up to maximum strains of  = 8.1% and  = 14.3% in experiments 1 and 2, respectively.The tensile loading rate is 100 mm∕min for both experiments, unloading rates differ (10 mm∕min in experiment 1 and 100 mm∕min in experiment 2).Simulations use the exact same strain histories to mimic experiments 1 and 2.

Fig. 5 .
Fig. 5. Relaxation test of a 2D hexagonal truss loaded up to  = 8.1% (left) and  = 14.3% (right) in experiments 1 and 2, respectively.For both experiments, the initial truss configuration before loading (point a ⃝) is compared to the deformed configuration at the indicated times (point b ⃝).Gray trusses represent the simulation results, while black lines were extracted from experiments.Optical images of the two samples are shown in Fig. 3.

Fig. 6 .
Fig. 6.Force vs. time of the 2D hexagonal trusses of experiments (i) 1 and (ii) 2. The viscoelastic stress relaxation of the 3D-printed base material leads to the typical relaxation behavior between b⃝ and c ⃝. Simulations based on the viscoelastic corotational beam element match experimental data to within a normalized root-mean-square error (NRMSE) of (i) ∼ 4.0% and (ii) ∼ 15.48%, measured between a ⃝ and d ⃝.

Fig. 7 .
Fig. 7. Magnified view of the force-time relaxation curve in experiment 1 from Fig. 6(i), comparing the results of FE calculations with different levels of beam refinement to the experimental data (shown as mean plus standard deviation, the latter being indicated by the region shaded in red).The legend reports the normalized root-mean-square error (NRMSE) during stress-relaxation between b ⃝ and c ⃝.

Fig. 8 .
Fig. 8.The twisted cube (made of an octet truss lattice) is shown both in its fully resolved discrete structural version and as an FEM-approximated continuous body, shown in the undeformed configuration with colors illustrating the displacement magnitude.

Fig. 9 .
Fig. 9. Simulated torque vs. time response, comparing the exact structural response of a cube-shaped truss made of 40 × 40 × 40 unit cells to that obtained from the two-scale model used in a 10 × 10 × 10 FE mesh.In both cases, the cube is initially twisted by 3 • ( a ⃝-b ⃝), then held fixed while stress relaxation occurs due to the viscoelastic base material ( b ⃝c ⃝), before being returned to the initial configuration ( c ⃝d ⃝).Results are presented for four types of truss topologies (bitruncated octahedron, octahedron, octet, and cube), all using the same linear viscoelastic base material based on the Prony series of TableC.3.All trusses have the same relative density (cf.Table1).

Fig. 11 .
Fig.11.Simulated rotation vs. time response of a twisted truss cube made of 30 × 30 × 30 unit cells of four different truss topologies (bitruncated octahedron, octahedron, octet, cube) in comparison to the results of the two-scale model used in a 10 × 10 × 10 brick mesh.The cube is initially twisted up to 3.7 • at a rate of 89 • ∕min and held at that twist angle for 2.5 s to allow the material to relax, and then released to return to the initial state through damped vibrations.All trusses have the same linear viscoelastic base material model with the parameters of TableC.3 and the same effective density according to Table1.For comparison, dashed lines illustrate the quasistatic relaxation response of each structure without considering inertial effects.

Fig
Fig. A.12 compares the resulting tensile forces 12(a), shear forces 12(b), and torsional moments 12(c) from both types of simulations for the three loading scenarios with particular focus on the relaxation behavior.In tension (without lateral constraints), the beam deforms homogeneously.In bending and torsion, we fully clamp the cross-sections at both ends for a more realistic deformation mode.As done in Fig.6, we evaluate the NRSME over the entire load history.The normal force shows convincing agreement for small axial strains (NRMSE ∼ 4%) but loses accuracy when the strain level increases (NRMSE > 10%).This effect is not observed under bending and torsion, where the accuracy remains high even up to large strains and large twists.We further study the influence of the slenderness ratio  on the accuracy of the linear viscoelastic corotational beam description.The deformation within each of the three load case scenarios is kept constant, while the beam diameter is increased and the NRSME between the beam response and the FE results is measured.As expected, Fig.A.13 shows a constant error in case of tension (since the applied deformation is homogeneous and the NRSME therefore independent of the diameter).In bending, the NRSME increases with slenderness ratio for both UL and TL settings.Lastly, the error for the twisted beam increases with slenderness ratio  on the accuracy of a twisted beam when using the updated-Lagrangian formulation (while being marginal and constant trend when using the total-Lagrangian setting).The presented comparison of the discrete beam calculation to a fully resolved FEM computation also underlines the efficiency of the proposed model concerning the required degrees of freedom and computational costs.On average, the solid FEM calculation with 6501 linear tetrahedral elements takes ∼ 350 times longer than the same calculation using a single corotational beam with the linear viscoelastic description (one full run is subdivided into 100 load steps).Especially when simulating large truss networks with several thousands of beams, we can achieve almost the same accuracy with the proposed beam model as with an expensive solid FEM calculation.

Fig
Fig.A.12.Net normal force , transverse force , and torsional moment  vs. time of for a cantilever beam of  = 10mm and diameter  = 1mm.All plots compare the response of the presented corotational beam element to a detailed FEM simulation using 6501 linear tetrahedral solid elements in Abaqus(Smith, 2009), based on updated-Lagrangian (UL) and total-Lagrangian (TL) schemes.The deformation is ramped up linearly over 2.5s, then kept constant for 5.0s, until the beam is returned to its initial configuration (analogous to the hexagonal truss study).Tables indicate the NRSME between results from the corotational beam description and the 3D FE simulation.(c) includes a visualization of the fully resolved beam modeled using Abaqus.

Fig. A. 13 .
Fig. A.13. NRSME between the corotational beam formulation and the 3D FE model for an increasing slenderness ratio .The applied deformation is kept constant for all simulations.In Abaqus simulations, the average element size is adjusted such as to prevent ℎ-refinement when thickening the beam.The NRSME between the normal force for tension, the transverse force for bending, and the torsional moment for twisting is measured over the entire course of the simulation, shown in Fig.A.12.FE results were obtained within updated-Lagrangian (UL) and total-Lagrangian (TL) settings.

Fig. B. 14 .
Fig. B.14. Macroscale finite element with the same angular velocity  applied at each node, whose kinetic energy is matched with that of a beam-based RUC, whose nodes are deformed by the same angular velocity.

Fig. D. 16 .
Fig. D.16.Stretching strains and curvature values within individual struts of the hexagonal truss, as obtained from discrete viscoelastic beam simulations of experiment 1 (left) and experiment 2 (right).

Table 1
Geometric parameters for four different 3D truss topologies (where    denotes the side length of the unit cell), which result in the same relative density (i.e., the total volume   occupied by beams within a unit cell of volume    , not considering the multiple-counting of volume at beam junctions) while maintaining a slenderness ratio of  = ∕ ≤ 0.1.All struts within each topology possess the same length  and diameter .
1)(and analogously for the bulk modulus () with parameters   ), which admits the extraction of all dimensionless Prony series contributions   (and   ), as summarized in TableA.2.

Table A .2
Normalized shear and bulk moduli and associated relaxation times of the Prony series calibrated for the TPE base material.