Concurrent Topology Optimization of Multi-Scale Composite Structures Subjected to Dynamic Loads in the Time Domain

: This paper presents a concurrent topology optimization of multi-scale composite structures subjected to general time-dependent loads for minimizing dynamic compliance. A three-ﬁeld density-based method is adopted to implement the concurrent topological design, with macroscopic effective properties of the microstructure evaluated through energy-based homogenization method (EBHM). Transient response is obtained from the two-scale ﬁnite element analysis with the HHT-α approach as an implicit time integration procedure. Design sensitivities are formulated employing the adjoint variable method (AVM) based on two main philosophies: “discretize-then-differentiate” and “differentiate-then-discretize” approaches, respectively. The method of moving asymptotes is adopted to update the design variables at two scales. Several benchmark examples are presented to demonstrate that the “discretize-then-differentiate” AVM attains consistent sensitivities in an inherent manner such that the resulting optimal topology is more efﬁcient when compared with the “differentiate-then-discretize” AVM. Moreover, the potential of the proposed method for concurrent dynamic topology optimization problems under general time-dependent loads is also highlighted.


Introduction
Additive manufacturing process enables the fabrication of structures in the light of an expected macrostructure layout along with underlying microstructures. This offers significant design space for designers to create lighter and more efficient structures. Concurrent topology optimization provides a rigorous mathematical framework for seeking optimized material distribution at macro and micro scales to achieve superior structural performances. Therefore, they are of great interest for exploring multi-scale modeling and design methodology in this exciting field [1][2][3].
The two-scale concurrent topology optimization framework simultaneously optimizes two sets of design variables representing respective layout of the macrostructure and periodic unit cell. This framework is widely applied to two-scale hierarchical structural design issues, such as static compliance [4][5][6], eigenfrequency [7][8][9], structural modal damping ratio [10], as well as thermomechanical behavior [11,12]. Bai et al. [4] introduced a two-step Helmholtz filtering/projection scheme to describe the shell interface, whereby a multi-scale topology optimization model for shell-infill structure is developed for minimizing the static compliance. Gangwar et al. [6] presented a concurrent material and a structure design framework considering shape and orientation of various phrases in a hierarchical system across multiple various length scales. Xiao et al. [7] designed graded lattice sandwich structures in terms of maximal natural frequency through multi-scale topology optimization, which is employed to integrate the optimization of thickness of two solid face-sheets and layout of lattice cells into a core layer. Zhang et al. [8] extended the work of Xiao et al. [7] to inhomogeneous cellular structures for maximizing the eigenfrequencies of desired modes based on mode-tracking strategy. Hu et al. [9] performed the multi-scale topology optimization of coated structures with multiple layers of graded lattice infill for maximization of the fundamental eigenfrequency. Ni et al. [10] proposed an optimization strategy to maximize the structural damping performance, where the damping material layout and its microstructural configuration are concurrently optimized. Ali et al. [11] formulated the concurrent multi-scale and multiphysics topology optimization for minimization of the thermal and mechanical compliances. Zhou et al. [12] designed lightweight channel-cooling cellular structures with eminent heat barrier and load-carrying capacity via metamodel-assisted concurrent multi-scale and multi-material topology optimization. For a comprehensive review on concurrent multi-scale topology optimization, one can refer to the published literature [13].
Despite this, certain challenges still remain in some efficient cumbersome sensitivity analysis and dynamic response analysis across multiple scalers for hierarchical structures under dynamic load. Concurrent topology optimization for dynamic response was investigated in both the frequency domain [14][15][16][17][18][19][20] and the time domain [21,22]. This work concentrates on a transient response optimization problem for minimizing the dynamic compliance of multi-scale composite structures under general time-dependent load. Millions of design variables for transient problems of multi-scale structures pose great significance to efficient sensitivity analysis when gradient-based topology optimization algorithm is implemented. Therefore, the adjoint variable method (AVM) is essential for sensitivity analysis. There are two dominant philosophies to implement the AVM in terms of the order of discretization and differentiation regarding the time variable, i.e., differentiatethen-discretize method and discretize-then-differentiate approach. Zhao et al. [22] adopted the AVM based on a differentiate-then-discretize approach to conduct the sensitivity analysis for transient concurrent topology optimization of two-scale hierarchical structures. Majority of investigations adopted the differentiate-then-discretize approach for linear transient problems due to its relative simplicity in formulation and implementation [22][23][24][25][26]. Nevertheless, Jensen et al. [27], Zhang et al. [28] and Ding et al. [29] demonstrated that the differentiate-then-discretize AVM can cause consistency errors representing differences between the calculated and accurate sensitivities through investigating a single DOF damping system. Alternatively, AVM based on a discretize-then-differentiate approach can diminish resulting consistency errors associated with the differentiate-then-discretize approach. Giraldo-Londono et al. [30] proposed a transient topology optimization implementation of an elastodynamic system employing the discretize-then-differentiate AVM, whereafter their work was further extended to local stress-constrained topology optimization problem with arbitrary dynamic loads [31]. Other studies, such as microstructural layout optimization of viscoelastical component under time-dependent loading and transient thermomechanical coupling problems have also been based on the differentiate-then-discretize AVM [32,33]. Recently, Kristiansen [34] developed a completely parallel framework to address the largescale transient topology optimization employing the fully discretized adjoint sensitivity analysis in [35]. Nevertheless, to the author's knowledge, very few investigations on multiscale concurrent topology optimization adopting the differentiate-then-discretize AVM are focused on linear transient problems due to comparatively cumbersome sensitivity analysis.
This work intends to construct an efficient two-scale concurrent topology optimization framework for minimizing the dynamic compliance of composite structures under transient loading. A three-field density-based method is exploited for multi-scale concurrent topology optimization to achieve material-structure integrated designs. The major contributions of this study consists of three aspects: (1) to formulate an efficient sensitivity computation for transient response optimization of two-scale hierarchical structures; (2) to demonstrate and discuss some findings in concurrent topology optimization aiming at the dynamic compliance minimization in the context of linear transient problems; and (3) to indicate the capabilities of the proposed concurrent topology optimization approach to design composite structures suffering from general transient loads.
The remainder of this paper is organized as follows. Section 2 briefly reviews the problem formulation of concurrent topology optimization for minimizing dynamic compliance of two-scale composite structures in the time domain. We present the HHT-α method in Section 2, followed by the adjoint sensitivity analysis via the discretize-then-differentiate approach in Section 3. Next, the inconsistent sensitivity via the differentiate-then-discretize approach is formulated in Section 4. Section 6 explains that the order of differentiation and discretization plays a critical role in the consistency of adjoint sensitivity analysis, and demonstrates the potential of the proposed approach to address a wide variety of concurrent topology optimization problems under general transient loading, with four numerical examples. Finally, the conclusions of this work are presented in Section 7.

Concurrent Topology Optimization for Dynamic Compliance Minimization
The concurrent topology optimization framework is presented to simultaneously achieve the optimal macrostructure and material microstructure for minimal dynamic compliance in the time domain. Material microstructure is assumed to be uniform in the composition of a macrostructure for convenient manufacturing. This framework is briefly outlined to comprehend the fundamental procedure of performing concurrent topology optimization in this section.

Three-Field Density-Based Approach
We adopt the three-field density-based approach [36,37] to guarantee clear topologies in two scales. Two sets of design variables are separately defined, namely macroscopic design density in structural design domain and microscopic design density in a unit cell. Each design variable ranges from 0 to 1. To diminish the chessboard pattern and meshindependence, the original design variables are regulated with a smooth regularization filter [38] and expressed as follows: where Φ i is the neighboring set of elements within a specified filter radius R in the macroscopic design domain that have a center located at the centroid of the ith element and Ψ i is the neighboring set of elements within a specified filter radius r in the unit cell that have a center located at the centroid of the jth element. v mac k is the volume of element k in the macroscopic design domain and v mic l is the volume of element l in the unit cell. The weighting factors w mac ki and w mic lj are defined using a linearly decaying function: where x and y denote the center position of elements in both macro and micro design domains, respectively. To achieve the clear black-white design, Wang et al. [39] modified the linearly filtered design densities in Equations (1) and (2) employing a threshold projection function: where the physical design variables, ξ i and η j , use the Ersatz parameters much less than one denoted by ξ min and η min , respectively, to inhibit numerical instabilities of the stiffness and mass matrices when ξ → 0 and η → 0 . β mac and β mic are exploited to regulate the aggressiveness of the projection function. ξ th and η th are the threshold density specified as 0.5 in this work.

Numerical Homogenization
To attain the clear configuration at both scales, the material interpolation schemes with penalization are employed. At the microscale, the modulus matrix of an element within the cellular microstructure is interpolated via SIMP [40]. At the macro-scale, the modulus matrix of an element within the macrostructure with porous material is interpolated with RAMP [41].
where D B is the elastic constitutive matrix of base material and D H is the effective macroscopic constitutive matrix, which is computed as follows: where I denotes a unit matrix, b denotes the strain matrix at the microscale and u m denotes the unknown displacement field excited by the unit test strains in the microstructural domain Ω m . The resultant displacement matrix u m is obtained through resolving the following unit cell equilibrium problem with periodic boundary conditions: where the stiffness matrix k mic is given by the following: To prohibit the local eigenmodes occurring in regions with low densities, the polynomial function, as suggested by [42], is selected to penalize the macroscopic element stiffness matrix via the RAMP model: where the penalization exponent p is set to be 3. In addition, the effective mass density of corresponding periodic cellular material is represented as follows: where ρ B is the physical mass density of the base material. The energy-based homogenization method (EBHM) was employed to calculate the macroscopic effective properties of the porous material [43].

Formulation of Compliance Minimization
When a two-scale hierarchical structure is excited by a transient external load, the finite element equation used to solve the boundary value problem for this elastodynamic system is expressed as follows: where M, C, and K represent the global mass, damping, and stiffness matrices, respectively.
u t , and u t are the respective acceleration, velocity, and displacement vectors in response to the force vector f t at time step t. N is the number of analysis steps. The global mass and stiffness matrices are assembled using the penalized macroscopic element matrix: where N is the matrix of shape functions and B is the first derivative of N.
We employ the Rayleigh damping to compute the damping matrix as linear combination of mass and stiffness matrices, such that where α r and β r are the respective mass and stiffness proportional damping coefficients, which are assumed to be design-independent in this work. This work aims to minimize the dynamic compliance for a two-scale hierarchical structure with the limited available amount of material in the time domain. Mathematically, we formulate this two-scale concurrent topology optimization problem as follows: where f (ξ, η, u(t)) is the concerned objective function, V mac and V mic are the respective volumes of macroscopic and microscopic design domains. ς and ϑ are the volume fraction upper bounds associated with macroscopic and microscopic constraints of G 1 and G 2 , respectively. Here, the macroscopic design domain is discretized into N mac elements, while the unit cell is discretized into N mic elements.

HHT-α Method
We apply the HHT-α method, a well-developed implicit time integration scheme, to solve the second-order initial value problems stated as Equation (14). Due to an unconditional stability along with a second-order convergence [44,45], the HHT-α method have been used for linear and nonlinear structural dynamic analysis [46,47]. The HHT-α method is characteristic of superior numerical dispersion and energy dissipation by introducing a parameter α into the Newmark method to control the numerical damping. Accordingly, the motion Equation (14) representing the dynamic equilibrium is modified as follows: The HHT-α method adopts finite difference relationships from the Newmark-β method and hence the recursive formula of displacement and velocity is determined with the following: where the Newmark parameters β and γ are constants which control the integration accuracy and stability, respectively, by satisfying the following relationship: By substitution of Equations (22) and (23) into Equation (21), the time-discretized motion equation in residual form is derived as follows: where Following a standard HHT-α scheme, we can obtain the dynamic response at each time step. We resolve Equation (25) for .. u t and thereupon compute u t and . u t by applying the Newmark-β Formulas (22) and (23), respectively. As for .. u 0 , by assuming . u 0 and u 0 to be design-independent, it can be computed using the following residual equation:

Adjoint Sensitivity Analysis Using Discretize-Then-Differentiate
We apply the discretize-then-differentiate AVM to construct the corresponding adjoint equation on the discretized elastodynamic system in space and time. The standard AVM sensitivity analysis is performed following two essential procedure. First, some residual equations are added into the objective function to develop an augmented function. Then, this augmented function is differentiated and the adjoint variables are derived by vanishing the derivative terms of state variables regarding design variables.
In terms of the chain rule, the sensitivities of both the objective and constraint functions with respect to the original design variables can be calculated as follows: The sensitivity of f with respect to the arbitrary design variable x(ξ i ,η i ) is also written as follows: In order to facilitate the sensitivity analysis, we transform Equations (22) and (23) into the following residual form: Sequentially, we add adjoint variables λ t , µ t and ζ t and rewrite Equation (38) as follows: From Equations (39) and (40), it is obvious that ∂P t /∂x = 0 and ∂Q t /∂x = 0. Due to the design-independence of the initial conditions, ∂u 0 /∂x = 0 and ∂ . u 0 /∂x = 0. We employ these simplifications and eliminate all implicit terms including ∂u/∂x, ∂ . u/∂x and ∂ .. u/∂x in Equation (41), such that the following adjoint equations can be obtained: By substituting the residual Equations of (25), (29), (39) and (40) into the adjoint Equations of (42) and (43), we obtain the solution of the adjoint problem as follows: Using the adjoint solution from Equations (44)-(47), we rewrite Equation (38) as follows:

Sensitivity Analysis for Design Variables at the Macroscale
Provided that the concurrent optimization problem (20) applies macrostructural density relevant information via the stiffness interpolation function, E = [E i ] = g ξ i , and the volume interpolation function, V = [V i ] = ξ i , it facilitates recasting the sensitivity information of macroscopic design variables according to these fields. Therefore, we compute the sensitivity of f with respect to ξ by chain rule as follows: where the sensitivities of f regarding the macroscopic element volume fractions and stiffness parameters can be attained as demonstrated in Equation (48).
The terms, ∂R t /∂E i and ∂R t /∂V i , are evaluated in terms of Equations (25) and (29), respectively. There is a case for t = 0.
where subscript (i, t) denotes the field vector of element i at time step t and subscript (i, t − 1) denotes the field vector at time step t − 1. From Equation (12), the partial derivative of E i with respect to ξ i is computed as follows: As such, the sensitivity of the objective function regarding the macroscopic design variables can be obtained by substituting Equations (50)-(56) into Equation (49), where the adjoint variables are solved using Equations (44)-(47).

Sensitivity Analysis for Design Variables at the Microscale
Due to the effective material properties as a bridge between macro and microstructures, it is convenient to obtain the sensitivity information for the microscale design variables in the light of these homogenized parameters. For transient response problems, the sensitivity of f regarding the microscale design variables is recast via chain rule as follows: where The sensitivity of f with respect to the effective material properties can be attained from Equation (48), i.e., where ∂ f /∂D H = 0 and ∂ f /∂ρ H = 0, according to the objective function as shown in Equation (20). Similarly, the partial derivatives, ∂R t /∂D H and ∂R t /∂ρ H , are evaluated using Equations (25) and (29), and for i = 0: for t = 1, · · · , N:

Solution Procedure
The flowchart of the proposed concurrent topology optimization for multi-scale structures is depicted in Figure 1.

Solution Procedure
The flowchart of the proposed concurrent topology optimiza structures is depicted in Figure 1. This procedure launches through inputting the FEM information material properties and boundary conditions) and the optimization projection parameters, filter radius and penalty parameter), followed of design variables. Then, on the basis of the current design variable mass density and the constitutive matrix are obtained via EBHM. Th of the multi-scale structure is computed using the HHT-α method w function and constraints are directly calculated. Subsequently, the adj ysis is performed based on the discretize-then-differentiate approach. This procedure launches through inputting the FEM information (i.e., the mesh, base material properties and boundary conditions) and the optimization parameters (i.e., the projection parameters, filter radius and penalty parameter), followed by the initialization of design variables. Then, on the basis of the current design variables, the homogenized mass density and the constitutive matrix are obtained via EBHM. The transient response of the multi-scale structure is computed using the HHT-α method whereby the objective function and constraints are directly calculated. Subsequently, the adjoint sensitivity analysis is performed based on the discretize-then-differentiate approach. Finally, the Method of Moving Asymptotes (MMA) [48] is employed to update the design variables. This optimization process is terminated once a certain convergence criterion is met.

Adjoint Sensitivity Analysis Using Differentiate-Then-Discretize
The differentiate-then-discretize AVM constructs the adjoint equation in a semidiscretized dynamic system on the basis of spatial discrete and time continuous field variables, and subsequently the transient response is evaluated at each time step. We rewrite an objective function Φ in the following integral form: where J is the duration of the dynamic event and t is the continuous time variable.
We introduce the motion Equation (14) into Φ and thereby obtain the sensitivity Φ by standard AVM: where the prime denotes differentiation regarding the design variables and λ denotes the smooth adjoint variable. Through twice integrating-by-parts, we rearrange Φ as follows: where we employ the assumption that the external load, as well as the initial condition, is design-independent for simplification. To remove the response derivatives u ' (J) and To transform the adjoint problem into the initial value problem, we use a variable transformation t = τ(s) = J − s and then construct a composite function Λ satisfying Λ(s) = λ(τ(s)). Accordingly, we rewrite Equation (68) by transforming all the terms including u ' and . u ' :

Λ(τ(s)) + KΛ(τ(s)) ds
To annul all the terms containing u ' and . u ' , we formulate the adjoint variable Λ as follows: where the sensitivity is simplified as follows: where * denotes the convolution operator. Following the obtained displacement and velocity, u t and . u t , we approximate the original objective function employing the rectangular formula: Based on the discretized adjoint variables Λ n solution from Equation (71), the sensitivity of objective function is approximated as follows: In virtue of the order of differentiation and discretization, this method is featured as differentiate-then-discretize in that we first differentiate the augmented objective function to achieve Equation (72) and subsequently implement the time discretization to achieve Equation (74). This approach is seemingly elegant since the resultant adjoint transient problem is similar to the primal problem. Nevertheless, the method encounters the notably inconsistent sensitivity, as indicated in the following numerical examples. Since the resultant optimal configuration is based on the objective function sensitivity, gradient-based topology optimization demands the precise sensitivity information to design variables. We examine the efficiency of both discretize-then-differentiate and differentiate-then-discretize approaches for AVM sensitivity analysis by comparing them with the sensitivity evaluated through the finite difference method (FDM).

Numerical Examples
This section offers four benchmark cases to validate the proposed approach: a cantilever beam, a clamped beam, a support structure, a building and a simply supported 3D structure. We compare the two-scale optimal results obtained from Zhao et al. [22] based on the differentiate-then-discretize AVM with those from this manuscript, based on the discretize-then-differentiate AVM in the given four examples. For all numerical examples, we adopt the damping coefficient α = 0.05 and determine the Newmark constants β and γ by employing the formulas β = (1 + α) 2 /4 and γ = (1 + 2α)/2, for at least second-order accuracy and unconditional stability, respectively. Moreover, in every example, we first verify the validness of the discretize-then-differentiate method for AVM sensitivity analysis, and then investigate the influence of loading parameters on the optimum solution using the transient concurrent topology optimization based on the discretize-then-differentiate AVM. All the programs in four benchmark cases are written with the available version of MATLAB 2021.

Cantilever Beam Design under Half-Cycle Sinusoidal Load
As depicted in Figure 2, a cantilever beam is subjected to a concentrated half-sine load vertically exerted at the midpoint of right free edge. The geometrical dimension of the cantilever beam is as follows: length L = 8 m, height H = 4 m and thickness h = 0.01 m. For a composite structure with uniform microstructure, its Young's modulus is 200 GPa, Poisson's ratio is 0.3 and mass density is 7800 kg/m 3 . The Rayleigh damping parameters α r and β r are assumed to be 10 s −1 and 1 × 10 −5 s, respectively. The macroscopic and microscopic design domains are discretized into 5000 and 2500 four-node square quadrilateral elements, respectively. The maximal volume fraction for the macrostructure is specified to be 50%, and that for the unit cell is defined as 50%. To solve this problem, we adopt the input parameters listed in Table 1. α r and β r are assumed to be 10 s −1 and 1 × 10 −5 s, respectively. The macroscopic microscopic design domains are discretized into 5000 and 2500 four-node square qua lateral elements, respectively. The maximal volume fraction for the macrostructur specified to be 50%, and that for the unit cell is defined as 50%. To solve this problem adopt the input parameters listed in Table 1.    Table 2 compares the design sensitivity between two AVM approaches, discretize-thendifferentiate and differentiate-then-discretize, and FDM for this cantilever beam problem with a load application duration of t f = 0.05 s. We demonstrate the consistency error, namely the relative difference normalized by the exact sensitivity through FDM. It can be found that the sensitivities obtained with discretize-then-differentiate are significantly consistent with those obtained through FDM. However, the differentiate-then-discretize AVM induces significant inconsistent sensitivities. Figure 3 presents the iteration histories of the objective function and the constraint, and the optimized solution obtained using these two AVM-based sensitivity analysis techniques with a load application duration of 0.05 s. Obviously, the optimized configuration via discretize-then-differentiate is more favorable due to a lower value of the objective function. Thus, these results show that the order of differentiation and discretization has obvious effect on the consistency errors, which in turn can produce the inefficient optimum design. histories of the objective function and the constraint, and the optimized solution obtaine using these two AVM-based sensitivity analysis techniques with a load application dur tion of 0.05 s. Obviously, the optimized configuration via discretize-then-differentiate more favorable due to a lower value of the objective function. Thus, these results sho that the order of differentiation and discretization has obvious effect on the consisten errors, which in turn can produce the inefficient optimum design.

DH =
Microstructure and elastic matrix   Figure 4 shows the iterative histories during concurrent optimization and the resulting optimal designs for the load application duration of t f = 0.03, 0.01 s. It is seen that the optimized topologies at two scales for t f = 0.05 s ( Figure 3b) and t f = 0.03 s (Figure 4a) are nearly identical to each other. But in the case of t f = 0.01 s (Figure 4b), the optimal configurations greatly distinguish from the counterparts obtained for larger load application duration. When t f = 0.01s, substantial porous material is placed near the free edge of this beam, which produces inertial force to offset the short impulsive loading, which is is subsequently verified with the dynamic response as plotted in Figure 4. Also, the optimized macrostructure links the large mass distributed towards the free edge to the bracing ends by two horizontal beam-like members, which contribute to reducing the vertical bending of the beam. Figure 5 depicts the vertical displacement history at the load application point and the transient dynamic compliance of the beam. These results demonstrate that the transient responses for t f = 0.05 s resemble those for t f = 0.03 s, while they are significantly different for t f = 0.01 s. Particularly for t f = 0.01 s, the optimal beam continuously deflects downward, although the external load gradually decreases following t f = 0.005 s, which attributes to the fact that the resulting inertial force is sufficiently large to drive the beam downward.

Clamped Beam Design under Half-Cycle Cosine Load
In this example, we consider a beam fixed at both ends and excited via a concentrated half-cosine load vertically at the center of the bottom edge, as shown in Figure 6. The macroscopic design domain has length L = 12 m, height H = 2 m and thickness h = 0.01 m. We adopt a linear elastic material with a Young's modulus of 200 GPa, Poisson's ratio of 0.3 and mass density of 7800 kg/m 3 . The macroscopic design domain and the unit cell are discretized by respective 5000 and 2500 bilinear square elements. The volume fraction limits for both macrostructure and unit cell are set to be 0.5. The Rayleigh damping parameters α r and β r are the same as those in the first example. To solve this problem, we adopt the input parameters listed in Table 3. Table 3. Input parameters used to solve the clamped beam problem.

Parameter Value
Simulation   Figure 4 shows the iterative histories during concurrent optimization and the resu ing optimal designs for the load application duration of = f 0.03,0.01s t . It is seen that optimized topologies at two scales for = f 0.05s t (Figure 3b) and , substantial porous material is placed near the free edge this beam, which produces inertial force to offset the short impulsive loading, which i subsequently verified with the dynamic response as plotted in Figure 4. Also, the op mized macrostructure links the large mass distributed towards the free edge to the br ing ends by two horizontal beam-like members, which contribute to reducing the verti bending of the beam.   attributes to the fact that the resulting inertial force is sufficiently large to drive the be

Clamped Beam Design under Half-Cycle Cosine Load
In this example, we consider a beam fixed at both ends and excited via a concentrat half-cosine load vertically at the center of the bottom edge, as shown in Figure 6. T macroscopic design domain has length L = 12 m, height H = 2 m and thickness h = 0.01 We adopt a linear elastic material with a Young's modulus of 200 GPa, Poisson's ratio 0.3 and mass density of 7800 kg/m 3 . The macroscopic design domain and the unit cell a discretized by respective 5000 and 2500 bilinear square elements. The volume fracti limits for both macrostructure and unit cell are set to be 0.5. The Rayleigh damping p rameters α r and β r are the same as those in the first example. To solve this proble we adopt the input parameters listed in Table 3.

Clamped Beam Design under Half-Cycle Cosine Load
In this example, we consider a beam fixed at both ends and excited via a concentrated half-cosine load vertically at the center of the bottom edge, as shown in Figure 6. The macroscopic design domain has length L = 12 m, height H = 2 m and thickness h = 0.01 m. We adopt a linear elastic material with a Young's modulus of 200 GPa, Poisson's ratio of 0.3 and mass density of 7800 kg/m 3 . The macroscopic design domain and the unit cell are discretized by respective 5000 and 2500 bilinear square elements. The volume fraction limits for both macrostructure and unit cell are set to be 0.5. The Rayleigh damping parameters α r and β r are the same as those in the first example. To solve this problem, we adopt the input parameters listed in Table 3.

Clamped Beam Design under Half-Cycle Cosine Load
In this example, we consider a beam fixed at both ends and excited via a concentrated half-cosine load vertically at the center of the bottom edge, as shown in Figure 6. The macroscopic design domain has length L = 12 m, height H = 2 m and thickness h = 0.01 m. We adopt a linear elastic material with a Young's modulus of 200 GPa, Poisson's ratio of 0.3 and mass density of 7800 kg/m 3 . The macroscopic design domain and the unit cell are discretized by respective 5000 and 2500 bilinear square elements. The volume fraction limits for both macrostructure and unit cell are set to be 0.5. The Rayleigh damping parameters α r and β r are the same as those in the first example. To solve this problem, we adopt the input parameters listed in Table 3.

Clamped Beam Design under Half-Cycle Cosine Load
In this example, we consider a beam fixed at both ends and excited via a concentrat half-cosine load vertically at the center of the bottom edge, as shown in Figure 6. T macroscopic design domain has length L = 12 m, height H = 2 m and thickness h = 0.01 We adopt a linear elastic material with a Young's modulus of 200 GPa, Poisson's ratio 0.3 and mass density of 7800 kg/m 3 . The macroscopic design domain and the unit cell a discretized by respective 5000 and 2500 bilinear square elements. The volume fracti limits for both macrostructure and unit cell are set to be 0.5. The Rayleigh damping p rameters α r and β r are the same as those in the first example. To solve this proble we adopt the input parameters listed in Table 3.    Table 4 compares the relative errors of the sensitivity obtained through the discretizethen-differentiate AVM with those of the sensitivity obtained with the differentiate-thendiscretize AVM for this clamped beam problem for a load application duration of t f = 0.5 s. The former achieves consistent sensitivities with the FDM, while the latter results in obvious consistency errors. This comparison affirms the efficiency of the discretize-thendifferentiate AVM for dynamic problems in the time domain. To verify the discretize-thendifferentiate AVM for transient concurrent topology optimization, we apply this approach to solve the clamped beam problem and carry out a comparison of the optimized solution with those obtained via the differentiate-then-discretize AVM. These results, as illustrated in Figure 7, reveal that concurrent topology optimization based on the discretize-thendifferentiate AVM is more efficient for the transient problem due to lower optimum value of the dynamic compliance. croscopic element thickness 0.01 m mber of elements discretized in macroscopic design domain 5000 mber of elements discretized in microscopic design domain 2500 Table 4 compares the relative errors of the sensitivity obtained through the discretiz then-differentiate AVM with those of the sensitivity obtained with the differentiate-the discretize AVM for this clamped beam problem for a load application duration = f 0.5s t . The former achieves consistent sensitivities with the FDM, while the latter sults in obvious consistency errors. This comparison affirms the efficiency of the disc tize-then-differentiate AVM for dynamic problems in the time domain. To verify the d cretize-then-differentiate AVM for transient concurrent topology optimization, we app this approach to solve the clamped beam problem and carry out a comparison of the o timized solution with those obtained via the differentiate-then-discretize AVM. These sults, as illustrated in Figure 7, reveal that concurrent topology optimization based on t discretize-then-differentiate AVM is more efficient for the transient problem due to low optimum value of the dynamic compliance.    Figure 8 shows the convergence history and the optimal design for t f = 0.05 s. As seen from the results in Figures 7b and 8, the optimal topologies are highly dependent on t f . For short-term dynamic load, the optimizer assigns less porous material within the neighborhood of load application point and instead adds two beam-like members. That is favorable to endure the increased local deflection near the load application point, which arise as a result of the augmentation of dynamic influence. Figure 8 shows the convergence history and the optimal design for = f 0.05s t . A seen from the results in Figures 7b and 8, the optimal topologies are highly dependent o f t . For short-term dynamic load, the optimizer assigns less porous material within t neighborhood of load application point and instead adds two beam-like members. That favorable to endure the increased local deflection near the load application point, whi arise as a result of the augmentation of dynamic influence.      neighborhood of load application point and instead adds two beam-like members. That favorable to endure the increased local deflection near the load application point, whi arise as a result of the augmentation of dynamic influence.

Clamped Beam Design under Half-Cycle Cosine Load
In this example, we consider a beam fixed at both ends and excited via a concentrated half-cosine load vertically at the center of the bottom edge, as shown in Figure 6. The macroscopic design domain has length L = 12 m, height H = 2 m and thickness h = 0.01 m. We adopt a linear elastic material with a Young's modulus of 200 GPa, Poisson's ratio of 0.3 and mass density of 7800 kg/m 3 . The macroscopic design domain and the unit cell are discretized by respective 5000 and 2500 bilinear square elements. The volume fraction limits for both macrostructure and unit cell are set to be 0.5. The Rayleigh damping parameters α r and β r are the same as those in the first example. To solve this problem, we adopt the input parameters listed in Table 3.

Clamped Beam Design under Half-Cycle Cosine Load
In this example, we consider a beam fixed at both ends and excited via a concentrated half-cosine load vertically at the center of the bottom edge, as shown in Figure 6. The macroscopic design domain has length L = 12 m, height H = 2 m and thickness h = 0.01 m. We adopt a linear elastic material with a Young's modulus of 200 GPa, Poisson's ratio of 0.3 and mass density of 7800 kg/m 3 . The macroscopic design domain and the unit cell are discretized by respective 5000 and 2500 bilinear square elements. The volume fraction limits for both macrostructure and unit cell are set to be 0.5. The Rayleigh damping parameters α r and β r are the same as those in the first example. To solve this problem, we adopt the input parameters listed in Table 3.

Support Structure Design under Rotating Load
As shown in Figure 10, we use a square structure fixed at th to a rotating load with a specified constant amplitude and angula of upper free edge. The square domain has length L = 3 m and base material has a Young's modulus of 70 GPa, a Poisson's ratio of 7800 kg/m 3 . The macroscopic design domain and the unit cell a tive 5000 and 2500 bilinear square elements. The volume fraction ture and the unit cell are defined as 0.3 and 0.5, respectively. Th rameters α r and β r are assumed to be 50 s −1 and 3 × 10 −5 s, re problem, we adopt the input parameters listed in Table 5.

Support Structure Design under Rotating Load
As shown in Figure 10, we use a square structure fixed at the bottom edge, subjected to a rotating load with a specified constant amplitude and angular frequency at the center of upper free edge. The square domain has length L = 3 m and thickness h = 0.05 m. The base material has a Young's modulus of 70 GPa, a Poisson's ratio of 0.3, and a mass density of 7800 kg/m 3 . The macroscopic design domain and the unit cell are discretized by respective 5000 and 2500 bilinear square elements. The volume fraction limits of the macrostructure and the unit cell are defined as 0.3 and 0.5, respectively. The Rayleigh damping parameters α r and β r are assumed to be 50 s −1 and 3 × 10 −5 s, respectively. To solve this problem, we adopt the input parameters listed in Table 5. Table 5. Input parameters used to solve the support structure problem.

Parameter Value
Simulation time 10π/ω, ω = 100π and 25π rad/s To demonstrate the consistency of adjoint sensitivity analysis for transient concurrent topology optimization, we plot the relative error of the two sensitivities obtained with both differentiate-then-discretize and discretize-then-differentiate through examining this support structure design under a rotating load. These results with angular frequency ω = 100π rad/s, as illustrated in Table 6, confirm that the latter can ensure consistent sensitivities despite more cumbersome implementation. In gradient-based topology optimization, an accurate sensitivity analysis is requisite for the exact optimal solution. As a consequence, the optimized design based on the discretize-then-differentiate approach is necessary to be more effective due to high accuracy in sensitivity computation. Figure 11 demonstrates that the objective function converges to the smaller value acquired with discretize-then-differentiate than the counterpart acquired via differentiate-then-discretize. As such, we prefer the former for a transient multi-scale topology optimization problem. In order to study the influence of angular frequency on the final design for this support structure, we present an additional optimal design for ω = 25π rad/s, as shown in Figure 12. The first design (Figure 11b) adds an extra lateral resistant system in its macroscopic topology to diminish the structural lateral motion, whereas the second (Figure 12) is just composed of two rod-like members in its macroscopic topology. These two designs have a similar microscopic topology. Table 6. Comparison of design sensitivity and optimum for the support structure problem.

Sensitivity Analysis Method
Peak Relative Error (%)  (Figure 11b) adds an extra lateral resistant system in its macrosco topology to diminish the structural lateral motion, whereas the second ( Figure 12) is j composed of two rod-like members in its macroscopic topology. These two designs ha a similar microscopic topology.   Figure 12. Iterative history (left) and optimized topologies obtained (right) for a support structure with ω = 25π rad/s. Figure 13 presents the time history of horizontal displacement at the load application point and transient dynamic compliance for the two optimum designs demonstrated in Figures 11b and 12. The results indicate that the dynamic effect happen through the initial time steps, followed by vibration attenuation owing to damping dissipation. Furthermore, as is expected, the optimal design obtained for ω = 100π rad/s produce the lower vibrational level than the counterpart obtained for ω = 25π rad/s due to the additional lateral resistant system. The results indicate that the dynamic effect happen through the init time steps, followed by vibration attenuation owing to damping dissipation. Furthermo as is expected, the optimal design obtained for ω π = 100 rad / s produce the lower vib tional level than the counterpart obtained for ω π = 25 rad / s due to the additional late resistant system. denotes the dynamic compliance and denotes the displacement along the rotating load, respectively.

Building Design under Ground Excitation
This example aims to design a building under a time-varying ground acceleration a sinusoidal function. Figure 14 states this optimization problem with the initial config ration, ground acceleration as well as specified volume constraint at two scales.

Clamped Beam Design under Half-Cycle Cosine Load
In this example, we consider a beam fixed at both ends and excited via a concentrated half-cosine load vertically at the center of the bottom edge, as shown in Figure 6. The macroscopic design domain has length L = 12 m, height H = 2 m and thickness h = 0.01 m. We adopt a linear elastic material with a Young's modulus of 200 GPa, Poisson's ratio of 0.3 and mass density of 7800 kg/m 3 . The macroscopic design domain and the unit cell are discretized by respective 5000 and 2500 bilinear square elements. The volume fraction limits for both macrostructure and unit cell are set to be 0.5. The Rayleigh damping parameters α r and β r are the same as those in the first example. To solve this problem, we adopt the input parameters listed in Table 3.

Clamped Beam Design under Half-Cycle Cosine Load
In this example, we consider a beam fixed at both ends and excited via a concentrated half-cosine load vertically at the center of the bottom edge, as shown in Figure 6. The macroscopic design domain has length L = 12 m, height H = 2 m and thickness h = 0.01 m. We adopt a linear elastic material with a Young's modulus of 200 GPa, Poisson's ratio of 0.3 and mass density of 7800 kg/m 3 . The macroscopic design domain and the unit cell are discretized by respective 5000 and 2500 bilinear square elements. The volume fraction limits for both macrostructure and unit cell are set to be 0.5. The Rayleigh damping parameters α r and β r are the same as those in the first example. To solve this problem, we adopt the input parameters listed in Table 3.

Building Design under Ground Excitation
This example aims to design a building under a time-varying ground acceleration in a sinusoidal function. Figure 14 states this optimization problem with the initial configuration, ground acceleration as well as specified volume constraint at two scales. The building with length L = 75 m, height H = 75 m and thickness h = 0.05 m is clamped at the bottom and a lumped mass mc at the center of top edge is placed. The Young's modulus, Poisson's ratio and mass density of the base material are 35 Gpa, 0.25 and 2400 kg/m 3 , respectively. The macroscopic design domain and the unit cell are meshed into respective 10,000 and 2500 square bilateral elements, where volume fraction limits at the two scales are prescribed as 0.5. The Rayleigh damping parameters α r and β r are assumed to be 2 s −1 and 2 × 10 −6 s, respectively. Note that when considering ground accelerations, we replace the external load f with −m c a g I in Equation (14). In this example, the frequency of ground acceleration is supposed to be 2.5π rad/s. To solve this problem, we adopt the input parameters listed in Table 7.
ratio and mass density of the base material are 35 Gpa, 0.25 and 2400 kg/m 3 , resp The macroscopic design domain and the unit cell are meshed into respective 10 2500 square bilateral elements, where volume fraction limits at the two scales scribed as 0.5. The Rayleigh damping parameters α r and β r are assumed to be 2 × 10 −6 s, respectively. Note that when considering ground accelerations, we rep external load f with − c g m a I in Equation (14). In this example, the frequency o acceleration is supposed to be π 2.5 rad/s . To solve this problem, we adopt the i rameters listed in Table 7. Similarly, we first review the consistency of adjoint sensitivity analysis for transient concurrent topology optimization and then demonstrate the influence of sensitivity approximation on the final topology with this building design. These results in sensitivity calculation with a lumped mass m c = 0.3 × 10 6 kg, as listed in Table 8, indicating that the discretize-then-differentiate AVM can present consistent sensitivity due to high accuracy in nature. However, the differentiate-then-discretize AVM inherently generates inconsistent sensitivities. Consequently, we can obtain a more efficient multi-scale topology optimized via discretize-then-differentiate, as demonstrated in Figure 15. Figure 15 shows the optimal topologies obtained for 2.5π rad/s and m c = 0.6 × 10 6 kg. As is seen from the results in Figures 15b and 16, the optimum design is greatly susceptible to the lumped mass magnitude. The cross bars conjoined to the lumped mass are slightly thicker with increasing lumped mass. This is due to the larger inertial loads transferred from the lumped mass to the building when m c is increasing. Additionally, merely a lateral resistant system develops on the upper end of the building for small m c in Figure 16a, while an additive lateral resistant system develops at the bottom for large m c in Figures 15b and 16b, which is in favor of incremental inertial forces' transfer to the supports.  To comprehend the dynamic behavior of the building underground excitation along both horizontal and vertical directions, we plot the dynamic response of the optimum designs for various lumped mass, as illustrated in Figure 17. As observed from these results, the vertical displacement at the load application point is much larger than the horizontal counterpart due to the lateral resistant system regardless of the magnitude of the lumped mass. Note that with increasing lumped mass, the resultant vertical displacements at the load application point increase in the amplitude, such that the corresponding dynamic compliances became slightly larger. This inertial effect obviously influences the optimal topology, which cannot be apprehended with static optimization formulations.   To comprehend the dynamic behavior of the building underground excitation alo both horizontal and vertical directions, we plot the dynamic response of the optim designs for various lumped mass, as illustrated in Figure 17. As observed from these

Simply Supported 3D Structure
This example optimizes a 3D structure to examine the capability of the presented algorithm for large-scale transient topology optimization. As shown in Figure 18, this design domain has the following dimensions: length L = 4.5 m, height H = 0.75 m and thickness h = 0.5 m. This structure is supported at the bottom four corners under the same transient load as the first example. The Young's modulus, Poisson's ratio and mass density of base material are 200 Gpa, 0.3 and 7800 kg/m 3 , respectively. The macroscopic design domain is discretized with 13,500 eight-nodes brick elements and the unit cell with 8000 eight-nodes brick elements. The volume fraction limits at the two scales are prescribed as 0.5. The Rayleigh damping parameters α r and β r are assumed to be 10 s −1 and 2 × 10 −5 s, respectively. Table 9 offers all the adopted input data to solve the problem.
sults, the vertical displacement at the load application point is much larger than the ho zontal counterpart due to the lateral resistant system regardless of the magnitude of t lumped mass. Note that with increasing lumped mass, the resultant vertical displa ments at the load application point increase in the amplitude, such that the correspondi dynamic compliances became slightly larger. This inertial effect obviously influences t optimal topology, which cannot be apprehended with static optimization formulations

Simply Supported 3D Structure
This example optimizes a 3D structure to examine the capability of the present algorithm for large-scale transient topology optimization. As shown in Figure 18, this d sign domain has the following dimensions: length L = 4.5 m, height H = 0.75 m and thic ness h = 0.5 m. This structure is supported at the bottom four corners under the sam transient load as the first example. The Young's modulus, Poisson's ratio and mass de sity of base material are 200 Gpa, 0.3 and 7800 kg/m 3 , respectively. The macroscopic desi domain is discretized with 13,500 eight-nodes brick elements and the unit cell with 80 eight-nodes brick elements. The volume fraction limits at the two scales are prescribed 0.5. The Rayleigh damping parameters α r and β r are assumed to be 10 s −1 and 2 × 1 s, respectively. Table 9 offers all the adopted input data to solve the problem.     Figure 19 depicts the final designs at macro/micro scales. Compared with the 2D structure, the design space is enlarged by incorporating more freedom and a hollow pattern is generated in the middle domain. For a unit cell, the main microscopic structural members have coincident orientations with the corresponding macroscopic structural counterparts. This is favorable to transfer the load from the loading point to the constrained points. This numerical result demonstrates that the proposed approach has the potential to handle the optimization problem of 3D structures. In the future work, a fully parallelized MPI framework for multi-scale transient topology optimization is proposed to efficiently solve the large-scale transient lattice optimization problems on the basis of [34].   Figure 19 depicts the final designs at macro/micro scales. Compared with the 2D structure, the design space is enlarged by incorporating more freedom and a hollow pat tern is generated in the middle domain. For a unit cell, the main microscopic structura members have coincident orientations with the corresponding macroscopic structura counterparts. This is favorable to transfer the load from the loading point to the con strained points. This numerical result demonstrates that the proposed approach has th potential to handle the optimization problem of 3D structures. In the future work, a fully parallelized MPI framework for multi-scale transient topology optimization is proposed to efficiently solve the large-scale transient lattice optimization problems on the basis o [34].

Conclusions
This paper develops an efficient concurrent topological design approach for improv ing the dynamic performance of composite structures. According to the homogenized properties calculated via EBHM, the multi-scale dynamic finite element analysis is accom plished in the composite structure subjected to an impact load with the HHT-α method Figure 19. Optimized macroscale (left) and microscale (right) designs for simply supported 3D structure.

Conclusions
This paper develops an efficient concurrent topological design approach for improving the dynamic performance of composite structures. According to the homogenized properties calculated via EBHM, the multi-scale dynamic finite element analysis is accomplished in the composite structure subjected to an impact load with the HHT-α method. Two adjoint sensitivity analysis schemes, differentiate-then-discretize and discretize-then-differentiate, are developed to evaluate the derivatives of dynamic responses regarding design variables at two scales. The consistency errors in the sensitivity calculations obtained from both adjoint sensitivity analysis schemes are compared to analyze how the inconsistent sensitivities influence the optimal solution for linear structural dynamic problems.
The popular AVM based on the differentiate-then-discretize approach encounters significant consistency errors in the sensitivity evaluation as demonstrated using the numerical examples. Alternatively, the discretize-then-differentiate AVM tackles this inconsistent sensitivity problem and achieves the effective optimal solution, whereby the multi-scale topology optimization problems associated with transient response are efficiently resolved. We consider arbitrary loading situations with varying amplitudes, directions, and application durations besides ground acceleration, such that the proposed approach can resolve a wide variety of transient concurrent topology optimization problems. It is noted that the in-ertial force can play a significant role in the final optimal design at both macrostructure and microstructure levels, particularly when the composite structure suffers from the impact load imposed at a fast rate of speed. In future work, we extend the proposed concurrent topology optimization formulation to multi-material design of composite structures with non-uniform microstructures at macro and micro levels. Furthermore, the clustering-based approach grouping the microscopic unit cells based on a physical quantity, is introduced to implement the multi-scale topology optimization for a considerable reduction in computational cost.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to lab privacy.

Conflicts of Interest:
The authors declare no conflict of interest.