Topology optimization for designing periodic microstructures based on finite strain viscoplasticity

This paper presents a topology optimization framework for designing periodic viscoplastic microstructures under finite deformation. To demonstrate the framework, microstructures with tailored macroscopic mechanical properties, e.g., maximum viscoplastic energy absorption and prescribed zero contraction, are designed. The simulated macroscopic properties are obtained via homogenization wherein the unit cell constitutive model is based on finite strain isotropic hardening viscoplasticity. To solve the coupled equilibrium and constitutive equations, a nested Newton method is used together with an adaptive time-stepping scheme. A well-posed topology optimization problem is formulated by restriction using filtration which is implemented via a periodic version of the Helmholtz partial differential equation filter. The optimization problem is iteratively solved with the method of moving asymptotes, where the path-dependent sensitivities are derived using the adjoint method. The applicability of the framework is demonstrated by optimizing several two-dimensional continuum composites exposed to a wide range of macroscopic strains.


Introduction
Topology optimization is a powerful tool that enables engineers to enhance structural performance. Since the work by Bendsøe and Kikuchi N. (1988), topology optimization has rapidly evolved and been applied to design a variety of physical systems (Deaton and Grandhi 2014;Bendsoe and Sigmund 2013). Specifically, topological design methods based on the inverse homogenization approach by Sigmund (1994)  Center for Design and Optimization, Lawrence Livermore National Laboratory, Livermore, CA, USA composite materials with enhanced properties, by optimizing the material microstructure. Examples of these include linear elastic materials with negative Poisson's ratio (Sigmund 1995), negative thermal expansion (Sigmund and Torquato 1997), and extreme bulk modulus (Huang et al. 2011), and also viscoelastic materials with extreme bulk modulus (Andreasen et al. 2014).
In the design of such architected composite materials, a common simplification in the formulation is to use linear elasticity, which implies that composite materials that exhibit irreversible processes and are exposed to large deformations cannot be designed. While linear assumptions suffice for many composite design optimization applications wherein the microstructure is subjected to moderate macroscopic strain levels, they fail to accurately model the microstructural material response when it is subjected to larger macroscopic strains. To model this microstructural response, nonlinear finite strain theory should be incorporated into the composite topology optimization.
The evaluation of the homogenized material properties becomes demanding when accounting for response nonlinearities, due to the need for iterative nonlinear analyses (Yvonnet et al. 2009). Incorporating nonlinear simulations into topology optimization problems is a formidable task.
Even the efficient parallel programming computations, as presented in Nakshatrala et al. (2013), for the design of nonlinear elastic composites, and the model reduction technique computations, as presented in Fritzen et al. (2016), for designing multiscale elastoviscoplastic structures, remain computationally intensive tasks.
In the works by  and Wallin and Tortorelli (2020), alternative nonlinear homogenization methods have been proposed to design architected materials. These homogenization methods use numerical tensile experiments for calculating the effective uniaxial material response and use topology optimization to achieve desired macroscopic properties. In this way, the computational effort required for the analysis is significantly reduced, since just a few desired properties are targeted.
Of the developed topology optimization methods, few consider structures exhibiting inelastic material response. Contributions incorporating small strain elastoplastic formulations include the work by Maute et al. (1998), where topology optimization is used to optimize structural ductility. Conceptual designs of energy absorbers for crashworthiness, modeled with 2D-beam elements, were established by Pedersen (2004). Protective systems with maximum energy dissipation for continuum structures subject to impact loading have also been designed by Nakshatrala and Tortorelli (2015), where a transient elastoplastic topology optimization formulation is used. Elastoplastic material modeling has further been used in Bogomolny and Amir (2012) for optimizing steel-reinforced concrete structures, in Kato et al. (2015) for optimizing composite structures for maximum energy absorption, and more recently in Wallin et al. (2016) for maximizing the energy absorbtion of structures under finite deformation. Also, in our previous paper , we combine dynamic rate-dependent finite strain effects with topology optimization to design maximal energy absorbing structures. In the design of such energy absorbing structures, e.g., bumpers in automotive engineering, it may also be desirable to control the volume change, since devices that operate in confined spaces will be extremely stiff if the deformation is associated with a volume increase.
To date, most topology optimization frameworks that incorporate nonlinearities focus on macroscopic design; research on composite material design via inverse homogenization that incorporates inelastic constituents is scarce. Examples of such studies include the recent works by Chen et al. (2018) and Alberdi and Khandelwal (2019) which assume infinitesimal deformation. In this paper, we generate optimal topologies of periodic microstructures for maximum energy dissipation with possible prescribed zero contraction constraints under finite deformation. To do this, we incorporate the large strain isotropic hardening viscoplasticity constitutive model, proposed by Simo and Miehe (1992), into a topology optimization framework. Macroscopic nonlinear material properties are evaluated by numerical tensile and shear tests of a single unit cell that is subjected to periodic boundary conditions. Our framework builds upon the works of Ivarsson et al. (2018), as rate-dependent finite strain effects are included, but for simplicity, we neglect inertia. Contrary to our previous work, where structural design was considered, we here study microstructures optimized for maximum energy absorption. Thus, the present work allows each point in a macroscopic structure to be tailored to its specific loading rate, load path, and load magnitude it experiences. In particular, we investigate the influence of nonproportional loading and distortion constraints. Ultimately we quantify the optimization cost and constraint functions via the response from a finite strain viscoplastic boundary value problem which is solved using a nested Newton approach, wherein the momentum balance equation is coupled with local constitutive equations. We discretize the kinematics using a total Lagrangian formulation and use the backward Euler scheme to integrate the constitutive equations.
To formulate a well-posed topology optimization problem, we use restriction wherein we constrain the feature length scale in the volume fraction field by mollifying it via a periodic version of the Helmholtz partial differential equation (PDE) filter (Lazarov and Sigmund 2011). The optimization problem is solved by the method of moving asymptotes (MMA), which requires gradients to form the convex MMA approximation. For path-dependent material response, it is challenging to compute response sensitivities. This is because the sensitivities depend on the response history, as opposed to, e.g., elasticity, where they do not. To compute the gradients, we use the adjoint sensitivity analysis because the design variables in the optimization problem greatly outnumber the constraints. Several studies (Bogomolny and Amir 2012;Wallin et al. 2016;Zhang et al. 2017;Amir 2017;Ivarsson et al. 2018;Zhang and Khandelwal 2019) of topology optimization for path-dependent problems have successfully adopted the adjoint sensitivity recipe presented in Michaleris et al. (1994). This requires the evaluation and storage of the primal response trajectory, and the solution of a terminal-valued adjoint sensitivity problem. We implement this adjoint method to compute the sensitivities.
The remainder of this paper is organized as follows: In Section 2, we present the kinematics and governing mechanical balance equations. Section 3 summarizes the constitutive model and Section 4 presents the numerical solution procedure and the periodic boundary condition implementation. Section 5 describes the regularization and material interpolation schemes. This is followed by the optimization formulation and the sensitivity analysis given in Section 6. The paper is closed by presenting several numerical examples of optimized composite designs and drawing conclusions.

Preliminaries
For nonlinear kinematics, the current, i.e., deformed configuration needs to be distinguished from the reference, i.e., undeformed configuration o . The location of a material particle, which initially at time t = 0 is located at X ∈ o , is described by the smooth motion ϕ such that its position at time t > 0 is x = ϕ(X, t). The local deformation surrounding X is described by the deformation gradient F = ∇ o ϕ = 1 + ∇ o u, where ∇ o is the gradient with respect to the material coordinates, X, and u = ϕ − X is the displacement field. The volumetric change is described by the Jacobian J = det (F ) > 0, where the inequality ensures the mapping ϕ is locally unique. The local deformation is also described by the right Cauchy-Green deformation tensor C = F T F , and the Green-Lagrange strain tensor E = 1 2 (C −1). Since the constitutive model is formulated in an incremental format, we also introduce the spatial velocity gradient l =Ḟ F −1 , whereḞ = dF /dt.
In absence of body forces and inertia, the momentum balance laws, i.e., equilibrium equation is given by where S is the symmetric second Piola-Kirchhoff stress tensor, related to the Cauchy stress σ via S = J F −1 σ F −T . However, when formulating the constitutive model, we will utilize the Kirchhoff stress τ that is connected to S and σ through τ = F SF T = J σ . The boundary conditions associated with (1) are where N o is the outward unit normal vector to the boundary ∂ o of the reference body o and the prescribed tractiont and displacementū are defined on the two complimentary parts of ∂ o , ∂ ot and ∂ ou , respectively. To derive the finite element formulation, we use the principal of virtual work, which can be stated as finding a smooth u that satisfies u =ū on ∂ ou and for all smooth admissible virtual displacements δu that equal zero on ∂ ou . The Lagrangian virtual strain δE, introduced in (3), is defined as δE = 1 2 F T ∇ o δu+ (∇ o δu) T F .

Constitutive model
Following Simo and Miehe (1992) and Ivarsson et al. (2018), the material model is based on isotropic hyperelasticity and isotropic hardening viscoplasticity. For finite strains, the infinitesimal deformation theory's additive split of the strain tensor is replaced by the finite deformation multiplicative split F = F e F vp , where F e and F vp are the elastic and viscoplastic parts, respectively.
The isotropic material response is modeled in terms of the left Cauchy-Green tensor b e = F e (F e ) T , and its isochoric part,b e = (J e ) −2/3 b e , where J e = det (F e ).
The hyperelastic response is governed by the neo-Hookean strain energy function (see, e.g., Simo and Hughes (1998)) where κ and μ are the bulk and shear moduli. Under the hyperelastic assumption, the Kirchhoff stress is obtained via τ = 2(∂w e /∂b e )b e . The rate-dependent evolution laws are obtained by maximizing the regularized dissipation whereẇ vp = τ : − 1 2 Lb e b e−1 is the viscoplastic power per unit volume in the reference configuration, and Lb e = b e − lb e − b e l T is the Lie derivative of b e . In (5), α is the internal variable that defines the isotropic hardening, K defines the drag stress, η ∈ (0, ∞) is the material viscosity, and the Macaulay brackets · are defined such that x = x for x ≥ 0 and x = 0 for x < 0. The von Mises-type dynamic yield function (cf. Ristinmaa and Ottosen 2000) f separates elastic response, f < 0, from viscoplastic response, f ≥ 0. In the above (6), τ dev is the deviatoric part of the Kirchhoff stress and is the current yield stress, where σ y0 is the initial yield stress. The drag stress evolves via the nonlinear exponential hardening model where H, Y ∞ , and δ are the linear hardening modulus, saturation stress, and saturation exponent, respectively.
Maximizing the regularized dissipation (5) with respect to the stresses τ and K provides the evolution equations where n = τ dev /||τ dev || is the normal to the yield surface, and λ = 1 η f is the magnitude of the viscoplastic load increment.
The total amount of absorbed viscoplastic work, W vp , is expressed as where T f is the load duration.

Numerical solution procedure
The boundary value problem defined in (3) is spatially discretized by the finite element method in Section 4.1, and the constitutive equations are discretized with the backward Euler scheme in Section 4.2. Section 4.3 presents the imposition of the periodic boundary conditions and the nested Newton solution procedure. For further details of the finite element discretization, see Belytschko et al. (2000).

Discretization of the equilibrium equations
To obtain the FE-formulation, we approximate the displacement in each finite element e o using the element shape function matrix, N e , according to u(X, t n ) ≈ N e (X)û (n) e , whereû (n) e is the element nodal displacement vector at the discrete time step t n . The virtual displacement is likewise interpolated as δu(X) = N e (X)δû e , i.e., the Galerkin approach is employed. The element Lagrangian virtual strain is expressed in Voigt notation as δE e (X, t n ) = B e X,û (n) e δû e , where the explicit expression of B e is given in Wallin et al. (2018). We also use the compact internal state variable notation w where b e and all other tensors are expressed in Voigt notation.
The residual equilibrium (3) is ultimately discretized as In (12), we introduced the internal force vector F int = F int (û, w) and the external load vector F ext , which are defined by where A is the finite element assembly operator, S = S(û, w) and ξ ∈ [0, 1] is the load parameter. For simplicity, we drop the superscript (·) (n) and assume that all quantities are evaluated at time step t n unless otherwise stated.

Integration of constitutive equations
The residual (12) is evaluated element-wise using Gauss integration. As such, the constitutive equations need to be integrated in time at each Gauss point of each element in the mesh. To this end, the viscoplastic evolution (9) and (10) are combined with (6) and (7) and discretized in time using the fully implicit Euler scheme. Ultimately we obtain one vector and one scalar local Gauss point-i residual equation.

Periodic boundary conditions
We now evaluate the homogenized material properties of the unit cell microstructure, i.e., composite. As part of this process, Bloch-type periodic boundary conditions are applied to a single unit cell, which implies that the displacement and traction fields on opposite boundaries are periodic, modulo a uniform term and antiperiodic, respectively. Without loss of generality, we restrict ourselves to a two-dimensional (2D) formulation. The periodicity is based on master-slave elimination wherein the displacement on the boundaries ∂ Fig. 1a) are separated into their master and slave displacements, u + kl and u − kl . Here the first subscript k represents the displacement component, and the second subscript l represents the boundary, e.g., u + xy is the u x displacement component over the surface The differences between master and slave displacement components on opposing faces are, e.g., whereū xx is the uniform axial displacement difference enforced across the boundaries ∂ x+ o and ∂ x− o , cf. Fig. 1b. The slave degrees of freedom are explicitly eliminated from the nonlinear system (12) by expressing the nodal displacement vector aŝ u = Tǔ +ū xx I xx +ū xy I xy +ū yx I yx +ū yy I yy , whereǔ contains the nodal displacement components of the master and the interior (belonging to o \ ∂ o ) degrees of freedom, and the constant matrix T defines the connection between the master and slave degrees of freedom. In (20), the vectors I xx and I yx have unit entries respectively for the u x and u y degrees of freedoms corresponding to the nodes on the boundary ∂ x+ o and zero entries elsewhere. I xy and I yy are similarly defined.
For later purposes, (20) is decomposed into a dependent and independent part. To simplify the derivations, only the x uniaxial load case of (19) is considered. To do this, we assign the value ofū xx and decompose (20) aŝ where P = T I xy I yx I yy , andû * = ǔ Tū xyūyxūyy T contains the dependent displacements. To accommodate the addedū xy ,ū yx , andū yy degrees of freedom, the residuals (12) and (17) are supplemented with the null traction conditions such that the boundaries ∂ y± o are macroscopically traction free and ∂ x± o has zero shear traction, i.e., uniaxial conditions are mimicked.
The internal virtual work δW int = F T int δû together with periodicity of (21) and the null traction of (22) imply that where, e.g., R xy = I T xy F int . Equation 23 is used to define the compact residual a b Fig. 1 Illustration of the periodic boundary conditions. a Undeformed unit cell divided into four boundaries and b deformed (dotted lines) unit cell with prescribed and unknown gap displacements, and anti-periodic traction t The coupled nonlinear system R * = 0 and C i = 0 of (17) and (24) is solved using two nested Newton loops, wherein linearizations around the current (û , w) iterate are performed using a Schur's complement approach Simo and Miehe 1992). At each iteration, an inner Newton loop is first completed to solve C i (û (n) ,û (n−1) , w (n) , w (n−1) ) = 0 for w (n) for the given values ofû (n) ,û (n−1) , and w (n−1) . The outer loop then commences whence the linearized system is solved for dû * , where K * T = P T K T P follows from (21) and (24), in which K T is the consistent tangent stiffness matrix defined by where C contains the collection of the N gp Gauss point residuals C i , i = 1, 2, . . . , N gp . The use of the consistent tangent is crucial for preserving the quadratic rate of asymptotic convergence of the Newton iteration procedure (Simo and Taylor 1985). The successive updates of w (n) and u * (n) continue until convergence is achieved. The time is then incremented and the procedure begins anew with the initial guess for the displacementsû * (n) 0 . A common choice for the initial displacement guess is the previous converged solution, i.e.,û * (n) 0 =û * (n−1) . Based on our experience, this particular choice is not efficient and therefore, we use the initial guesŝ i.e.,û * (n) 0 is a scaling of the previous converged solution, which leads tô The first n = 1 time step uses the initial guessû Remark 1: In the numerical examples, we enforce orthotropic domain symmetry of the unit cell. For the x uniaxial case, this a priori ensuresū xy =ū yx = 0 and R xy = R yx = 0. As such, we only enforce the I T yy F ext = 0 reaction condition to evaluateū yy .
Remark 2: For a pure shear load case, we prescribeū xy > 0,ū yx = 0 and computeū xx andū yy to enforce I T xx F ext = 0 and I T yy F ext = 0.

Regularization and material interpolation
In this topology optimization study, we use restriction to formulate a well-posed problem. The continuous nondimensional volume fraction field c(X) describes the amount of material for reference point X, where locations filled with material are defined by c(X) = 1 and void locations are identified by c(X) = 0. The continuous volume fraction field c is a function of the design field φ(X) ∈ [0, 1] which is discretized to be piecewise uniform over the elements e o such that φ(X) = φ e for X ∈ e o . In this relationship, c is obtained by smoothing φ so as to define a well-posed topology optimization problem. The smoothing uses the Helmholtz partial differential equation filtering technique, cf., e.g., Lazarov and Sigmund (2011), wherein where l o controls the length scale, i.e., smoothing and Δ o is the Laplacian with respect to the material coordinates X. The filtering technique has previously been used with homogeneous Neumann boundary conditions, ∇ o c ·N o = 0 to conserve volume. In this work, periodicity is enforced on c, which implies that ∇ o c · N o is antiperiodic and hence we also conserve volume. We solve (29) using a standard finite element formulation. The filtered continuous volume fraction field c and the weight function δc are interpolated using the element shape function vector, i.e., c(X) = c N e (X)c e within finite element e o , and the resulting linear system is solved for c, i.e., the nodal values of c. In (30), . .] T is the design variable vector. Since the finite element discretization is fixed, the matrices K f and F f are computed once. In order to enforce periodicity, similar to (20), the constraint is included where the connectivity matrix T f connects the master nodes, stored in c * , to the slave nodes. Inserting (31) into (30), and premultiplying by T T f , yields the resulting linear system where Upon computing c * from (32), the volume fraction field is obtained from (31) and the finite element interpolation c(X) = c N e (X)c e .
To approach a distinct solid-void design, we use a Heaviside thresholding technique (Guest et al. 2004;) and material penalization. First, the continuous filtered volume fraction c is thresholded using the smoothed Heaviside step functionH , to define the physical volume fraction field The parameters β H and ω are such that lim β H →∞H (c) = u s (c − ω), where u s is the unit step function. Second, the material parameters κ, μ, H, σ y0 , and Y ∞ for the constituent material are penalized via the RAMP material interpolation scheme (Stolpe and Svanberg 2001). Using the physical volume fractionĉ of (33), the material interpolation function for each parameter χ is defined as where the parameter q controls the level of penalization, χ 0 is the nominal material parameter value, and δ 0 = 10 −4 is a small residual value used in this ersatz material model to avoid numerical problems which otherwise arise if, e.g., the material stiffness is zero.
To summarize, the design variables φ e are first filtered by solving (32) for c * , then the thresholded physical periodic volume fraction fieldĉ is computed via (31) and (33), and finally the penalized material parameters are evaluated through the application of (34) for χ 0 = {κ, μ, H, σ y0 , Y ∞ }.

Topology optimization
The objective considered in this topology optimization is to design a material microstructure with maximum viscoplastic energy absorption, W vp , cf. (18), subject to displacement constraints. One displacement constraint enforced during uniaxial loadingū xx prescribes the transverse macroscopic displacement componentū yy . This Poisson's ratio-type constraint is imposed by requiring and using the implicit Euler scheme to discretize (35). As seen above, gū yy < δū ensures the root-mean-square deviation (RMSD) ofū yy from the prescribed target trajectoryū * yy is less than a small value δū taken as 3 · 10 −3ū xx . Similarly,ū xx can be constrained for a prescribed uniaxial loadingū yy .

Optimization formulation
The topology optimization problem is to find a design φ by solving O : whereV is the maximum allowable volume in the reference configuration and the mappinĝ which follows from (31)-(33) and the interpolation c = c N e c e , within each element. In the above reduced space optimization problem, we treat the equilibrium (24) and local constitutive (17) residuals as implicit constraints, i.e., in (36) we require whereû =û(û * (f (φ))) and w = w(f (φ)).
Remark 3: The optimization problem (36) allows for maximization of viscoplastic energy absorption subject to a maximum volume constraint by letting δū → ∞.

Sensitivity analysis
The optimization problem (36) is solved with the MMA scheme, cf. Svanberg (1987), using the standard MMA parameters as described in Svanberg (2002). The algorithm uses gradients of the objective and constraint functions to form a convex approximation of the optimization problem. We therefore perform a sensitivity analysis in each optimization iteration, and as the number of design variables greatly outnumber the constraints, we use the adjoint method presented in Michaleris et al. (1994). The sensitivity analyses for the objective function W vp and the constraint gū yy are performed in a similar manner by considering the generalized response function The sensitivities with respect to the design variable φ are obtained by multiple applications of the chain rule, i.e., cf. (37), where ∂c ∂c * = T f and ∂c * ∂φ = K * −1 f T T f F f which follows from (31) and (32), respectively. In (41), we introduced D(·) D( ) to denote the total differentiation of (·) with respect to ( ). Differentiation of with respect to the filtered nodal volume fractions c gives where we use the equality ∂û (n) /∂û * (n) = P . The implicitly defined sensitivities Dû * (n) /Dc and Dw (n) /Dc are annihilated by augmenting a series of adjoint equations to the functional . To this end, the adjoint vectors λ * (n) and γ (n) are used to define the equivalent functional where we emphasize that the augmented functional˜ equals , since (38) and (39) hold. Differentiating the augmented function (43) with respect to c and applying the chain rule yields The above is divided into an explicit, D˜ E /Dc, and implicit, D˜ I /Dc, part as D˜ /Dc = D˜ E /Dc + D˜ I /Dc, where the explicit part is and (the transpose of) the implicit part is The implicitly defined derivatives Dû * (N) /Dc and Dw (N) /Dc in (46) are annihilated by making judicious choices of the adjoint vectors. We first require that λ (1) and γ (1) satisfy Or, by inserting (48) into (47) and some rearrangement, satisfy where K * (N) T = P T K (N) T P is the tangent operator at the terminal time step t N . Once λ * (1) is obtained from (49), the adjoint variable γ (1) is obtained from (50). In the remaining time steps of the adjoint analysis, we determine the λ * (n) and γ (n) to annihilate Dw (N −n+1) Dc and Dû * (N −n+1) Dc for n = 2, . . . , N, by solving Alternatively, substituting (52) into (51), rearranging and using (26), we solve is the adjoint load vector at step t N−n+1 . Equation 53 is solved for λ * (n) and subsequently used in (54) to obtain γ (n) . Once the adjoint vectors λ * (n) and γ (n) are computed, the sensitivity D /Dc = D˜ E /Dc is obtained from (45), and D /Dφ is computed via (41).
The last two terms of the sensitivity D˜ E /Dc of (45) are computed in the same manner for both the constraint = gū yy and for the objective function = W vp . Using (13), (24), ∂ĉ/∂c =H and ∂c/∂c = c N T e within e o , the first term of (45) is expressed as (59) where λ e contains the e o element components of λ = P λ * , and ∂S (n) which follows from (4), τ = 2(∂w e /∂b e )b e , τ dev = μb e,dev and (34). The last term of the sensitivity (45) is

Numerical examples
In this section, the presented topology optimization framework is applied to design several solid-void composite materials for maximal energy absorption. The material parameters for the full material are those of steel and are presented in Ivarsson et al. (2018).
The primal viscoplasticity problem is solved with an adaptive time stepping scheme such that where N newt is the number of Newton iterations that is required to attain convergence at the current time step t n , cf., e.g., Borgqvist and Wallin (2013), and the parameters f iter and v are chosen as f iter = 4 and v = 0.4. The convergence criterion is such that |R * T Δû * | < 10 −10 Nmm and diverging time steps are restarted with Δt n+1 ← Δt n+1 /2. The initial time increment is Δt 1 = T f /500 and the time increments used to perform the sensitivity analysis are consistent with those used to solve the primal problem. Although the analysis is for two dimensions, we use fully integrated 8-node tri-linear brick elements to discretize the 1 mm × 1 mm × 0.01 mm unit cell with 80 × 80 × 1 elements. Plane strain boundary conditions are enforced on the z−displacement components. Because of this plane strain condition, the z-faces of the designs are necessarily periodic.
The Heaviside threshold parameter is ω = 0.5 and a continuation strategy is used for β H wherein its initial value β o = 10 −10 is used for the first 80 MMA optimization iterations, whereafter it is increased by increments of 1 every n β design updates until the maximum value of β H = 10 is reached. In the first two examples n β = 10, and for the optimizations including displacement constraints, the continuation is slower with n β = 20. The RAMP penalization parameter is q = 10 for κ, μ, H, and Y ∞ . To avoid large viscoplastic strain in low-density elements, the penalization parameter for the initial flow stress σ y0 is chosen such that q σ y0 = 2 < q. The maximum volume fraction isV = 0.4 and the filter length scale parameter is l o = 1.25/80 mm.
In this study, 4-fold orthotropic domain symmetry of the design is imposed (case a), cf. Fig. 2a. This greatly reduces the number of design variables. Moreover, in the analysis, e.g., for uniaxial loading (1) the periodic boundary conditions are replaced by the usual domain boundary conditions over the quarter domain model and (2) zero shear forces are ensured (see, e.g., Wallin and Tortorelli (2020)), i.e., we necessarily have R xy = R yx = 0. Thus, we save computational effort as we need not consider the I T xy F ext = 0 and I T yx F ext = 0 equations. We also impose cubic 8-fold domain symmetry (case b), cf. Fig. 2b. This further reduces the number of design variables; however, to avoid the imposition of complicated boundary constraints, we analyze this as if it has 4-fold symmetry. The latter case b ensures microstructures that perform identically when subjected to uniaxial loadsū xx andū yy of same magnitude and thus, only one uniaxial load case is simulated. The initial design must be nonuniform to avoid initially uniform sensitivities that inhibit the optimization (Alberdi and Khandelwal 2019;Träff et al. 2019). The optimization algorithm is therefore initiated such that all design variables are set to φ e = 0.4, followed by a perturbation δ. We consider two different perturbations; (a) a value δ 4 is added to the four center elements, i.e., φ e = 0.4 + δ 4 , and (b) a different random number perturbation is added to all elements, 1 such that φ e = 0.4 + δ rnd × (rand(0, 1) − 0.5).

Constitutive model
Prior to the optimization, we demonstrate the rate effects associated with the constitutive model. The demonstration is performed using a single element model of volume fraction field φ e = 1, by simulating the uniaxial response for u xx = 0.05 mm with symmetry boundary conditions of Fig. 2a. Figure 3 shows the macroscopic force component I T xx F int = F xx , normalized by a factor of 1/(σ y0 A o ), where A o = 0.01 mm 2 is the area of the loaded unit cell x-face in the reference configuration, versus loadū xx for the three different macroscopic engineering strain ratesε xx = 10 −15 s −1 ,ε xx = 100 s −1 , andε xx = 500 s −1 . The figure clearly shows the influence of the strain rate effects on the microstructural response.

Test problem: two simple load cases
We start by designing materials subjected to a slow loading rateε = 10 −15 s −1 , without any displacement constraints. We do not enforce domain symmetry in this initial example and we consider two load cases: (a) a uniaxial tensile load case whereū xx = 0.01 mm is prescribed and we solve for u xy ,ū yx ,ū yy , and (b) a pure shear load case whereū xy = 0.01 mm andū yx = 0 mm are prescribed and we solve for u xx andū yy . Figure 4 shows two optimized microstructural designs, both initiated with δ 4 = 0.6, i.e., the four center elements are assigned the value 1. The yellow and blue regions represent the base steel material and void, respectively, and we emphasize that the physical volume fraction fieldĉ is plotted. As expected, the two structures dissipating maximum energy consist of bars aligned with the maximum principal stress axes. Clearly, these designs are periodic, and although not explicitly enforced, they respectively exhibit 4-fold and 8-fold domain symmetries.

Biaxial tensile loads
In many practical applications, it is of great importance that structures can withstand more complex load cases than those of uniaxial tension and pure shear. Therefore, we consider a biaxial loading scenario,ū xx =ū yy = 0.01 mm again withε = 10 −15 s −1 , cf. Fig. 5. In this study, we use 4-fold domain symmetry and the corresponding symmetry boundary conditions of Fig. 2a. We present optimized microstructures obtained using three different initial designs. Figure 6 a and b show two designs optimized with δ 4 = 0.1 and δ 4 = −0.1, respectively. The topology of Design 6a consists of perpendicular bars aligned in the loading directions, whereas Design 6b has a more complex microstructure with both rectangular and elliptical material regions. This difference is reflected in the performance as Design 6a absorbs 1.0 % more viscoplastic energy than Design 6b, (cf. Fig. 7), i.e., Design 6b is a local maxima. We also note that Design 6a absorbs approximately 123 % more specific viscoplastic energy than a fully solid unit cell in whichĉ = 1.
Next, the optimization is performed using a random initial 8-fold symmetric perturbation δ rnd = 0.1. The optimized design, cf. Fig. 6c, consists of an octagonal frame with connected bars aligned in the loading directions. The energy absorption of this design is also slightly (1.1 %) lower than that of Design 6a.
Designs 6a-c are alternate solutions to the optimization problem. And as shown in Fig. 7, the resulting energy absorption capabilities are not equal. In the upper part of this figure, we also see that the convergence towards a distinct solid-void design is slower for Designs 6b-c as compared with Design 6a. It is also seen that only minor design changes are observed when the thresholding parameter β H is increased after optimization iteration 80, but the viscoplastic work increases as the intermediate values ofĉ are thresholded closer towards 0 and 1.
We conclude that that the initial design choice is crucial, and that several local maxima to this optimization problem exist.

Varying biaxial load paths
Next, we demonstrate the importance of the path-dependent material model for material design under biaxial load cases. We consider the same final load stateū xx =ū yy = 0.01 mm at t = T f as for the previous optimization case, but vary the load path. Specifically, we use the same initial design δ 4 = 0.1 and compare the optimized Design 6a, obtained with the loading conditionsū xx (t) =ū yy (t) = t/T f · 0.01 mm, to several other loading scenarios withū xx (t) =ū yy (t) for 0 < t < T f . Figure 8 shows the resulting optimized designs with corresponding prescribed trajectoriesū xx (t) andū yy (t) which are split into two piecewise linear trajectories where the parameter a values corresponding to each design appear in Table 1. As seen in Fig. 8, different loading trajectories generate different optimized microstructures. This study clearly demonstrates the importance of including the

Displacement constraints
In the following examples, we design microstructures with nearly zero transverse contraction when subjected to uniaxial tensile loading, e.g., for theū xx loading, the constraint a b c Fig. 9 Undeformed optimized designs (left) and array of 3 × 3 unit cells (right) under the loadū xx = 0.01 mm gū yy is included in the optimizations with the target trajectoryū * yy = 0 ∀ t n ∈ [0, T f ]. To avoid constraining both the transverse and z-component displacements, plane stress conditions are used, i.e., periodicity is not enforced on the z-faces. An initial design with δ rnd = 0.002 is used. For the cases where viscous effects are studied, i.e., whenε = 10 −15 s −1 , the loading rate is ramped up froṁ ε = 0.5 s −1 to the specified final value during the first two optimization iterations. This avoids the possibility that the initial design response is nearly elastic which we found otherwise hinders the optimization.

Influence of strain rate and load magnitude
In this example, we design material microstructures that perform equally under two uniaxial loadsū xx andū yy by imposing the 8-fold domain symmetry illustrated in Fig. 2b. Thus, the computational effort is eased as only the single uniaxial load caseū xx needs be simulated. The effect of the load magnitude and loading rate on the microstructural topology is demonstrated by performing five different optimizations (Figs. 9 and 11). In the first three Designs 9a-9c are optimized under the moderate loadsū xx = 0.01 mm, with slowε xx = 10 −15 s −1 , intermediateε xx = 100 s −1 , and fastε xx = 500 s −1 load rates (cf. Section 7.1 and Fig. 3).
We compare the performance of these three designs by exposing them to the same load magnitudeū xx = 0.01 mm and fast loading rateε xx = 500 s −1 , and measure their corresponding transverse contraction −ū yy , cf. Fig. 10. In this figure we see that, although Designs 9a and 9c have a similar microstructure, Design 9a violates the contraction constraint (35), cf. the solid-triangle line in Fig. 10. Design 9b, corresponding to the solid-circle line of Fig. 10, also contracts more than the allowed tolerance δū = 3 · 10 −3ū xx = 3 · 10 −5 mm. Figure 10 also shows that the contraction of all three designs is linearly increasing for loadsū xx less than approximately 0.002 mm but is nonlinear for higher loadsū xx > 0.002 mm. Clearly, the linear/nonlinear contraction versus load is a consequence of elastic and viscoplastic material response.
The last two optimizations are performed for a higher load u xx = 0.05 mm, at slowε xx = 10 −15 s −1 and fastε xx = 500 s −1 loading rates. The optimized designs are shown in Fig. 11. Although no obvious similarities can be seen from the optimized unit cells (cf. the left column of Fig. 11), these two optimizations generate two microstructures with similar topology as seen in the right column of Fig. 11. Both Designs 11a and 11b satisfy the displacement constraint as shown by the center solid-circle and solid-triangle lines of Fig. 12a.
We further analyze the performance of the two Fig. 11 designs by plotting their corresponding macroscopic forcedisplacement curves, cf. Fig. 12b. The two center solidcircle and solid-triangle lines correspond to the macroscopic force-displacement responses of Designs 11a and 11b, respectively. As expected, Design 11b which is optimized a b Fig. 11 Undeformed optimized designs (left) and array of 3 × 3 unit cells (right) under the loadū xx = 0.05 mm for a higher loading rate is exposed to higher macroscopic forces.
Next, we simulate the responses of Designs 11a and 11b at the loadū xx = 0.05 mm, but with interchanged loading rates, i.e., Designs 11a and 11b are subjected to the load rateṡ ε = 500 s −1 andε = 10 −15 s −1 , respectively. Results of these simulations are displayed with red lines in Fig. 12 a and b.
The bottom two curves of Fig. 12b present the forcedisplacement responses of Designs 11a and 11b when subjected to the slow loading rate. As expected, Design 11b which is not optimized for this load rate performs worse than Design 11a; it supports less force and violates the displacement constraint, in fact it has negative transverse contraction. Surprisingly, Design 11a supports more force a b Fig. 12 Performance of Designs 11a and 11b. a contraction −ū yy versus loadū xx and b normalized macroscopic force F xx /(σ y0 A o ) versus loadū xx . The red lines with filled markers correspond to responses obtained from simulations with interchanged loading rates and hardens more than Design 11b when it is exposed to a fast loading rate. However, it violates the displacement constraint. These simulations show that more contraction leads to higher viscoplastic energy absorption. They also demonstrate that it is necessary to include rate effects in the optimization of viscoplastic microstructures.
We also compare the performance of Designs 9a and 11a by exposing them to the same loadū xx = 0.05 mm and slowε xx = 10 −15 s −1 load rate. Figure 13 a shows that Design 9a which is optimized for the moderate load level u xx = 0.01 mm performs slightly better than Design 11a which is optimized for theū xx = 0.05 mm load, i.e., in this example, the energy absorption capability is not improved when optimizing for high load magnitudes. a b Fig. 13 a Normalized viscoplastic energy absorption versus loadū xx and b −ū yy versusū xx , for Designs 11a and 9a. Design 9a, optimized forū xx = 0.01 mm, absorbs 2.7 % more viscoplastic energy but also has 6.5 times higher RMSD ofū yy from the targetū * yy = 0 during theū xx = 0.05 mm load than Design 11a (which is optimized for u xx = 0.05 mm) However, Design 11a fails to satisfy the displacement constraint when subjected to the larger load, cf. Fig. 13b, which shows strain-dependent contraction in theū xx ∈ [0.01, 0.05] mm interval, in contrast to Design 11a that exhibits nearly zero contraction throughout the load interval. Similar to the results presented in Fig. 12, this comparison shows that higher contraction leads to increased viscoplastic energy absorption.

Multiple uniaxial tensile loads: influence of load magnitude
In this final example, we optimize the unit cell under multiple load cases. Specifically, both longitudinal and transversal load cases are simulated using 4-fold domain symmetry with symmetry boundary conditions (cf. Fig. 2a). We include both the gū yy and gū xx constraints in the optimization.
The objective is to maximize the geometric mean of the viscoplastic energy absorption, (W  Fig. 14. The microstructural Design 14a, optimized for the moderate load magnitude 0.01 mm, consists of a rectangular pattern of material with slight offsets to satisfy the near zero contraction constraints, whereas Designs 14b and 14c, optimized for the higher load magnitudes 0.05 mm and 0.10 mm, respectively, consist of more complex topologies. All three designs satisfy the displacement constraints.
The macroscopic force-displacement responses of Designs 14a-14c, plotted in Fig. 15, show that Design  14a supports higher forces and thus can dissipate more viscoplastic energy than Designs 14b and 14c in the load interval [0, 0.01] mm, for both theū xx andū yy loads. It is also seen that Design 14b supports higher forces than Design 14c in the load interval [0, 0.05] mm for both loads.
Finally, we compare Designs 14a and 14c, optimized for the load magnitudes 0.01 mm and 0.10 mm, respectively, by exposing them to the same high load interval [0, 0.10] mm. Figure 16 shows pronounced difference in the contraction versus load, e.g., −ū yy versusū xx between Designs 14a and 14c. Indeed, Design 14a has nearly zero contraction in the [0, 0.01] mm load interval for which it was designed but shows a considerable strain-dependent contraction in the (0.01, 0.10] mm load interval. Figure 17 shows the normalized macroscopic forcedisplacement responses for Designs 14a and 14c. Here it is seen that Design 14c, optimized for the high load magnitude 0.10 mm, softens for loadsū xx ,ū yy > 0.08 mm, whereas Designs 14a exhibits macroscopic hardening. Similar to the previous studies, this example demonstrates that higher contraction (cf. Fig. 16) leads to more supportive structures (cf. Fig. 17) and thus higher viscoplastic energy absorption.

Conclusions
We have established a topology optimization framework for designing periodic microstructural composite materials for maximum viscoplastic energy absorption. Constraints have also been imposed on designs such that, under uniaxial tensile loads, their macroscopic transverse contraction is close to zero. The optimization problems are solved using the MMA scheme and sensitivities required to drive the algorithm are obtained using the adjoint method. The effects of strain rate, load magnitude, and path-dependency are shown by optimizing designs under several complex macroscopic load cases. It is clear that these effects are of great importance on the optimized microstructural design. The present approach allows each material point in a macroscopic structure to be tailored to the loading rate, load path, and load magnitude it experiences.
The unit cell that characterizes the design domain has been discretized with 3D finite elements. This formulation enables a smooth transition to large-scale 3D problems by relaxing of the plane strain and plane stress assumptions. In future work, we will apply parallel programming to the framework to efficiently solve such large-scale problems.

Appendix: Sensitivity verification
To verify correct implementation of the adjoint sensitivity analysis presented in Section 6.2, we compare the sensitivities obtained by the adjoint method with numerical sensitivities obtained by the central difference method (CDM). We use 4-fold domain symmetry and the corresponding symmetry boundary conditions of Fig. 2a to discretize a quadrant of the unit cell by 5 × 5 × 1 tri-linear brick elements, and subject it to biaxial loading conditions (cf. Fig. 5) with load path illustrated in Fig. 8e. To ensure large deformations and viscoplastic response, the final prescribed displacementsū xx andū yy are increased to 0.05 mm and the load duration is decreased to T f = 0.01 s; this is done by scaling the Fig. 8 load path. Each design variable value is given by the random number generator through φ e = rand(0, 1), and for the PDE filtering, we set l o = 0.25 mm. The adaptive time-stepping is omitted for this analysis that instead uses 40 equally spaced load increments. Figure 18 shows that the sensitivities obtained by the adjoint (solid, blue curve) and the central difference method (dashed, red curve) are almost equal; the relative error between the two methods is smaller than the perturbation size Δh = 10 −6 as shown by the dotted line; hence, the derivation and implementation of the adjoint sensitivities are correct.