Thermodynamically consistent concurrent material and structure optimization of elastoplastic multiphase hierarchical systems

The concept of concurrent material and structure optimization aims at alleviating the computational discovery of optimum microstructure configurations in multiphase hierarchical systems, whose macroscale behavior is governed by their microstructure composition that can evolve over multiple length scales from a few micrometers to centimeters. It is based on the split of the multiscale optimization problem into two nested sub-problems, one at the macroscale (structure) and the other at the microscales (material). In this paper, we establish a novel formulation of concurrent material and structure optimization for multiphase hierarchical systems with elastoplastic constituents at the material scales. Exploiting the thermomechanical foundations of elastoplasticity, we reformulate the material optimization problem based on the maximum plastic dissipation principle such that it assumes the format of an elastoplastic constitutive law and can be efficiently solved via modified return mapping algorithms. We integrate continuum micromechanics based estimates of the stiffness and the yield criterion into the formulation, which opens the door to a computationally feasible treatment of the material optimization problem. To demonstrate the accuracy and robustness of our framework, we define new benchmark tests with several material scales that, for the first time, become computationally feasible. We argue that our formulation naturally extends to multiscale optimization under further path-dependent effects such as viscoplasticity or multiscale fracture and damage.


Introduction
Multiphase hierarchical systems apply the concept of microheterogeneity repetitively across a hierarchy of well-separated length scales: composite microstructures at a smaller scale form the base constituents for new microstructures at the next larger scale.This principle constitutes the backbone of virtually all biological materials, enabling them to combine various functional properties at different length scales with favorable mechanical properties at the macroscale through evolutionary mechanisms (Wegst et al. 2015;Zheng et al. 2014;Fratzl and Weinkamer 2007;Ritchie et al. 2009;Bhushan 2009;Egan et al. 2015).In other words, biological materials adapt their form (or shape/structure) against the dynamic external environment and improve the microstructure architecture, fulfilling the local needs imposed by physiological, phylogenetic, and reproductive constraints (Wolff 1986;Gibson 2012;Gao et al. 2003).A rational understanding of microstructure interdependencies across hierarchical scales on the macroscale properties helps pave the way forward to many engineering applications involving biological materials such as the genetic tailoring of crops (Brulé et al. 2016;McCann et al. 2014), bone remodeling (Rodrigues et al. 1999;Blanchard et al. 2016), and the fabrication of bioinspired engineering materials (Wegst et al. 2015;Holstov et al. 2015).
Multiscale modeling of hierarchical materials in conjunction with structural optimization methods constitutes a promising pathway to elucidate optimum microstructure configurations in multiphase hierarchical systems.In this context, recently developed concurrent multiscale analysis 195 Page 2 of 31 and topology optimization methods (Xia andBreitkopf 2014, 2015;Rodrigues et al. 2002;Coelho et al. 2008;Nakshatrala et al. 2013;Da et al. 2017;Zhang et al. 2018) naturally fit to the dual optimization of structure (form) and material (microstructure architecture).The idea is to decompose the multiscale problem into two nested sub-problems, one at the macroscale (structure) and the other at the microscale (material).At each macroscale material point, the microscale sub-problem provides a locally optimal material response and can be interpreted as the reformulation of a material constitutive law for the macroscale structure optimization problem.The material optimization problem is typically solved within an FE 2 type computational homog- enization framework (Feyel and Chaboche 2000;Fish 2013) at each macroscale Gauss point.Other approaches not based on scale separation or periodicity such as Nguyen and Schillinger (2019) also seem possible.We note that throughout the article, we will use the term multiphase hierarchical system for the combined representation of the multiphase hierarchical material and the macrostructure domain that habitats it.
The base constituents in multiphase hierarchical systems often exhibit elastoplastic material properties resulting in a path-dependent macroscale mechanical response and dissipation-driven self-adapting mechanisms.In the case of path-dependent problems, efficient methods for the computational optimization of multiphase hierarchical systems are still in its infancy (Da 2019).In the context of fracture resistance and damage, recent contributions have proposed structure optimization methods optimizing the inclusion characteristics in matrix-inclusion type multiphase materials for path-dependent objective functions (Xia et al. 2018;Li et al. 2021;Kato 2010;Kato and Ramm 2013;Hilchenbach and Ramm 2015).The corresponding optimization formulations, however, remain in the format of a monoscale design, where all the morphological design parameters of the material are represented at the structure scale.In this case, the structure scale discretization is dictated by the smallest length scale of constituents, making these methods computationally prohibitive for multiphase hierarchical systems with several well-separated length scales.To the best of our knowledge, no work has been reported so far in the literature on a decomposed concurrent material and structure optimization formulation for path-dependent problems involving elastoplastic multiphase hierarchical systems.
The first major roadblock in the context of elastoplastic behavior across hierarchical scales is the non-trivial problem decomposition into material and structure optimization subproblems.The current state of concurrent material and structure optimization methods focuses only on end-compliance type optimization problems with an overall linear elastic response at both the material and structure levels.Interested readers are referred to a recent review article by Wu et al. (2021) that extensively covers existing approaches for designing hierarchical structures for linear end-compliance minimization problems.The variational structure of the displacement-based formulation of end-compliance optimization corresponds to a saddle point problem with respect to the admissible set for design variables and the space of kinematically admissible displacements (Lipton 1994).With pointwise definitions of material design variables, the saddle point property enables a natural decomposition into the material and structure optimization subproblems (Jog et al. 1994).The equivalent variational structure for combined analysis and optimization of path-dependent problems that consider the complete deformation process is non-trivial and yet to be investigated, which is also one of the conclusions in the current review article (Wu et al. 2021).
The second critical roadblock is the computational cost for multiscale analysis through computational homogenization that for multiphase hierarchical systems grows exponentially with each scale characterization (Yuan and Fish 2009;Le et al. 2015;Liu et al. 2016;Bessa et al. 2017).Adding the topology optimization at the structure level results in even higher computational cost, since it requires solving many multiscale problems for different realizations of the structure topology during a typical optimization algorithm.This drawback limits existing approaches to small two-scales problems, even in the simplest case of hierarchical materials with linear elastic constituents.Continuum micromechanics provides a rigorous framework to derive analytical estimates of macroscale elastoplastic properties (Zaoui 2002;Suquet 2014;Morin et al. 2017) and has been successfully applied to describe many multiphase hierarchical systems such as plant, wood, bone, and cementitious materials (Gangwar and Schillinger 2019;Gangwar et al. 2021;Hofstetter et al. 2005;Hellmich et al. 2004;Fritsch et al. 2009;Pichler and Hellmich 2011).In the context of concurrent material and structure optimization, we recently integrated continuum micromechanics based homogenized estimates in end-compliance optimization problems, which rendered our framework computationally tractable for multiphase hierarchical systems with several material length scales (Gangwar and Schillinger 2021).
In this article, we establish, for the first time, a thermodynamically consistent formulation of concurrent material and structure optimization, including suitable sub-problem formulations for multiphase hierarchical systems with elastoplastic constituents at the material scales.The structure optimization problem addresses the macroscale density distribution, while the material optimization problem at each material point seeks the optimized macroscale response with respect to microscale (material) design variables.In particular, we reformulate the material optimization problem based on the maximum plastic dissipation principle such that it assumes the format of an elastoplastic constitutive law that can be efficiently solved via modified return mapping algorithms.To focus on the key concepts of the decomposed formulation, we make a few assumptions on the macroscale behavior of the inelastic hierarchical materials including isothermal processes, the existence of associative flow rule, and rate-independent ideal elastoplastic response.
We express the homogenized stiffness and yield criterion as a function of material design variables within a continuum micromechanics framework, which enables the computationally tractable treatment of our optimization formulation.In particular, we focus on a quadratic stress average micromechanical approach for estimating homogenized yield criterion, which has been successfully used in modeling the limit strength of a broad range of hierarchical materials such as metal-matrix composite, cement-mortar, wood, and crop stems (Suquet 1997;Pichler and Hellmich 2011;Hofstetter et al. 2008;Gangwar et al. 2021).
Our article is organized as follows: In Sect.2, we briefly review the relevant thermomechanical principles of elastoplasticity along with multiscaling concepts in continuum micromechanics, which form the basis of our further developments.In Sect.3, we formulate the path-dependent stiffness maximization problem, decomposing material and structure optimization sub-problems for elastoplastic multiphase hierarchical systems.We then describe its discretization within the framework of the finite element method.In Sect.4, we develop an algorithmic procedure for the material optimization problem based on the maximum plastic dissipation principle.In Sect.5, we consolidate all our developments in an algorithmic framework and provide pertinent implementation details.Finally, we verify our framework with benchmark problems in Sect.6.

Thermomechanical formulation of elastoplasticity
We start by briefly reviewing the basic principles of elastoplasticity from a thermodynamics viewpoint, including the mechanical work identity, the notion of plastic dissipation from the second law of thermodynamics, and the derivation of the constitutive equations for associative plasticity reflecting on the principle of maximum plastic dissipation.We primarily follow the exposition of Simo and Hughes (2006).

The mechanical work identity
We consider a macroscale initial boundary value problem defined on a domain Ω and restrict our attention to a time interval [0, T] .The position of a material point in the domain Ω is denoted by x .The macroscale density at a material point x is denoted as (x) .The domain Ω is subjected to a traction t(t) at the Neumann boundary Γ N and the prescribed displace- ments ūE (t) at the Dirichlet boundary Γ D with a body force b(x, t) , where t ∈ [0, T] .Then, the macroscale displacement field ū(x, t) at a material point x and at time t ∈ [0, T] is a mapping ū ∶ Ω × [0, T] → ℝ 3 .We define the corresponding velocity and strain fields at (x, t) ∈ Ω × [0, T] as where sym(◻) represents the symmetric part of a second- order tensor.
With the kinematically admissible velocity field v(x, t) and the macroscale stress field (x, t) , the mechanical work iden- tity is where This identity directly comes from the application of the principle of virtual power (PVP) with a specific choice of the velocity field as the test function (Simo and Hughes 2006).The PVP is a fundamental exposition and can be applied for various applications ranging from, for instance, modeling DNA macromolecules to elastic foundations (Kalliauer et al. 2020;Höller et al. 2019).Please refer to Germain for systemic application of PVP to derive fundamental equations in continuum mechanics (Germain 1973).

Remark 1
We emphasize that we chose ( ◻) notation for the macroscale displacement ū , and the given boundary condi- tions ūE and t .Later in this article, the displacement solution and boundary conditions drive the optimization algorithm, and, therefore, this choice directly relates with the introduced notations for the sought solutions for the optimized multiscale configurations.

Constitutive relations from the second law of thermodynamics
The constitutive relations between the macroscale stress (x, t) and the macroscale displacement field ū(x, t) (through (2) (3) 195 Page 4 of 31 the macroscale strains E(x, t) ) close the global governing equations stated in (2) and (3).The second law of thermodynamics governs the form of these constitutive relations, and we summarize the important results in the following.First, we decompose the macroscale strain tensor E into an elastic and plastic part assuming small strains, denoted by E e and E p , as We introduce the notion of internal potential energy and dissipation within the context of elastoplasticity.We define the internal energy of the system as where Ψ(E e ) is the Helmholtz free energy density defined in terms of the stored elastic energy function W and the contributions from hardening effects.In this presentation, we consider the case of perfect plasticity, which implies that the contribution from hardening is zero and Ψ = W.
Next, we look at the difference between the stress power P int ( , v) and the rate of change of the internal energy V int , which we denote by D mech .Assuming isothermal conditions, the Clausius-Duhem version of the second law of thermodynamics (Truesdell and Noll 2004;Tadmor et al. 2012) follows as We identify D mech as the total instantaneous mechanical dis- sipation in the domain Ω at time t ∈ [0, T] , which is always non-negative.
Inserting the definitions of P int ( , v) and V int from (3) and (5) and using the strain decomposition (4), we arrive at where ( ◻) denotes the material time derivative of a quantity.The rate of elastic and plastic part of the macroscale strain tensor lie in the space of second-order symmetric tensors , that is Ėe , Ėp ∈ .Any permissible value of Ėe , Ėp ∈ will lead to a total strain rate, thanks to the additive decomposition Ė = Ėe + Ėp ∈ , which can describe a plausible kinematic process.The principle of thermodynamic determinism requires that (7) remains valid for any kinematic process, which implies The first equation is a typical elastic constitutive relation, where the stress is defined as the derivative of the free (4) ) E e and ∶ Ėp ≥ 0.
energy function with respect to the elastic part of the strain tensor.The second equation ( 8) represents the irreversible nature of an elastoplastic process implying that the dissipation energy is always non-negative.This relation constrains the possible stress states a material can undergo and indicates that the stress depends on the rate of the plastic part of the strain tensor.

The principle of maximum plastic dissipation
The principle of maximum plastic dissipation is a cornerstone of the mathematical formulation of associative plasticity.In the following, we derive the material constitutive equations for perfect plasticity from the viewpoint of this principle.We first assume a yield criterion ( ) , with ∈ denoting any possible stress state.Its zero isosurface is the usually convex yield surface that encloses the space of admissible stresses For a given plastic strain E p ∈ , we define the plastic dis- sipation D p at a material point for perfect plasticity as where ∈ Σ now denotes an admissible stress state.
In the local form, the principle of maximum plastic dissipation states that, for a given plastic strain E p ∈ , the plas- tic dissipation D p attains its maximum for the actual stress tensor among all possible stresses ∈ Σ .Mathematically, the principle is The classical formulation of associative plasticity (flow rule, loading/unloading conditions) directly follows from this principle.To this end, we first transform the maximization principle into a minimization problem by changing the sign of the objective function.Next, we transform the constraint minimization problem into an unconstrained problem by introducing the cone of Lagrange multipliers where L 2 denotes the space of all square integrable functions.The corresponding Lagrangian functional The solution to (11) is then given by a point ( , ) ∈ × satisfying the Karush-Kuhn-Tucker optimality conditions for the Lagrangian functional (13).The conditions are The first equation in ( 14) is the associated flow rule, often also called the normality of the flow rule.The second and third equations in ( 14) are the classical loading/unloading conditions.The only requirement for these equations to hold uniquely is the convexity of the elastic range Σ .A sufficient condition for this requirement is the convexity of the yield criterion function ( ) .We will exploit these aspects later on for devising solution strategies for our optimization framework.

Multiscaling concepts in continuum micromechanics
Continuum micromechanics forms a rigorous foundation for the (semi-)analytical estimation of homogenized stiffness and strength properties of materials with hierarchical microstructures.Here, we state the key results that are relevant in the context of this article.For a detailed review, interested readers are referred to the presentations given in Zaoui (2002), Suquet (1997).

Estimation of homogenized elastic properties
The goal of continuum micromechanics is to estimate the homogenized response of a representative volume element (RVE) filled with microheterogeneous material.For the existence of such an RVE, a minimal requirement is that the characteristic length, d, of the considered inhomogeneities and deformation mechanisms is much smaller than the size, l, of the RVE.Moreover, l must be much smaller than the characteristic length scale of the variation in the loading on the macroscale structure, L. Therefore, a proper scale separation implies In each phase r of the RVE, the average microscopic stress r , the average microscopic strain r , and the phase stiffness r are linked as: r = r ∶ r .The kinematic compatibil- ity for the homogeneous strain boundary conditions for the RVE relates the macroscale strain tensor E with the volume average of microscopic strains r .Similarly, the equilibrated microscopic stresses r and the macroscale stress tensor fulfill the volume average relation following the homogeneous stress boundary conditions.With r as the volume fraction of the phase r, these relations are ( 14) ≥ 0, ( ) ≤ 0, and  ( ) = 0.
A link between the macroscale strain E and the average microscopic strain r of phase r is established with a fourth order concentration tensor r as A comparison of the macroscale constitutive relation = ℂ ∶ E with ( 16) and ( 17) yields the homogenized stiff- ness ℂ in terms of the volume fraction, stiffness, and con- centration tensor of constituent phases as It is clear from (18) that the estimation of the concentration tensors entails the homogenized stiffness ℂ .The simplest choice for is to assume a uniform strain state throughout the RVE, that is = , where is a fourth-order sym- metric unit tensor.This choice leads to the well-established Voigt mixture rule for homogenized stiffness.However, the Voigt rule does not consider any other statistical information beyond the volume fraction of phases.We note that the Voigt rule is often applied for "homogenization/interpolation" between a solid material and voids in conjunction with relaxation to ill-defined 0-1 type problems in topology optimization (Bendsøe and Sigmund 1999;Allaire and Aubry 1999).
The estimation of the concentration tensor r based on Eshelby's matrix-inclusion solutions can incorporate the volume fraction, the shape of phases, and their interaction with each other.Eshelby's matrix-inclusion problem relates strains in an ellipsoidal inclusion perfectly bonded with the surrounded homogeneous infinite elastic matrix to the applied homogeneous strains at infinity.Following Zaoui (2002), the estimation of r from the matrix-inclusion problem entails the homogenized stiffness expression as Here, the Hill tensor ℙ 0 r characterizes the morphology of the inclusion phase r and its interaction with the surrounding reference matrix with stiffness tensor ℂ 0 .The Hill tensor ℙ 0 r depends on the morphology, that is, the shape and orientation of the inclusion phase as well as the stiffness tensor of (17) r = r ∶ E. .
195 Page 6 of 31 the reference matrix.The analytical expressions for ℙ 0 r can be found in (Laws 1977(Laws , 1985;;Masson 2008).With (19), the homogenized stiffness of the RVE can be expressed as a function of constituent phase characteristics.

Estimation of homogenized elastic limit strength
A macroscale RVE reaches the elastic limit state when any one of the constituents in the RVE yields.Let us focus on the weakest constituent phase, denoted by index r = w .We assume that its elastic limit behavior is described by the yield criterion where * w is the effective stress measure in the weak phase w.Moreover, we assume that it is the only constituent that exhibits inelastic behavior.The effective stress or "stress peaks" in phase w can then be estimated with the secondorder moment of the stress field in this phase, which is the quadratic stress average over the phase volume V w Suquet (1997), expressed as Next we assume that w is a scalar deviatoric stress-based yield criterion such as the von Mises criterion with known yield strength Y w , bulk modulus w and effective volume function φw for the weak phase w.Following (Suquet 1997), with the admissible stress and homogenized stiffness ℂ , the weak phase criterion w translates to the macroscopic yield criterion as We emphasize that ( 22) is of the for m of = √ ∶ ∶ − R ≤ 0 that represents the general quad- ratic form of classical rate-independent plasticity models.In the context of this article, it implies that the elastic domain defined by ( 22) satisfies two critical geometric properties.These properties are (1) the convexity of the elastic domain and (2) the degree-one homogeneity of the yield criterion.These properties are very important for developing solution algorithms for our optimization formulation, and we will recall them later in subsequent sections.

A framework for path-dependent concurrent material and structure optimization
In this section, we formulate a concurrent material and structure optimization method maximizing the path-dependent stiffness for elastoplastic multiphase hierarchical systems.
We then focus on the finite element discretization of this formulation.

Thermodynamically consistent formulation
We start by looking at a representative problem illustrated in (22)

Global optimization problem with micromechanical design variables
We introduce the definition of macroscale density (x) and microstructure characterization m(x, t) .We assume that the macroscale density (x) is fixed with respect to time, while m(x, t) is a function of loading history representing a local adaption of microstructure with time.The set m(x, t) contains the geometric and mechanical characterization of phases that span multiple well-separated microscales, consisting of volume fraction, material properties, shape, and orientation of the different phases in the hierarchical system.The homogenized material constitutive relations defined by the plastic dissipation D p and the Helmholtz free energy Ψ in Sect.2.1 depend on the macroscale density (x) and the microstructure characterization field m(x, t) .Therefore, the design vector is [ (x), m(x, t)] T .A typical objective is to maximize the structural stiffness for the path-dependent nonlinear structure designs.It translates as the maximization of the total mechanical work expended in the course of a deformation process (Fritzen et al. 2016).Assuming a quasi-static case with no inertial effects, the total mechanical work f w in the considered time interval [0, T] follows directly from the mechanical work identity (2) and the definition of stress power (6) as Utilizing the definitions of D mech , V int and P ext (v) from Sect.2.1, we arrive at We note that the (pseudo-)time t represents the loading history.
Augmenting D p and Ψ with (x) and m(x, t) , we set up our optimization problem using the introduced definition of total mechanical work f w from (24) as ( 23) (24) A ad and E ad define the set of admissible design variables at the macro-and microscales, respectively, with possible design constraints.The second part of this equation is the continuous version of the objective function proposed in Fritzen et al. (2016).On the one hand, the velocity field v(x, t) and, thus, the displacement field ū(x, t) implicitly depend on (x) and m(x, t) .On the other hand, the first part of ( 25) is an explicit expression in terms of the design variables (x) and m(x, t) with known macroscale strains E and E p that arise from the solution of the global equilibrium equations.The first part of ( 25) is the basis for the development and mathematical analysis of our optimization formulation.
Later on in this article, we will come back to the second part of ( 25) to motivate sensitivity calculations.

Definition of the sets of admissible design variables
The admissible set A ad seeks a limit on the total material mass M req available for design.Mathematically, it can be defined as where min and max are the bounds on the macroscale material density .Without loss of generality, the definition of the admissible set E ad is illustrated via the representative multiphase hierarchical system shown in Fig. 1.We observe a wellseparated three-scale hierarchical system with three base constituent materials denoted as Material A, B, and C with densities A , B , and C , respectively.Material A and B are linear elastic, and Material C exhibits a perfectly elastoplastic response.At a material point P, the volume fraction of Material B and C at the lowermost scale are B and C such that B + C = 1 .Material B forms spherical inclusions in the matrix of Material C at this scale.The homogenized material from this scale forms the matrix M that hosts the inclusions of Material A with the orientation A and elongation A at the mesoscale.The density of the matrix M is simply M = ( B B + C C ) .The volume fraction of Material A and matrix M are A and M with We note that this set can be arbitrarily extended if necessary.
With these definitions, the microscale design admissible set E ad for the representative multiphase hierarchical system shown in Fig. 1 follows as ( 26) at their respected scales.Furthermore, the elongation ratio of the inclusions of Material A is bounded by max .We note that the bounds are constant and do not depend on the loading history or the material position in the domain.We emphasize again that the multiscale configuration of Fig. 1 is used for illustration, but the underlying representation is easily generalized to cover any other multiphase hierarchical system.

Decomposition into material and structure optimization problems
The definition of the admissible set E ad at a material point x only depends on the macroscale density (x) of this point.It implies that the definition of E ad is pointwise, and we can thus rewrite the statement (25) as We also assume that the macroscale density (x) is fixed with respect to the loading history.Therefore, we are allowed to swap the integral and maximization operations in (28), hence The statement (29) allows us to decompose the optimization problem into two sub-problems.The outer "structure" optimization problem is ( 27) We also explicitly write the equivalent statement in terms of the total external work done in the deformation process from ( 25), which we will exploit in the subsequent sections for discretization and sensitivity calculations.In this statement, m(x, t) optimizes the following sub-problem or "material" optimization problem: The macroscale density (x) dictates the construction of the admissible space E ad , and, therefore, we take it out from the definitions of D p and Ψ in ( 31) and consider it in E ad .In the second line, we rewrite the definition of the plastic dissipation D p with the help of the principle of maximum plastic dissipation discussed in Sect.2.1.3.The constitutive relations in (31) defined through the Helmholtz free energy Ψ and the plastic dissipation D p remain to be discussed.For linearized elasticity, the stored elastic energy function W takes a quadratic form in the elastic part of the strain tensor E e = E − E p .The homogenized elasticity tensor ℂ is a function of the microstructure char- acterization field set m(x, t) ∈ E ad ( (x)) .For the perfect plasticity case, that is Ψ = W , the quadratic form follows as The elastic constitutive equation from ( 8) entails the following stress-strain relationship: Similarly, the elastoplastic material constitutive equations through the maximum plastic dissipation principle stated in ( 11) is augmented to include (x) and m(x, t) as where the definition of the elastic closure Σ in terms of (x) and m(x, t) is The homogenized stiffness ℂ(m(x, t)) and the homogenized yield criterion ( , m(x, t)) can be estimated as a function of microstructure variables for instance via the continuum micromechanics principles outlined in Sect.2.2.
Remark 2 Please note that, in this article, we assume perfect plasticity for the definitions of the Helmholtz free energy Ψ and the plastic dissipation D p appearing in the decomposed optimization formulation stated in ( 30) and ( 31).For modeling plasticity with hardening within the decomposed formulation, the Helmholtz free energy Ψ in (32) and the plastic dissipation D p in ( 35) can be augmented to include harden- ing contributions by utilizing continuum micromechanicsbased homogenization schemes, such as outlined in Fritsch et al. (2009); Morin et al. (2017) for the example of bones.Throughout this article, however, we keep the assumption of perfectly associated plasticity to focus on the foundations of the decomposed formulation and the general ideas for its computational treatment.

Interpretation as an inelastic constitutive law
The combination of ( 30) and ( 31) constitutes the concurrent material and structure optimization formulation.The maximization problem in (30) seeks the optimal material distribution (x) in the domain Ω .For a given material distri- bution (x) , the optimization problem (31) finds the optimal microstructure configuration maximizing the stress/deformation power for the known macroscale strains at each material point x .Both statements are coupled through the macroscale strains and, therefore, through the displacement field solution ū(x, t) that satisfies the global equilibrium equations.This interdependency makes the global equilibrium a constitutively nonlinear problem analogous to a typical initial boundary value problem with an inelastic constitutive law.Therefore, we propose to interpret the material optimization problem as a reformulated elastoplastic constitutive law that provides the locally optimal material response with respect to the external loading history.The microstructure variable m(x, t) can be thought of as an "internal state variable" anal- ogous to any path-dependent history variable encountered in elastoplasticity formulations.This interpretation will be used in Sect. 4 for devising the optimization algorithm for the material optimization problem.

Finite element discretization
In the next step, we discretize our concurrent material and structure optimization formulation within the context of the finite element method.We use vector-matrix notation, ( 35) consistent with the standard finite element discretization of the initial boundary value problem introduced in Simo and Hughes (2006), to represent the introduced quantities in the global equilibrium equations.

Model definitions in the discrete setting
We start by dividing the time interval [0, T] into n load parti- tions and split the domain Ω into N e finite elements: Here, Ω j is the domain of element j, and each element is equipped with N gp Gauss quadrature points.We focus on a typical (quasi-)time interval [t n , t n+1 ] with known equili- brated state at time step t n .In the context of this work, we choose standard nodal finite elements with Lagrange basis functions that can be assembled to approximate the macroscale displacements and strains at load increment (n + 1) over the complete domain: where N is the (assembled) displacement interpolation oper- ator and B is the (assembled) strain-displacement opera- tor (Hughes 2000).In the sense of the standard Galerkin method, we use the same finite element basis functions for the representation of the solution and the test functions (Hughes 2000).

Remark 3
We emphasize that from here on, displacementtype vector quantities with subscript n or (n + 1) such as u n+1 denote the vector of unknown displacement-type coef- ficients at the corresponding load increment.Tensor quantities with subscript n or (n + 1) such as macroscale strains E n+1 or stresses n+1 denote the corresponding tensor fields approximated in terms of the corresponding displacement finite element solution at the corresponding load increment.
The design vector [ (x), m(x, t)] T for our example mul- tiscale configuration in Fig. 1 can now be defined in this discrete setting as [ , m] T , where where x ∈ {1, .., N gp }, j ∈ {1, .., N e }, n ∈ {0, .., n load − 1}.The macroscale density j is assumed to be constant in each element and load increment, with j being the element index.The microstructure design variable set m is defined at each (macroscale) Gauss point and load increment with m n+1 as the microstructure characterization set at load increment (n + 1) .The microstructure configuration m x,j n+1 at a Gauss point x inside element j at load increment (n + 1) consists of volume fraction The macroscale plastic strain E p n is known from the equilibrated solution state at load step n.We note that we write these constitutive relations in tensor notation given its direct relation with continuum micromechanics principles stated in Sect.2.2.

Discrete form of the structure optimization problem
With the introduced definitions, we can write the discrete form of the material and structure optimization formulation ( 30) and (31).In this work, we employ the trapezoidal rule for the numerical evaluation of the integrals over the (quasi-) time domain, which is second-order accurate with respect to the (quasi-)time step size.The discrete version of the structure optimization problem (30) then becomes In accordance with the notation introduced above in the context of a finite element discretization, f ext n+1 is the external force vector, ūn+1 is the converged vector of the macroscale nodal displacements, and Δ ūn+1 ∶= ūn+1 − ūn is the incre- ment of the displacement vector in load increment (n + 1) .The force residual rn+1 is calculated utilizing the optimal microstructure configuration mn+1 in load increment (n + 1) .We note that the optimal microstructure configuration mn+1 and ūn+1 are dependent on each other justifying the choice of ( ◻) notations introduced in Sect.2.1.M( ) is the total mass of the occupying domain, and j and |Ω j | are the density and volume of element j.
The total mechanical work f w is the discrete version of the second statement in (30), which is equivalent to the objective function proposed in Fritzen et al. (2016).This essentially is the area under the characteristic force-displacement curve approximated with the trapezoidal rule, as illustrated in Fig. 2. The first condition in (42) ensures that the global equilibrium is satisfied in all load steps.The second and third conditions of (42) are the discrete definitions of the macroscale admissible design variable set A ad .The total (42) available mass M req can be expressed in terms of the frac- tion M frac with respect to the mass when the densest material occupies the complete domain.
The force residual rn+1 at load increment (n + 1) is defined as where w x contains the Gauss point weight and the deter- minant of the Jacobian matrix for element j.We observe that the microstructure design variable set m is implicitly accounted for by the residual definitions in each load increment.The macroscale stress n+1 at each Gauss point is evaluated by solving the nonlinear elastoplastic constitutive relations ( 39) and ( 41) with known microstructure configuration mx,j n+1 that solves the material optimization problem detailed in the following subsection.Therefore, the global equilibrium equation ( 43) is nonlinear and requires iterative solution approaches such as the Newton-Raphson incremental procedure (Simo and Hughes 2006;de Souza Neto et al. 2011).
Structure optimization involving inelastic material models is potentially ill-posed in a force-controlled setting (Swan and Kosaka 1997; Huang and Xie 2008;Schwarz et al. 2001;Maute et al. 1998;Cho and Jung 2003).Therefore, we only consider displacement-controlled loading in this article, incrementally applied through prescribed displacements ūE (t) , see Sect.2.2.1.In a displacement-controlled setting, f ext n+1 represents the discretized form of the loading potential resulting from the non-zero displacement boundary conditions.This assumption also simplifies the sensitivity calculations of the objective function with respect to the design variables for the optimization algorithms, which we will discuss in Sect.5.1.

Discrete form of the material optimization problem
For a given material distribution and the macroscale displacement solution vector ūn+1 , the material optimization prob- lem for the Gauss point x inside element j for load increment (n + 1) follows from (31) as The first part of this equation directly comes from the incremental form of the principle of maximum plastic dissipation outlined in (41).Similarly, the second part is the incremental form of the Helmholtz free energy rate defined in (31).( 43) We emphasize that all quantities at load increment n are known, and, therefore, Ψ(E n − E p n ) does not play any role in this maximization problem.The optimized configuration mx,j n+1 is sought in the microscale design variable space m ] with constraints defini- tions that follow from the admissible set E ad (27).From (44), we can rewrite the material optimization statement as including all constraints defined through the stress admissible set Σ n+1 and microscale design admissible set E ad .The first two conditions are essentially the elastoplastic constitutive equations relating the macroscale stress with the macroscale strains via (39) and the constraint on the macroscale stress defined through the homogenized yield criterion (41).The third condition is the definition of the Helmholtz free energy in terms of microscale design variable m x,j n+1 . The rest of the conditions follow in a straightforward manner from the constraints definitions in E ad .The solution of (45) at each Gauss point in each load increment yields the optimized microstructure configuration set m.
We emphasize that in contrast to the equivalent strain energy maximization for concurrent optimization problems involving overall linear elastic multiphase hierarchical systems (Xia and Breitkopf 2014;Gangwar and Schillinger 2021), finding the solution to the material optimization problem (45) is not straightforward.Both the macroscale plastic strain E p n+1 and the optimized microstructure mx,j n+1 are unknown.Intuitively, the material optimization problem maximizes the area under the homogenized elastoplastic stress-strain curve for each material point.Multiple stress-strain curves are available at each load increment, defined by the the microscale design variable m x,j n+1 . This inter-dependency couples the history variable E p n+1 with m x,j n+1 , necessitating a challenging novel algorithmic treatment to tackle this maximization problem. (45) 195 Page 12 of 31

Algorithmic treatment of the material optimization problem
For the algorithmic treatment, we interpret the material optimization problem as a reformulated constitutive law at each material point that provides a locally optimal mechanical response to the loading history.This interpretation allows us to treat the microscale design variable m x,j n+1 as an additional internal state variable within the context of the classical formulation of plasticity.With this interpretation, we first exploit the principle of maximum plastic dissipation to motivate a solution strategy for the material optimization problem.We then cast this strategy into an algorithmic procedure that assumes the format of a typical return map algorithm for the integration of elastoplastic constitutive equations.Finally, we leverage continuum micromechanics and the associated homogenized elastoplastic constitutive relations, which enable further simplifications that make our framework computationally feasible.

The principle of maximum plastic dissipation revisited
As explained above, the first part of the material optimization problem ( 44) is the incremental statement of the maximum plastic dissipation principle.This part defines the interaction between the next stress state n+1 and the optimized microstructure state mx,j n+1 through the homogenized yield criterion ( , m x,j n+1 ) .Focusing on this part only, we can combine both statements in a single one as where This maximization problem seeks the macroscale stress state n+1 and a solution mx,j n+1 within the modified admissible space definition Σ n+1 .The solution mx,j n+1 restricts the search space for the solution mx,j n+1 of the original material optimization problem (45), which we will further detail in the subsequent discussion.We note that the interpretation of the microscale design variable m x,j n+1 as internal state variable naturally arises from these statements.
We define a Lagrangian functional that converts the constraint optimization problem (46a) into an unconstrained problem following (13): where is in the cone of Lagrange multipliers defined through (12).The solution to (46a) is given by a point ( n+1 , mx,j n+1 , Δ n+1 ) that satisfies the Karush-Kuhn-Tucker optimality conditions for (47).The conditions entail The general structure of ( 48) is similar to the typical local constitutive equations for plasticity (flow rule, loading/ unloading conditions) as described in ( 14).Equation ( 48) 2 represents the evolution of microstructure state in a particular load increment (n + 1) .All these equations together form a coupled nonlinear system that requires a special computational treatment.The solution mx,j n+1 from (48) provides important insights into the interaction of the plastic update and the microstructure update.In the following, we will can utilize these insights to design an algorithmic framework for solving the original material optimization problem (44).

Algorithmic procedure in the form of return map algorithms
Analogous to the elastic-plastic operator split formulas for inelastic constitutive equations, we define a trial elastic state by freezing the plastic flow and microstructure evolution state during the current load increment.It implies that the macroscale plastic strain and optimal microstructure configuration state in the current load increment are known and equal to that of the previous load increment, that is With these assumptions, the trial elastic state is Page 13 of 31 195 Here, E e,tr n+1 , tr n+1 and tr n+1 denote the macroscale trial elastic strain, trial elastic stress, and trial yield criterion, respectively.Further, the convexity of the homogenized yield criterion leads to the following important property Simo and Hughes (2006) where n+1 is the homogenized yield criterion calculated at the macroscale stress n+1 and the optimal microstructure configuration mx,j n+1 after load increment (n + 1).

Elastic trial state and plastic vs. microstructure updates
The trial state (49) and property (50) in combination with (48) lead to three important cases defining the restrictions posed by the solution mx,j n+1 on the search space for the original material optimization problem (45).In the following, we will discuss all three cases in detail.
Case 1 If tr n+1 < 0 , then property (50) entails n+1 < 0 , implying a purely elastic step.In this case, the solution ( 49) .Therefore, with zero dissipation, the solution mx,j n+1 of the material optimization problem reduces to the strain energy maximization that follows from (45) as In conclusion, the solution mx,j n+1 does not pose any restrictions to the search space E ad for the optimized microstruc- ture configuration mx,j n+1 in this case.To support our discussion of the remaining two cases, Fig. 3 presents a geometric interpretation in a typical returnmapping context.We note that in Fig. 3, the gray region represents the family of available homogenized yield criterion envelops ( , m x,j n+1 ) = 0 at each load increment, defined by the set of microscale design variables m x,j n+1 . Case 2 If tr n+1 > 0 , and if it is possible to find at least one microstructure configuration mx,j n+1 ∈ E ad ( j ) such that (51) n+1 ∶= ( tr n+1 , mx,j n+1 ) ≤ 0 , property (50) indicates that tr(2) n+1 ≥ n+1 ⟹ n+1 < 0 .The solution is possible via the microstructure update provided that new microstructure state mx,j n+1 follows the constraint ( tr n+1 , mx,j n+1 ) ≤ tr(2) n+1 ≤ 0. We call this case adaption to elastic state through microstructure evolution.We observe in Fig. 3a that due to the current microstructure state mx,j n (dashed line), the material goes into plastic state; however, the material adapts itself to fall back to the elastic state by updating the microstructure state to mx,j n+1 (dotted line), while maximizing the total strain energy.
Again, equation ( 48) leads to E p n+1 = E p n and Δ n+1 = 0 following the discussion in Case 1.The solution mx,j n+1 follows as The constraint in this problem ensures that the state remains elastic and can be interpreted as a restriction on the search space for mx,j n+1 posed by the feasible solutions of (48).We write in the strain space to emphasize that the strain state is known and problem ( 52  to Δ n+1 > 0 from (48) 1 , and the microstructure evolution condition from (48) 2 entails ∕ m x,j n+1 = 0 .It means that the microstructure configuration remains unchanged, that is mx,j n+1 = mx,j n .It implies that the solution mx,j n+1 of (48) restricts the search space for mx,j n+1 to a single point, that is mx,j n .With mx,j n+1 = mx,j n , the rest of the relations in (48) reduces to the typical elastoplastic constitutive equations with known stiffness and yield criterion.This is illustrated in Fig. 3b, where we see that no elastic update is possible for the current trial state; therefore, the microstructure state remains unchanged, and the next stress state n+1 is solved with a standard return map algorithm such as the closest point projection algorithm (Simo and Taylor 1985).

An analogue to the elastoplastic return map algorithm
For an intuitive understanding, Fig. 4 presents a graphical solution of the material optimization problem for the example of a one-dimensional linear elastic-perfectly plastic model.The gray region in these graphs represents the family of stress-strain curves for different microstructure design configurations m x,j n+1 . The material state (stress, strain, microscale configuration) at load increment n is known, and the macroscale strain E n+1 at load increment (n + 1) is given.Typical strain increments are infinitesi- mal, and increments are large for illustration only.The next material state from the material optimization problem warrants that the area increment (red shaded region in the graphs) is maximized.
The trial elastic stress state tr n+1 assumes the microstructure state and the macroscale plastic strain in this load increment are equal to the previous load increment.Case 1 results in a purely elastic update as shown in Fig. 4a.Except for the first load increment, this leads to a trivial solution for linear-elastic perfectly plastic models with the same microstructure configuration, that is mx,j n+1 = mx,j n .Case 2 in Fig. 4b is of particular interest.The trial stress tr n+1 predicts a plastic update.However, it is possible to find material configurations mx,j n+1 such that ( tr n+1 , mx,j n+1 ) ≤ 0 .The material adapts itself by falling back onto the elastic state via an appropriate update of the microscale configuration mx,j n+1 (denoted with the red dashed line), maximizing the total strain energy following (52).In Case 3, no material configuration allows an elastic state for the trial stress tr n+1 .Therefore, the material configuration remains unchanged, and the stress-strain state is updated through a return-mapping/closest point projection algorithm as shown in Fig. 4c.
We cast these cases into an algorithmic frame analogous to a standard elastoplastic return map algorithm.We summarize the result in the following box.

A special choice: continuum micromechanics based homogenization
In this article, we focus on homogenized yield criterion based on the quadratic stress average that can be obtained within a continuum micromechanics framework as reviewed in Sect.2.2.2.To illustrate the simplifications that can be achieved with this choice, we consider again the one-dimensional example from Fig. 4. Using the homogenized yield criterion based on continuum micromechanics, we arrive at Fig. 5 that graphically illustrates the associated simplifications in the algorithmic procedure.
In particular, the microscale configuration corresponding to the maximum stiffness (maximum strain energy density) also results in the maximum strength properties for 195 Page 16 of 31 the homogenized response, which implies that adaption to elastic state through microstructure evolution (Case 2) cannot occur here.In the following, we will provide a more detailed account of these simplifications.

Special properties induced through this choice
For illustration purposes, we fall back to our initial example of a representative multiphase hierarchical system defined in Fig. 1.We recall that Material C at the microscale is a perfectly elastoplastic material that follows the von Mises failure criterion with yield strength Y C and bulk modulus C .For this example, we can then write the homogenized yield criterion given in ( 22) as a function of the stress and microscale design variable where φC is the equivalent volume fraction of Material C computed as φC = (1 − . For a detailed derivation of the homogenized stiffness ℂ(m x,j n+1 ) , we refer inter- ested readers to Appendix 1 in Gangwar and Schillinger (2021).These estimates hold the following three properties that form the basis for further simplifications in the algorithmic procedure for the material optimization: Property 1 The microscale configuration mx,j n+1 corresponding to the maximum stiffness (maximum strain energy density) from (51) also maximizes the homogenized strength response.Figure 5 graphically represents this property for the one-dimensional case with a family of possible stress-strain curves.Here, the microscale configuration corresponding to the higher linear-elastic slope leads to a higher limit strength for the homogenized response.Thus, the stress-strain curve for the configuration mx,j n+1 acts as an envelope for all the stress-strain curves defined by the possible microscale configurations m x,j n+1 . Utilizing the definitions , this property can be summarized as ( 53) x,j n+1 ) ≤ 0 is not possible for any microstructure configuration m x,j n+1 . Therefore, following our discussion in Sect.4.2, we can conclude that an adaption to the elastic state through microstructure evolution (Case 2) is inconceivable for the continuum micromechanics schemes outlined in this paper.
Property 2 An important conclusion from the previous section is that the microscale design update is possible in an elastic step only.The elastic part of macroscale strain E e n+1 ∶= E n+1 − E p n+1 at each Gauss point therefore entails the optimal material orientation θx,j A,n+1 for load increment (n + 1) .In the elastic step, the material optimization problem is essentially a strain energy maximization.The maximum strain energy is obtained for a general orthotropic material by aligning the material axis with the principal strain axes for the elastic strains (Jog et al. 1994;Pedersen 1989).
Property 3 If the external loading is monotonically increasing, the optimal material orientation θx,j A,n+1 is the only microscale variable that may change in each load increment.We denote the set of remaining microscale design variables as ] .The optimal configuration ml(x,j) n+1 for m l(x,j) n+1 remains unchanged throughout the loading history, that is ml(x,j) n+1 = ml(x,j) n ∀n = 1, 2, ..., n load − 1 .We pro- vide a proof of this property in A.

Simplification of the material optimization problem
The three special properties discussed above entail two important simplifications.First, adaption to the elastic state through microstructure evolution (Case 2) cannot occur.Second, except for the material orientation , the optimized material configuration remains unchanged throughout the loading history.This implies that the material optimization problem is solved for the first load increment only via the strain energy maximization (51) for the optimized configuration mx,j n+1 .Later, the optimized material orientation is updated for each load increment by aligning the material axis with the principal strain axes of the elastic part of the macroscale strain tensor E e n+1 .We would like to emphasize the crucial role of the continuum micromechanics based estimates for the concurrent material and structure optimization of multiphase hierarchical systems.These estimates render both objective functions and constraint definitions of the material optimization statement (45) and, therefore, the strain energy maximization (51) as "discretization-free" analytical or semi-analytical expressions.This reduced problem is a straightforward constraint optimization problem that can be solved with standard gradient-based methods.The solution to this material optimization problem is equivalent to solving a set of (n + p) nonlinear equations with (n + p) variables, where n and p are the total number of microscale design variables and the total number of equality constraints, respectively.Thus, the continuum micromechanics enables us to handle computational challenge as well as the complex constraint definitions in the optimization of multiphase hierarchical systems.For a detailed discussion on these aspects, please refer to our previous work (Gangwar and Schillinger 2021).
In conclusion, the total cost of solving all the material optimization problems is equivalent to the case of an endcompliance type optimization problem with a linear elastic response at the material scales.These simplifications result in an enormous reduction in computational effort, making our framework computationally tractable for the elastoplastic case.

Comments on computer implementation
In this section, we provide an overview of our optimization framework with a focus on essential computer implementation details.First, we derive the essential sensitivity calculations of the objective function f w with respect to the design variables for the structure optimization problem (42).We then briefly touch upon the optimality criteria method for updating the design variables in each structure optimization iteration, utilizing the computed sensitivities (Sigmund 2001).Finally, we consolidate all developments into a single algorithmic framework.

Sensitivity analysis
The format of the structure optimization problem ( 42) is equivalent to the topology optimization formulation for elastoplastic structures presented by Fritzen et al. (2016), Xia et al. (2017).They derive the sensitivities using the path-dependent adjoint method (Buhl et al. 2000;Cho and Jung 2003).In the following, we provide a sketch of the derivation and highlight the important results.For a detailed derivation, we refer interested readers to Fritzen et al. (2016), Xia et al. (2017).
The adjoint method begins with the construction of a Lagrangian function f * w that satisfies the zero residual constraints rn+1 and rn at (quasi-)time t n+1 and t n for each term of the trapezoidal rule stated in ( 42).With the Lagrange multipliers n+1 and n+1 that are of the same dimensions as the vector of unknowns ūn+1 , the Lagrangian function f * w follows as Since rn+1 and rn vanish at the equilibrium solution, the sen- sitivity of f * w is same as that of f w , implying that The derivative of f * w with respect to the design variable j follows from (55) as The derivative of rn+1 with respect to j is evaluated follow- ing the residual definition in (43).Substituting the definition of n+1 given in (39) into (43), the derivative expression becomes (55) 195 Page 18 of 31 K tan n+1 is the global finite element stiffness matrix of the mechanical system at the equilibrium of load step (n + 1) .We note that the quantities inside the square brackets in the ( 58) and second term of this equation are only computed for element j.For all remaining elements, n+1 does not depend on j .The second term is zero for all corresponding entries, maintaining dimensional consistency with the vector f ext n+1 .We observe that the sensitivities of f w as expressed in ( 57) and ( 58) require computationally extensive calculations of unknown derivatives.Therefore, our aim is to obtain the values of the Lagrange multipliers n+1 and n+1 in such a way that these unknown derivatives can be eliminated from the sensitivity expression.To this end, we classify the degrees of freedom (DOF) into essential (index E; associated with the Dirichlet boundary conditions) and free (index F; remaining).According to this classification, we can partition vectors and matrices as shown for the following generic objects v and M: Since the displacements ūE on the Dirichlet boundary Γ D are prescribed, they are independent of the current value of the optimization variable .This observation leads to at an arbitrary load step index q = 0, ..., n load − 1 .With dis- placement-controlled loading, the only possible non-zero entries in the global force vector f ext q are the reaction forces The relations ( 60) and ( 61) lead to an educated choices for the vectors n+1 and n+1 such that the unknown derivatives with respect to the design variables in ( 57) and ( 58) can be eliminated (see Fritzen et al. (2016); Xia et al. (2017) for details).The final expression for the sensitivity of the objective function f w with respect to the design variable j is With the prescribed displacement increments Δ ūE n+1 , the choice of Lagrange multipliers n+1 and n+1 that lead to the above expression is Set up load increments ∆ ūE n+1 for each loading step index n = 0, .., n load − 1; Compute the macroscale strain Set up the linear system: , and solve for δ ū(k+1) ; 18 Apply Newton correction to the displacements: We note that the history of kinematic state variables E and E p in (62) are known from the solution of the global equilib- rium equations at each load increment.For the representative multiscale configuration in Fig. 1, the derivative of the homogenized stiffness ℂ with respect to the element density j in (62) can be evaluated by means of the chain rule.Following Gangwar and Schillinger (2021), the expression is where x,j A,n+1 and x,j C,n+1 relate to j via (45).The partial derivatives of ℂ with respect to x,j A,n+1 and x,j C,n+1 are evaluated at the optimal microstructure configuration mx,j n+1 by using finite difference approximations.

Structure optimization scheme
We utilize the algorithmic procedure for the structure optimization that we outlined in our previous work Gangwar and Schillinger (2021) and that we briefly summarize in the following.First, we define sensitivity numbers to rank the element sensitivities that are used to update the macroscale design variables in each design iteration: To avoid mesh dependency and checkerboard patterns, the sensitivity numbers are first smoothed with a filtering scheme defined as where N j is the set of neighboring elements for which center- to-center distance Δ(j, j � ) to element j ′ is smaller than the filter radius r min .To improve convergence, the sensitivity numbers are further averaged with the sensitivity numbers of the previous design iteration as The ratio of sensitivity numbers and the mass constraint are combined to ( 63) (64) where Λ i is the Lagrange multiplier corresponding to the total material mass constraint in design update i, and is a damping parameter.The macroscale density is updated by means of the well-known optimality criteria method Sigmund ( 2001): To prevent a singular global stiffness matrix, the lower limit min on j is limited by a small value, set in our case to 0.001.The maximum possible element density max depends on the density of the constituents at the microscales and the prescribed bounds in (45). is a small move parameter that improves the stability, for instance by preventing multiple holes appearing and disappearing during optimization.The Lagrange multiplier Λ i is updated using the bisection method to satisfy the mass constraint.The design iterations stop when the density convergence criteria are met.

General algorithm
Algorithm 1 consolidates all the developments into an algorithmic framework.It mainly consists of three blocks.
The outer block represents macroscale structure optimization iterations with iteration index i, using the optimality criteria method detailed in (69).It stops when the relative change in macroscale density falls below the tolerance tol , and the converged solution is the optimum macroscale density ρ .For a given macroscale density distribution, the middle block solves the initial boundary value problem with known load increments Δ ūE n+1 at each load increment n.The global equilibrium for each load increment is solved with the Newton-Raphson method that uses the linearization of (43) (Simo and Hughes 2006).Here, (•) (k)  n+1 denotes the value of a particular variable (•) at the kth iteration at load step (n + 1) .The Newton-Raphson scheme stops when the norm of the residual force vector drops below a tolerance threshold tol , and we adopt tol = 10 −5 in this article.
The inner block solves the material optimization problem described in Sect. 4 at each Gauss point with prescribed state variables at this iteration stage for each load increment.For the schemes based on continuum micromechanics outlined in this paper, the material optimization problem is solved for the optimized configuration mx,j n+1 maximizing the strain energy via (51) at the first load increment only.Thereafter, the material orientation is updated for each load (68) increment by aligning the material axis with the principal strain axes of the elastic part of the macroscale strain tensor according to our discussion in Sect.4.3.The equations due to strain energy maximization can be solved with standard gradient-based constraint optimization methods such as the quasi-Newton method of Broyden, Fletcher, GoldFarB, and Shanno (BFGS) or sequential least squares programming (SLSQP) methods.

Numerical examples
In this section, we define two test examples with elastoplastic multiphase hierarchical material models that are suitable to illustrate the computational efficiency and validity of our path dependent concurrent material and structure optimization framework.First, we consider a standard cantilever type benchmark problem and modify its material definitions analogous to the multiscale configuration shown in Fig. 1.Later, we demonstrate the potential of our framework for biotailoring applications by solving a prototype problem that integrates a hierarchical material model for cereal stems Gangwar et al. (2021).

Problem description and hierarchical design
Figure 6 modifies the definition of the standard cantilever design problem to demonstrate the developed concepts in this article.The length and height of the macrostructure are 2.0 m and 1.0 m, respectively.The left edge is fixed, and a displacement loading of u * = 7.5 mm is prescribed at the central 10% of the right edge that we divide in six load steps with a constant load increment of Δ ūE = 1.25 mm.We discretize the macroscale structure with an 80 × 40 mesh of 4-node plane strain quadrilateral elements, resulting in a characteristic element size of l e = 25 mm and 3, 200 mac- roscale design variables.Each element contains four Gauss points, resulting in 80 × 40 × 4 = 12, 800 material optimiza- tion problems in each load step.
As illustrated in Fig. 6, we consider a hierarchical system that consists of Material A, B, and C at two different length scales.Their densities (in Kg/m 3 ) are A = 0 , B = 0.5 , and C = 1.0 , their Young's moduli (in GPa) are E A = 0.0 , E B = 0.5 , and E C = 1.0 , and Poisson's ratio of all constituents is 0.3.Material C is elastoplastic with yield strength 1 MPa.We assume that Material A forms cylindrical inclusions in the homogenized matrix of Material B and C. At each Gauss point, the material microstructure is parametrized by the volume fraction , the orientation , and the volume fraction C > h , where h is a small positive number.We restrict min to 0.001 and max to 0.799 to satisfy these requirements.The total amount of material mass available is restricted to 40% of the maximum possible mass.As an initial condition at the macroscale, we assume the maximum possible density max in each element.At the material level, we assume an initial microstructure configuration A = 0.0 , A = 0.0 , and C = 1.0 at each Gauss point.In each design update, we reduce the target mass fraction by 0.025 until we reach the specified mass fraction M frac = 0.4 .The move parameter and the damping parameter are set to 0.05 and 0.5.The filter radius r min is reduced linearly from r min = 20 l e to r min = 4 l e with design iterations for improv- ing the convergence of the structure optimization algorithm following Xia et al. (2017).
The structure optimization algorithm stops when the relative change in the macroscale density field falls below the tolerance tol = 10 −3 .Figure 7 illustrate the convergence of the macroscale design update.We notice that the algorithm takes 34 density updates to converge to the final design with converged objective function value 0.39022 N-mm. Figure 8a and b plot the optimized macroscale density and equivalent plastic strains overlaid on the density plot, respectively.The plastic strains are concentrated at the clamped end's boundaries, consequently pushing the material towards these regions.As highlighted in Fig. 8b, the sharp features near the clamped end in the optimal density distribution mimic the plastic front emphasizing its importance for the final design.
Figure 9 illustrates the optimized morphology at the mesoscale and the equivalent volume fraction of Material B and C from the lowermost scale.The yellow color in Fig. 9a represents the matrix material that results from the homogenization of the lowermost scale, and the blue color displays the volume fraction and orientation of Material A inclusions.The inclusions follow the direction of the largest principal stress.The equivalent volume fractions of Material A, B, and C at the macroscale are defined as: φA =  A , φB = (1 −  A )(1 −  C ) , and φC = (1 −  A ) C . Figure 9b displays the equivalent volume fraction of Material B and C at the macroscale for the final design, where we use 60% opacity for both.We can observe the regions dominated by Material B, C, and a mixing zone.The stiffer Material C is deposited in the regions anticipated to yield first, while Material B dominates the transition zone.
Figure 10 illustrates the evolution of the optimization process by plotting equivalent plastic strains overlaid on the density distribution and the equivalent volume fraction of Material B and C, all at selected design iterations.The evolution of the macroscale density and equivalent plastic strains shows that the design process attempts to attenuate the plastic front.In this process, the algorithm pushes more material towards the region close to the clamped boundary, delaying yielding in this region.Material C is the stiffest material among the constituents and exhibits elastoplastic behavior.Its evolution is heavily influenced by the plastic front, which leads for instance to sharp features in the macroscale density configuration.

Performance comparison with equivalent linear elastic design
We finally illustrate the impact of elastoplastic design on the structural performance by comparing with a corresponding design that only assumes linear elastic material response at the microscales and ignores any plasticity effects.To this end, Fig. 11 plots the optimized density distribution and the equivalent volume fraction of Material B and C for the linear elastic design that assumes purely elastic properties of Material C at the lowermost scale.Comparing these plots with Fig. 8 and 9b, we can find apparent differences in the optimized layouts.The plastic design places more material towards the clamped region with clear features imitating the plastic front, whereas these attributes are missing in the linear elastic design.
Figure 12 quantitatively compares the structural performance of the elastoplastic design over the corresponding linear elastic design.We subject the optimal configurations in both cases to the same displacement loading of u * = 7.5 mm, whereas Material C is elastoplastic for both configurations (see Fig. 6). Figure 12 demonstrates the load-displacement curves with the equivalent plastic strains plots for different load levels for both cases.Although the response at low load levels is practically identical, the load-displacement curves start to deviate from each other at the higher load levels, when elastoplastic behavior originating from material scales governs the overall response.The features in elastoplastic design highlighted in Fig. 8b plays a crucial role in attenuating the propagation of plastic front by delaying the plasticization of the hierarchical material system.Figure 12 demonstrates the differences in yielded regions and the magnitude of equivalent plastic strains in both cases, which clearly shows the superior behavior of the elastoplastic design.Moreover, the relative percentage gain in the structural performance for the elastoplastic design, defined as ((f w,nl − f w,lin )∕f w,lin ) × 100 , is 4% , where f w,nl and f w,lin are objective function values for the elastoplastic and equivalent linear elastic designs, respectively.This gain is expected to grow with the applied load level u * .Thus, we conclude that the linear elastic design and the plastic design are functioning differently, and it is important to consider plastic effects at different scales in multiphase hierarchical systems that are expected to develop dissipation-based energy absorption mechanisms against external impacts.

Towards predicting self-adapting mechanisms in plants
Biomaterials exhibit multiscale inelastic behavior and develop dissipation-based energy absorption mechanisms optimizing its hierarchical composition across material microscales against the external biophysical stimuli through natural evolution (Wegst et al. 2015;Fratzl and Weinkamer 2007;Bhushan 2009).A rational understanding of microstructure interdependencies with self-adapting mechanisms will pave the way towards many biotailoring applications with improved properties, for instance, in the context of the targeted breeding of agricultural crops (Brulé et al. 2016;Berry et al. 2004).A few studies have attempted the multiscale optimization of biological systems such as boneremodeling and bioinspired materials (Rodrigues et al. 1999;Coelho et al. 2008;Radman et al. 2013).Several roadblocks, however, such as high computational cost and non-trivial problem decomposition in the case of elastoplastic behavior have limited these approaches to simple linear elastic problems with no more than two scales.With the following prototype model, we demonstrate the potential of our optimization framework in overcoming these roadblocks for the computationally efficient modeling of self-adaption of biomaterials.
Crop stem materials organize themselves hierarchically across multiple length scales.The hierarchical scales in crops range from base constituents such as cellulose, hemicellulose, and lignin, to cell wall, functional tissues, crosssection, and structure scale node morphology levels.In our previous work, we experimentally profiled this hierarchical organization through microimaging (micro-CT, light microscopy, transmission electron microscopy) and chemical analysis, focusing on cereal stems (Gangwar et al. 2021), which we briefly summarize in Fig. 13.Exploiting this data, we developed and validated a continuum micromechanics model of cereal stem materials that accurately relates material composition with elastoplastic mechanical behavior across different scales.We provide all implementation information relevant in the scope of this article in 1.
Figure 14 summarizes the prototype model for the hierarchical optimization of a cereal node region, given that stem failure has been generally observed in this region.We take the length and height of the macrostructure domain as 10 mm and 5 mm, which is typical of a node dimension.The left edge is completely fixed, and the displacement in X direction on the right edge is fixed.The prescribed displacement in the Y direction on the right edge is u * = 0.4 mm.The displacement u * = 0.4 mm is divided in four load steps with increment Δ ūE = 0.1 mm.We discretize the macroscale structure with a 60 × 30 mesh of 4-node plane- strain quadrilateral elements.This macrostructure model definition is equivalent to a cereal node cross-section that undergoes combined shear and bending loads.
Following our multiscale material model, the stem cross-section consists of an outer-shell layer and a solidpith region.The primary functions of the outer-shell are non-mechanical, such as protecting against insects and regulating gas exchange.Thus, we only consider the solidpith region for hierarchical optimization.The microstructure design variables m x,j n+1 consist of the cell wall fraction par(x,j) wall,n+1 in the parenchyma, the fiber fraction x,j fib,n+1 in the vascular bundles, the vascular bundle fraction x,j vb,n+1 , and the orientation x,j n+1 of the anisotropy axis of the solid pith with respect to the global X direction (see Fig. 14).Lignin in the parenchyma cell wall material exhibits elastoplastic behavior.The parenchyma tissues and xylem-phloem vessels in the vascular bundles are also responsible for food storage and nutrient-water transport.We incorporate these biological constraints by adopting the bounds on the volume fractions that we measured through microimaging  Gangwar et al. (2021) in the material optimization problem.At the structure scale, the total amount of material is restricted by the reported average density, which can be interpreted as the limitation posed by the available biological energy required in the synthesis of biomass per unit of stem material.
Figures 15 and 16 illustrate the design evolution history and the final macroscale density with the equivalent plastic strains.The macroscale design algorithm takes 38 density updates to converge to the final design.In the optimal layout, the branches from the left internode converge to the central node region, and the branches of the right internode emerges, a morphology that was also observed in real plants through micro-CT images (Ghaffar and Fan 2015).The plastic strains are concentrated at the end and middle regions due to the anticipated high shear deformations.The optimal density layout puts material in areas to attenuate plastic fronts, which can also be observed in the design evolution history in Fig. 15.These observations again emphasize the role of the plastic front in the optimal structural layout.
Figure 17 plots the optimal microstructure configuration at different scales at the final load level.At the mesoscale, the material orientation follows the stress flow direction in the main branches, while the morphology is more complex in the central node region.We also plot the optimal configurations at lower material scales in the main branches.These results indicate the choice of a stronger solid-pith material for the optimal mechanical response.The predicted morphology is in qualitative agreement with what our collaborators in plant science have found via field experiments for lodging-resistant cereals (Gangwar et al. 2023).Based on these promising results, we believe that our optimization framework can help pave the way towards efficient and sustained biotailoring applications, supported by modeling and simulation.

Summary, conclusions, and outlook
In this article, we established rigorous theoretical foundations for an efficient concurrent material and structure optimization framework for multiphase hierarchical systems with elastoplastic constituents at the material scales.In particular, we developed an efficient solution strategy for the material optimization problem based on the maximum plastic dissipation principle in the format of a typical return map algorithm for an elastoplastic constitutive law.Finally, we integrated analytical expressions of the homogenized stiffness and the yield criterion that are derived via continuum micromechanics, enabling a computationally tractable implementation for elastoplastic multiphase hierarchical systems.
We verified the validity and efficiency of our framework with newly defined benchmark problems that, for the first time, is computationally feasible via our framework.It consists of a macroscale configuration in the form of a standard cantilever, but involves hierarchical material definitions with elastoplastic constituents at the microscale.The optimized macroscale and microstructure configurations computed via our framework demonstrated the importance of plasticity effects that originate from the material microscales in developing dissipation-based energy absorbing mechanisms.In addition, we applied our framework for investigating selfadapting mechanisms in cereal plant structures, outlining its potential for biotailoring applications.Our framework is a first attempt at a decomposed material and structure formulation that optimizes the path-dependent macroscale mechanical response of elastoplastic multiphase hierarchical systems.We would also like to point out the limitations of the presented framework.In this article, we restricted ourselves to the class of inelastic hierarchical materials, for which we can assume the existence of an associative flow rule, a rate-independent ideal elastoplastic response, and an isothermal process at the macroscale.We also confined the quadratic stress average based micromechanical scheme for estimating the homogenized yield criterion by an additional assumption that only one of the constituents in the hierarchical material exhibits inelastic behavior with a deviatoric stress-based yield criterion.These assumptions restrict the plausible elastoplastic failure mechanisms that originate from the material microscales.We anticipate that these mechanisms can be incorporated into our optimization framework, for instance by combining the variational procedure for incremental homogenization and transformation field analysis (TFA) (Dvorak and Benveniste 1992;Brassart et al. 2011).Moreover, the formulation can potentially be extended for other path-dependent problems, where nonlinear effects such as viscoplasticity, fracture, and damage originate from the material microscales.Thus, our framework helps push forward path-dependent concurrent material and structure optimization to consider nonlinearities exhibited at the microscales, with a number of potential applications, including multiscale additive manufacturing and architecting metamaterials (Meza et al. 2015;Sanders et al. 2021).

Appendix A Comments on continuum micromechanics-based simplifications
We present a proof of Property 3 introduced in Sect.4.3.First, we write the strain energy maximization expressions in index notation.We denote the orthonormal basis that corresponds to the global coordinate system as {e p } .We drop superscript (x, j) from variables for conciseness.Problem (51) at load increment n with E e n ∶= E n − E p n can be rewritten in index notations as The principle coordinate system, that is co-linear with the principal strain directions, uses the Roman indexed basis {ê i } .The transformation matrix Q pi between both systems depends on the optimal configuration θA,n from Property 2.  Utilizing Q pi and its orthogonal property, we define the ten- sor component transformations as where Êe ij(n) are the components of the elastic part of the macroscale strain tensor in the principle coordinate system.Also, Êe ij(n) = 0 if i ≠ j , and diagonal components ( i = j ) are the principle strain values.Similarly, Ĉijkl are the compo- nents of the stiffness tensor in the principal coordinate system.Using (71), we reformulate the maximization problem (70) in terms of m l n = [ A,n , A,n , C,n ] as With ml n as the solution to this maximization problem, we can rewrite (72) as We assume that the macroscale loading increases monotonically.Therefore, the elastic part of the macroscale tensor Êe ij(n+1) load increment (n + 1) in the principal coordinate system can be written in terms of the components Êe ij(n) with an appropriate scaling.Exploiting the definition of the Kronecker delta , we write Êe ij(n+1) in terms of scaling components a i as As the scaling components are non-negative, we augment the expression (73) and arrive at This is the expression for the strain energy maximization problem at load increment (n + 1) analogous to (72) or (73).
We emphasize that our definition of the admissible set E ad depends only on the macroscale density distribution and, therefore, remains unchanged throughout the loading history.In addition, Case 2 is not conceivable as laid out in our discussion of Property 1 in Sect.4.3.Therefore, it implies that the possible admissible solutions of the right hand side expression at increment (n + 1) in 75 are the same as those at increment n.

Appendix B Details on the cereal prototype model
We briefly summarize the key components for implementing the hierarchical optimization of the cereal prototype model in Sect.6.2.For further details on multiscale modeling of stiffness and strength of crop stem material within continuum micromechanics, we refer the interested reader to Section 3 in Gangwar et al. (2021).In the current model, the microscale design variables defined at each Gauss point at load step (n + 1) are the cell wall fraction par(x,j) wall,n+1 in the parenchyma, the fiber fraction x,j fib,n+1 in the vascular bundles, the vascular bundle fraction Here, ℂ par and ℂ vb are the homogenized stiffness tensors of the parenchyma tissue and the vascular bundle tissue in the solid-pith region, respectively.For the analytical expression of these estimates, interested readers are referred to Section 3.3 in Gangwar et al. (2021).The composition of all other RVEs in the multiscale material model of cereal stems is considered constant, and model parameters corresponding to the Gopher oat variety are used (see Appendix B in Gangwar et al. (2021)).T is a standard rotation matrix for tensor transformations.
Lignin exhibits elastoplastic material behavior at the constituent level, and the macroscale limit state point corresponds to the yielding of lignin.In this prototype model, we assume that only lignin in parenchyma cell wall material is elastoplastic (see Fig. 14).In this case, the macroscale homogenized yield criterion reads as in the parenchyma cell wall material is fixed as given in Gangwar et al. (2021).
With these definitions in hand, we write the material optimization problem for a Gauss point x inside element j for load increment (n + 1) following ( 45) as where j is the given macroscale dry density for the finite element with index j.The first three lines in the constraints definition represent the microstructure dependent constitutive equations.The fourth statement connects j with the microscale design variables via the rule of mixture.We adopt bounds for par(x,j) wall,n+1 , x,j fib,n+1 and x,j vb,n+1 that are [0.01,0.38] , [0.75, 0.90] , and [0.01, 0.16] , respectively.The upper bound reflects the measured microscale parameters reported in Gangwar et al. (2021).The strain energy maximization part in the algorithmic treatment of the material optimization problem is a constraint optimization problem with nonlinear equality constraint.We utilize the sequential least squares programming (SLSQP) method implemented in the SciPy library to solve this problem.The setup of the structure optimization problem for this prototype model resembles (42).The sensitivity analysis for the macroscale design updates follows from Sect.5.1.The derivative of the homogenized stiffness ℂ with respect to the element density j follows via the chain rule as The derivatives of ℂ with respect to the microscale design variables at the material level are evaluated by finite difference approximations.The move parameter and the damping parameter are set to 0.02 and 0.5.The filter radius r min is reduced linearly from r min = 15 l e to r min = 4 l e with design iterations.x,j fib,n+1 x,j fib,n+1 j + ℂ x,j vb,n+1 x,j vb,n+1 j .

Fig. 1 .
We assume a fixed reference domain Ω subjected to a traction t(t) at the Neumann boundary Γ N and the prescribed displacements ūE (t) at the Dirichlet boundary Γ D with a body force b(x, t) , where t ∈ [0, T] .

Fig. 1
Fig. 1 Sketch of a representative elastoplastic multiphase hierarchical system Page 8 of 31 Here, the volume fraction of Material A is bounded by min A and max A , and the volume fraction of Material C is bounded by min C and max C

195
Page 10 of 31 for Material A, and volume fraction x,j C,n+1 of Material C. We again emphasize that we use this definition of the multiscale configuration in Fig. 1 for illustration purposes, this procedure can easily be generalized to cover any other multiphase hierarchical system.The constitutive equation relating the macroscale stress n+1 with the macroscale strains E n+1 and E p n+1 at the Gauss point x follows from (33) as where the homogenized stiffness ℂ(m x,j n+1 ) is evaluated for the microstructure configuration m x,j n+1 .To derive the incremental form of the elastoplastic constitutive equations, the discrete version of the maximum plastic dissipation principle at the Gauss point x from (35) is where the admissible stresses lie in a set Σ n+1 defined by the homogenized yield criterion ( , m x,j n+1 ) evaluated at m x,j n+1 ∶

Fig. 2
Fig. 2 Total mechanical work f w in the course of the deformation process

Fig. 3
Fig. 3 Geometric illustration of solution strategy for the material optimization problem, based on the return-map algorithm interpretation n+1 ) = 0 does not have any solution.It implies that no microstructure state can solve (48) with the chosen trial elastic strain E e,tr n+1 .Therefore, only a plastic update is feasible, and E

Fig. 4
Fig. 4 Graphical solution of the material optimization problem in one dimension for the three possible cases

Fig. 5
Fig. 5 Simplifications in the algorithmic procedure for material optimization induced by continuum micromechanics estimates.The optimized material configuration mx,j n+1 remains unchanged with loading history

56
Initialize load increment couter n = 0; Initialize ū0 = 0 =⇒ E 0 = 0, and E p 0 for load step (n + 1) , which results in 38, 400 microscale design variables in each load step.The minimum volume fraction of Material A is set to min A = 0.2 .The existence of the homogenized yield criterion in (53) requires φC = (1 −  x,j A,n+1 )  x,j C,n+1 > 0 .It implies that the bounds  max A < 1 − h and  min

Fig. 7
Fig. 7 Convergence of objective function f w (in N-mm) and mass fraction with respect to number of macroscale design iterations

Fig. 8
Fig.8Macroscale density distribution and equivalent plastic strain distribution of the cantilever benchmark problem for a total prescribed displacement of u * = 7.5 mm.Highlighting circles demonstrates that the plastic front influences the optimal density distribution

Fig. 10
Fig. 10 Evolution of macroscale density configuration and equivalent plastic strains (rainbow colormap, ×10 −3 units ) and equivalent volume fractions of Material B and C

Fig. 11
Fig. 11 Final design of the cantilever benchmark in the equivalent linear elastic case (Material C purely elastic)

Fig. 15
Fig. 15 Convergence of the objective function f w (in 10 −3 N-mm) with respect to the number of design iterations, plotted along with snapshots of the macroscale density configuration and the equivalent plastic strains (in ×10 −3 units, rainbow colormap) ) e p ⊗ e q ∶ C pqrs (m n )e p ⊗ e q ⊗ e r ⊗ e s ∶ E e rs(n) e r ⊗ e s ) C pqrs (m n ) E e rs(n) .

Fig
Fig. 17 Optimal microstructure configuration for total prescribed displacement u * = 0.4 mm of the anisotropy axis of the solid-pith material.The macroscale homogenized stiffness tensor ℂ in the global coordinate system can be written as a function of the microscale design variables m inGangwar et al. (2021)  as

26 end Algorithm 1
k+1) ; Concurrent structure and material optimization framework for elastoplastic structures with multiphase hierarchical materials.
25Update density ρ i+1 using the optimality criteria algorithm; i++; Q pi e p ; Êe ij(n) = Q pi Q qj E e pq(n) ; C pqrs = Q pi Q qj Q rk Q sl Ĉijkl ,