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

Abstract

This paper presents a concurrent topology optimization of multi-scale composite structures subjected to general time-dependent loads for minimizing dynamic compliance. A three-field 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 finite 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 efficient 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.

1. 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., differentiate-then-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 large-scale transient topology optimization employing the fully discretized adjoint sensitivity analysis in [35]. Nevertheless, to the author’s knowledge, very few investigations on multi-scale 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.

2. 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.

2.1. 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 mesh-independence, the original design variables are regulated with a smooth regularization filter [38] and expressed as follows:

$${\overline{\xi }}_i=\frac{{\displaystyle {\sum }_{k\in {\Phi }_i}{w}_{ki}^{\mathrm{mac}}{v}_k^{\mathrm{mac}}{\xi }_k}}{{\displaystyle {\sum }_{k\in {\Phi }_i}{w}_{ki}^{\mathrm{mac}}{v}_k^{\mathrm{mac}}}}$$

(1)

$${\overline{\eta }}_j=\frac{{\displaystyle {\sum }_{l\in {\Psi }_j}{w}_{lj}^{\mathrm{mic}}{v}_l^{\mathrm{mic}}{\eta }_l}}{{\displaystyle {\sum }_{l\in {\Psi }_j}{w}_{lj}^{\mathrm{mic}}{v}_l^{\mathrm{mic}}}}$$

(2)

where ${\Phi }_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 ${\Psi }_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_k^{\mathrm{mac}}$ is the volume of element $k$ in the macroscopic design domain and $v_l^{\mathrm{mic}}$ is the volume of element $l$ in the unit cell. The weighting factors $w_{ki}^{\mathrm{mac}}$ and $w_{lj}^{\mathrm{mic}}$ are defined using a linearly decaying function:

$$w_{ki}^{\mathrm{mac}}=R-\Vert x_k-x_i\Vert$$

(3)

$$w_{lj}^{\mathrm{mic}}=r-\Vert y_l-y_j\Vert$$

(4)

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:

$${\tilde{\xi }}_i={\xi }_{\min }+\left(1-{\xi }_{\min }\right)\frac{\tanh \left({\beta }^{\mathrm{mac}}{\xi }_{\mathrm{th}}\right)+\tanh \left({\beta }^{\mathrm{mac}}\left({\overline{\xi }}_i-{\xi }_{\mathrm{th}}\right)\right)}{\tanh \left({\beta }^{\mathrm{mac}}{\xi }_{\mathrm{th}}\right)+\tanh \left({\beta }^{\mathrm{mac}}\left(1-{\xi }_{\mathrm{th}}\right)\right)}$$

(5)

$${\tilde{\eta }}_j={\eta }_{\min }+\left(1-{\eta }_{\min }\right)\frac{\tanh \left({\beta }^{\mathrm{mic}}{\eta }_{\mathrm{th}}\right)+\tanh \left({\beta }^{\mathrm{mic}}\left({\overline{\eta }}_j-{\eta }_{\mathrm{th}}\right)\right)}{\tanh \left({\beta }^{\mathrm{mic}}{\eta }_{\mathrm{th}}\right)+\tanh \left({\beta }^{\mathrm{mic}}\left(1-{\eta }_{\mathrm{th}}\right)\right)}$$

(6)

where the physical design variables, ${\tilde{\xi }}_i$ and ${\tilde{\eta }}_j$, use the Ersatz parameters much less than one denoted by ${\xi }_{\min }$ and ${\eta }_{\min }$, respectively, to inhibit numerical instabilities of the stiffness and mass matrices when $\overline{\xi }\to 0$ and $\overline{\eta }\to 0$. ${\beta }^{\mathrm{mac}}$ and ${\beta }^{\mathrm{mic}}$ are exploited to regulate the aggressiveness of the projection function. ${\xi }_{\mathrm{th}}$ and ${\eta }_{\mathrm{th}}$ are the threshold density specified as 0.5 in this work.

2.2. 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].

$$D_j^{\mathrm{mic}}={\tilde{\eta }}_j^p{D}^{\mathrm{B}}$$

(7)

$$D_i^{\mathrm{mac}}=g\left({\tilde{\xi }}_i\right)D^{\mathrm{H}}$$

(8)

where $D^{\mathrm{B}}$ is the elastic constitutive matrix of base material and $D^{\mathrm{H}}$ is the effective macroscopic constitutive matrix, which is computed as follows:

$$D^{\mathrm{H}}=\frac{1}{\left|{\Omega }_m\right|}{\displaystyle \underset{{\Omega }_m}{\int }{D}_j^{\mathrm{mic}}\left(I-b{u}_m\right)}\mathrm{d}{\Omega }_m$$

(9)

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 ${\Omega }_m$.

The resultant displacement matrix $u_m$ is obtained through resolving the following unit cell equilibrium problem with periodic boundary conditions:

$$k^{\mathrm{mic}}{u}_m={\displaystyle \underset{{\Omega }_m}{\int }{b}^{\mathrm{T}}{D}_j^{\mathrm{mic}}\mathrm{d}{\Omega }_m}$$

(10)

where the stiffness matrix $k^{\mathrm{mic}}$ is given by the following:

$$k^{\mathrm{mic}}={\displaystyle \underset{{\Omega }_m}{\int }{b}^{\mathrm{T}}{D}_j^{\mathrm{mic}}b\mathrm{d}{\Omega }_m}$$

(11)

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:

$$g\left({\tilde{\xi }}_i\right)=\left(15{\tilde{\xi }}_i^p+{\tilde{\xi }}_i\right)/16$$

(12)

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:

$${\rho }^{\mathrm{H}}=\frac{1}{\left|{\Omega }_m\right|}{\displaystyle \underset{{\Omega }_m}{\int }{\rho }^{\mathrm{B}}{\tilde{\eta }}_j}\mathrm{d}{\Omega }_{mj}$$

(13)

where ${\rho }^{\mathrm{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].

2.3. 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:

$$M{\ddot{u}}_t+C{\dot{u}}_t+K{u}_t=f_t\hspace{1em}(t=0,\cdots ,\overline{N})$$

(14)

where $M$, $C$, and $K$ represent the global mass, damping, and stiffness matrices, respectively. ${\ddot{u}}_t$, ${\dot{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. $\overline{N}$ is the number of analysis steps. The global mass and stiffness matrices are assembled using the penalized macroscopic element matrix:

$$K={\displaystyle \sum _{i=1}^{N^{\mathrm{mac}}}g\left({\tilde{\xi }}_i\right)k_i^0}$$

(15)

$$M={\displaystyle \sum _{i=1}^{N^{\mathrm{mac}}}{\tilde{\xi }}_i{m}_i^0}$$

(16)

where

$$k_i^0={\displaystyle {\int }_{{\Omega }_i}{B}^{\mathrm{T}}{D}^{\mathrm{H}}B}\mathrm{d}{\Omega }_i$$

(17)

$$m_i^0={\rho }^{\mathrm{H}}{\displaystyle {\int }_{{\Omega }_i}{N}^{\mathrm{T}}N}\mathrm{d}{\Omega }_i$$

(18)

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

$$C={\alpha }_{\mathrm{r}}M+{\beta }_{\mathrm{r}}K$$

(19)

where ${\alpha }_{\mathrm{r}}$ and ${\beta }_{\mathrm{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:

$$\begin{array}{ll}\underset{\xi ,\eta }{\min }\hfill & f\left(\xi ,\eta ,u(t)\right)={\displaystyle \sum _{t=0}^{\overline{N}}{f}_t^{\mathrm{T}}{u}_t}\hfill \\ s.t.\hfill & M{\ddot{u}}_t+C{\dot{u}}_t+K{u}_t=f_t(t=0,\cdots ,\overline{N})\hfill \\ \hfill & G_1=({\displaystyle \sum _{i=1}^{N^{\mathrm{mac}}}{\tilde{\xi }}_i{v}_i^{\mathrm{mac}}})/V^{\mathrm{mac}}-\varsigma \le 0\hfill \\ \hfill & G_2=({\displaystyle \sum _{j=1}^{N^{\mathrm{mic}}}{\tilde{\eta }}_j{v}_j^{\mathrm{mic}}})/V^{\mathrm{mic}}-\vartheta \le 0\hfill \\ \hfill & 0\le {\xi }_i\le 1,\hspace{1em}1\le i\le N^{\mathrm{mac}}\hfill \\ \hfill & 0\le {\eta }_j\le 1,\hspace{1em}1\le j\le N^{\mathrm{mic}}\hfill \end{array}$$

(20)

where $f\left(\xi ,\eta ,u(t)\right)$ is the concerned objective function, $V^{\mathrm{mac}}$ and $V^{\mathrm{mic}}$ are the respective volumes of macroscopic and microscopic design domains. $\varsigma$ and $\vartheta$ 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^{\mathrm{mac}}$ elements, while the unit cell is discretized into $N^{\mathrm{mic}}$ elements.

3. 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 $\alpha$ into the Newmark method to control the numerical damping. Accordingly, the motion Equation (14) representing the dynamic equilibrium is modified as follows:

$$\begin{array}{c}M{\ddot{u}}_t+(1-\alpha )C{\dot{u}}_t+\alpha C{\dot{u}}_{t-1}+(1-\alpha )K{u}_t+\alpha K{u}_{t-1}\\ =(1-\alpha )f_t+\alpha f_{t-1},\hspace{1em}t=1,\dots ,\overline{N}\end{array}$$

(21)

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:

$$u_t=u_{t-1}+\Delta t{\dot{u}}_{t-1}+\Delta t^2\left[\left(1/2-\beta \right){\ddot{u}}_{t-1}+\beta {\ddot{u}}_t\right]$$

(22)

$${\dot{u}}_t={\dot{u}}_{t-1}+\Delta t\left[\left(1-\gamma \right){\ddot{u}}_{t-1}+\gamma {\ddot{u}}_t\right]$$

(23)

where the Newmark parameters $\beta$ and $\gamma$ are constants which control the integration accuracy and stability, respectively, by satisfying the following relationship:

$$0\le \alpha \le 1/3,\hspace{1em}\beta =\left(1+{\alpha }^2\right)/4,\hspace{1em}\gamma =\left(1+2\alpha \right)/2$$

(24)

By substitution of Equations (22) and (23) into Equation (21), the time-discretized motion equation in residual form is derived as follows:

$$R_t=M_1{\ddot{u}}_t+M_0{\ddot{u}}_{t-1}+C_0{\dot{u}}_{t-1}+K{u}_{t-1}-(1-\alpha )f_t-\alpha f_{t-1}=0$$

(25)

where

$$M_1=M+(1-\alpha )\gamma \Delta tC+(1-\alpha )\beta \Delta t^{2}K$$

(26)

$$M_0=(1-\alpha )(1-\gamma )\Delta tC+(1-\alpha )\left(\frac{1}{2}-\beta \right)\Delta t^{2}K$$

(27)

$$C_0=C+(1-\alpha )\Delta tK$$

(28)

Following a standard HHT-α scheme, we can obtain the dynamic response at each time step. We resolve Equation (25) for ${\ddot{u}}_t$ and thereupon compute $u_t$ and ${\dot{u}}_t$ by applying the Newmark-β Formulas (22) and (23), respectively. As for ${\ddot{u}}_0$, by assuming ${\dot{u}}_0$ and $u_0$ to be design-independent, it can be computed using the following residual equation:

$$R_0=M{\ddot{u}}_0+C{\dot{u}}_0+K{u}_0-f_0=0$$

(29)

4. 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:

$$\frac{\partial f}{\partial {\xi }_i}={\displaystyle \sum _{k\in {\Phi }_i}\frac{\partial f}{\partial {\tilde{\xi }}_k}\frac{\partial {\tilde{\xi }}_k}{\partial {\overline{\xi }}_k}\frac{\partial {\overline{\xi }}_k}{\partial {\xi }_i}}$$

(30)

$$\frac{\partial G_1}{\partial {\xi }_i}={\displaystyle \sum _{k\in {\Phi }_i}\frac{\partial G_1}{\partial {\tilde{\xi }}_k}\frac{\partial {\tilde{\xi }}_k}{\partial {\overline{\xi }}_k}\frac{\partial {\overline{\xi }}_k}{\partial {\xi }_i}}$$

(31)

$$\frac{\partial f}{\partial {\eta }_i}={\displaystyle \sum _{l\in {\Psi }_j}\frac{\partial f}{\partial {\tilde{\eta }}_l}\frac{\partial {\tilde{\eta }}_l}{\partial {\overline{\eta }}_l}\frac{\partial {\overline{\eta }}_l}{\partial {\eta }_j}}$$

(32)

$$\frac{\partial G_2}{\partial {\eta }_i}={\displaystyle \sum _{l\in {\Psi }_j}\frac{\partial G_2}{\partial {\tilde{\eta }}_l}\frac{\partial {\tilde{\eta }}_l}{\partial {\overline{\eta }}_l}\frac{\partial {\overline{\eta }}_l}{\partial {\eta }_j}}$$

(33)

where

$$\frac{\partial {\tilde{\xi }}_k}{\partial {\overline{\xi }}_k}=\left(1-{\xi }_{\min }\right)\frac{{\beta }^{\mathrm{mac}}{(\mathrm{sech}({\beta }^{\mathrm{mac}}({\overline{\xi }}_k-{\xi }_{\mathrm{th}})))}^2}{\tanh ({\beta }^{\mathrm{mac}}{\xi }_{\mathrm{th}})+\tanh ({\beta }^{\mathrm{mac}}(1-{\xi }_{\mathrm{th}}))}$$

(34)

$$\frac{\partial {\tilde{\eta }}_l}{\partial {\overline{\eta }}_l}=\left(1-{\eta }_{\min }\right)\frac{{\beta }^{\mathrm{mic}}{(\mathrm{sech}({\beta }^{\mathrm{mic}}({\overline{\eta }}_l-{\eta }_{\mathrm{th}})))}^2}{\tanh ({\beta }^{\mathrm{mic}}{\eta }_{\mathrm{th}})+\tanh ({\beta }^{\mathrm{mic}}(1-{\eta }_{\mathrm{th}}))}$$

(35)

$$\partial {\overline{\xi }}_k/\partial {\xi }_i=w_{ki}{v}_i^{\mathrm{mac}}/{\displaystyle {\sum }_{i\in {\Phi }_k}{w}_{ki}{v}_i^{\mathrm{mac}}}$$

(36)

$$\partial {\overline{\eta }}_l/\partial {\eta }_j=w_{lj}{v}_j^{\mathrm{mic}}/{\displaystyle {\sum }_{j\in {\Psi }_l}{w}_{lj}{v}_j^{\mathrm{mic}}}$$

(37)

The sensitivity of $f$ with respect to the arbitrary design variable $x$(${\xi }_i$,${\eta }_i$) is also written as follows:

$$\frac{\mathrm{d}f}{\mathrm{d}x}=\frac{\partial f}{\partial x}+{\displaystyle \sum _{t=0}^{\overline{N}}\frac{\partial f}{\partial u_t}\frac{\partial u_t}{\partial x}}$$

(38)

In order to facilitate the sensitivity analysis, we transform Equations (22) and (23) into the following residual form:

$$P_t=-u_t+u_{t-1}+\Delta t{\dot{u}}_{t-1}+\Delta t^2[(\frac{1}{2}-\beta ){\ddot{u}}_{t-1}+\beta {\ddot{u}}_t]=0t=1,\dots ,\overline{N}$$

(39)

$$Q_t=-{\dot{u}}_t+{\dot{u}}_{t-1}+\Delta t[(1-\gamma ){\ddot{u}}_{t-1}+\gamma {\ddot{u}}_t]=0t=1,\dots ,\overline{N}$$

(40)

Sequentially, we add adjoint variables ${\mathsf{\lambda }}_t$, ${\mathsf{\mu }}_t$ and ${\mathsf{\zeta }}_t$ and rewrite Equation (38) as follows:

$$\frac{\mathrm{d}f}{\mathrm{d}x}=\frac{\partial f}{\partial x}+{\displaystyle \sum _{t=0}^{\overline{N}}\frac{\partial f}{\partial u_t}\frac{\partial u_t}{\partial x}}+{\displaystyle \sum _{t=0}^{\overline{N}}{\lambda }_t^{\mathrm{T}}\frac{\mathrm{d}{R}_t}{\mathrm{d}x}}+{\displaystyle \sum _{t=1}^{\overline{N}}{\mu }_t^{\mathrm{T}}\frac{\mathrm{d}{P}_t}{\mathrm{d}x}}+{\displaystyle \sum _{t=1}^{\overline{N}}{\zeta }_t^{\mathrm{T}}\frac{\mathrm{d}{Q}_t}{\mathrm{d}x}}$$

(41)

From Equations (39) and (40), it is obvious that $\partial P_t/\partial x=0$ and $\partial Q_t/\partial x=0$. Due to the design-independence of the initial conditions, $\partial u_0/\partial x=0$ and $\partial {\dot{u}}_0/\partial x=0$. We employ these simplifications and eliminate all implicit terms including $\partial u/\partial x$, $\partial \dot{u}/\partial x$ and $\partial \ddot{u}/\partial x$ in Equation (41), such that the following adjoint equations can be obtained:

$${\mathsf{\lambda }}_0^{\mathrm{T}}\frac{\partial R_0}{\partial {\ddot{u}}_0}+{\mathsf{\lambda }}_1^{\mathrm{T}}\frac{\partial R_1}{\partial {\ddot{u}}_0}+{\mathsf{\mu }}_1^{\mathrm{T}}\frac{\partial P_1}{\partial {\ddot{u}}_0}+{\mathsf{\zeta }}_1^{\mathrm{T}}\frac{\partial Q_1}{\partial {\ddot{u}}_0}=0$$

(42)

$$\{\begin{array}{l}{\displaystyle \sum _{\ell =1}^{\overline{N}}\left({\lambda }_{\ell }^{\mathrm{T}}\frac{\partial R_{\ell }}{\partial u_t}+{\mu }_{\ell }^{\mathrm{T}}\frac{\partial P_{\ell }}{\partial u_t}+{\zeta }_{\ell }^{\mathrm{T}}\frac{\partial Q_{\ell }}{\partial u_t}+\frac{\partial f}{\partial u_t}\right)=0}\\ {\displaystyle \sum _{\ell =1}^{\overline{N}}\left({\lambda }_{\ell }^{\mathrm{T}}\frac{\partial R_{\ell }}{\partial {\dot{u}}_t}+{\mu }_{\ell }^{\mathrm{T}}\frac{\partial P_{\ell }}{\partial {\dot{u}}_t}+{\zeta }_{\ell }^{\mathrm{T}}\frac{\partial Q_{\ell }}{\partial {\dot{u}}_t}\right)=0}\\ {\displaystyle \sum _{\ell =1}^{\overline{N}}\left({\lambda }_{\ell }^{\mathrm{T}}\frac{\partial R_{\ell }}{\partial {\ddot{u}}_t}+{\mu }_{\ell }^{\mathrm{T}}\frac{\partial P_{\ell }}{\partial {\ddot{u}}_t}+{\zeta }_{\ell }^{\mathrm{T}}\frac{\partial Q_{\ell }}{\partial {\ddot{u}}_t}\right)=0}\end{array}$$

(43)

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:

$${\mathsf{\mu }}_{\overline{N}}=\frac{\partial f}{\partial u_{\overline{N}}},\hspace{1em}{\mathsf{\zeta }}_{\overline{N}}=0,\hspace{1em}{M}_1{\mathsf{\lambda }}_{\overline{N}}=-\beta \Delta t^2{\mathsf{\mu }}_{\overline{N}}-\gamma \Delta t{\mathsf{\zeta }}_{\overline{N}}$$

(44)

$${\mathsf{\mu }}_{t-1}=\frac{\partial f}{\partial u_{t-1}}+K{\mathsf{\lambda }}_t+{\mathsf{\mu }}_t,\hspace{1em}{\mathsf{\zeta }}_{t-1}=C_0{\mathsf{\lambda }}_t+\Delta t{\mathsf{\mu }}_t+{\mathsf{\zeta }}_t$$

(45)

$$M_1{\mathsf{\lambda }}_{t-1}=M_0{\mathsf{\lambda }}_t-\Delta t^2\left[\beta {\mu }_{t-1}+\left(\frac{1}{2}-\beta \right){\mu }_t\right]-\Delta t\left[\gamma {\zeta }_{t-1}+\left(1-\gamma \right){\zeta }_t\right]$$

(46)

$$M{\mathsf{\lambda }}_0=M_0{\mathsf{\lambda }}_1-\left(\frac{1}{2}-\beta \right)\Delta t^2{\mathsf{\mu }}_1-\left(1-\gamma \right)\Delta t{\mathsf{\zeta }}_1$$

(47)

Using the adjoint solution from Equations (44)–(47), we rewrite Equation (38) as follows:

$$\frac{\mathrm{d}f}{\mathrm{d}x}=\frac{\partial f}{\partial x}+{\displaystyle \sum _{t=0}^{\overline{N}}{\lambda }_t^{\mathrm{T}}\frac{\partial R_t}{\partial x}}$$

(48)

4.1. 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=\left[E_i\right]=\left[g\left({\tilde{\xi }}_i\right)\right]$, and the volume interpolation function, $V=\left[V_i\right]=\left[{\tilde{\xi }}_i\right]$, it facilitates recasting the sensitivity information of macroscopic design variables according to these fields. Therefore, we compute the sensitivity of $f$ with respect to $\xi$ by chain rule as follows:

$$\frac{\mathrm{d}f}{\mathrm{d}\xi }=\frac{\mathrm{d}f}{\mathrm{d}E}\frac{\partial E}{\partial \xi }+\frac{\mathrm{d}f}{\mathrm{d}V}\frac{\partial V}{\partial \xi }$$

(49)

where the sensitivities of $f$ regarding the macroscopic element volume fractions and stiffness parameters can be attained as demonstrated in Equation (48).

$$\mathrm{d}f/\mathrm{d}{E}_i=\partial f/\partial E_i+{\displaystyle \sum _{t=0}^{\overline{N}}{\lambda }_t^{\mathrm{T}}\partial R_t/\partial E_i}$$

(50)

$$\mathrm{d}f/\mathrm{d}{V}_i=\partial f/\partial V_i+{\displaystyle \sum _{t=0}^{\overline{N}}{\lambda }_t^{\mathrm{T}}\partial R_t/\partial V_i}$$

(51)

The terms, $\partial R_t/\partial E_i$ and $\partial R_t/\partial V_i$, are evaluated in terms of Equations (25) and (29), respectively. There is a case for $t=0$.

$$\frac{\partial R_t}{\partial E_i}=\frac{\partial K}{\partial E_i}\left(u_0+{\beta }_{\mathrm{r}}{\dot{u}}_0\right)=k_i\left(u_{i,0}+{\beta }_{\mathrm{r}}{\dot{u}}_{i,0}\right)$$

(52)

$$\frac{\partial R_t}{\partial V_i}=\frac{\partial M}{\partial E_i}\left({\ddot{u}}_0+{\alpha }_{\mathrm{r}}{\dot{u}}_0\right)=m_i\left({\ddot{u}}_{i,0}+{\alpha }_{\mathrm{r}}{\dot{u}}_{i,0}\right)$$

(53)

and for $t=1,\cdots ,\overline{N}$

$$\begin{array}{ll}\frac{\partial R_t}{\partial E_i}& =\frac{\partial K}{\partial E_i}\left[\left(1-\alpha \right)\left(u_t+{\beta }_{\mathrm{r}}{\dot{u}}_t\right)+\alpha \left(u_{t-1}+{\beta }_r{\dot{u}}_{t-1}\right)\right]\\ & =k_i\left[\left(1-\alpha \right)\left(u_{i,t}+{\beta }_{\mathrm{r}}{\dot{u}}_{i,t}\right)+\alpha \left(u_{i,t-1}+{\beta }_r{\dot{u}}_{i,t-1}\right)\right]\end{array}$$

(54)

$$\begin{array}{ll}\frac{\partial R_t}{\partial V_i}& =\frac{\partial M}{\partial E_i}\left[{\ddot{u}}_t+{\alpha }_{\mathrm{r}}\left(\left(1-\alpha \right){\dot{u}}_t+\alpha {\dot{u}}_{t-1}\right)\right]\\ & =m_i\left[{\ddot{u}}_{i,t}+{\alpha }_{\mathrm{r}}\left(\left(1-\alpha \right){\dot{u}}_{i,t}+\alpha {\dot{u}}_{i,t-1}\right)\right]\end{array}$$

(55)

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 ${\tilde{\xi }}_i$ is computed as follows:

$$\frac{\partial E_i}{\partial {\tilde{\xi }}_i}=\frac{1}{16}\left(15p{\tilde{\xi }}_i^{p-1}+1\right)$$

(56)

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).

4.2. 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:

$$\frac{\mathrm{d}f}{\mathrm{d}{\eta }_j}=\frac{\mathrm{d}f}{\mathrm{d}{D}^{\mathrm{H}}}\frac{\partial D^{\mathrm{H}}}{\partial {\eta }_j}+\frac{\mathrm{d}f}{\mathrm{d}{\rho }^{\mathrm{H}}}\frac{\partial {\rho }^{\mathrm{H}}}{\partial {\eta }_j}$$

(57)

where

$$\frac{\partial D^{\mathrm{H}}}{\partial {\tilde{\eta }}_j}=\frac{p{\tilde{\eta }}_j^{p-1}}{\left|{\Omega }_m\right|}{\displaystyle \underset{{\Omega }_m}{\int }{D}_j^{\mathrm{mic}}\left(I-b{u}_m\right)}\mathrm{d}{\Omega }_m$$

(58)

$$\frac{\partial {\rho }^{\mathrm{H}}}{\partial {\tilde{\eta }}_j}=\frac{1}{\left|{\Omega }_m\right|}{\displaystyle \underset{{\Omega }_m}{\int }{\rho }^{\mathrm{B}}}\mathrm{d}{\Omega }_{mj}$$

(59)

The sensitivity of $f$ with respect to the effective material properties can be attained from Equation (48), i.e.,

$$\frac{\mathrm{d}f}{\mathrm{d}{D}^{\mathrm{H}}}=\frac{\partial f}{\partial D^{\mathrm{H}}}+{\displaystyle \sum _{t=0}^{\overline{N}}{\lambda }_t^{\mathrm{T}}\frac{\partial R_t}{\partial D^{\mathrm{H}}}}$$

(60)

$$\frac{\mathrm{d}f}{\mathrm{d}{\rho }^{\mathrm{H}}}=\frac{\partial f}{\partial {\rho }^{\mathrm{H}}}+{\displaystyle \sum _{t=0}^{\overline{N}}{\lambda }_t^{\mathrm{T}}\frac{\partial R_t}{\partial {\rho }^{\mathrm{H}}}}$$

(61)

where $\partial f/\partial D^{\mathrm{H}}=0$ and $\partial f/\partial {\rho }^{\mathrm{H}}=0$, according to the objective function as shown in Equation (20).

Similarly, the partial derivatives, $\partial R_t/\partial D^{\mathrm{H}}$ and $\partial R_t/\partial {\rho }^{\mathrm{H}}$, are evaluated using Equations (25) and (29), and for $i=0$:

$$\frac{\partial R_0}{\partial D^{\mathrm{H}}}=\frac{\partial K}{\partial D^{\mathrm{H}}}\left(u_0+{\beta }_{\mathrm{r}}{\dot{u}}_0\right)={\displaystyle \sum _{i=1}^{N^{\mathrm{mac}}}g\left({\tilde{\xi }}_i\right)\frac{\partial k_i^0}{\partial D^{\mathrm{H}}}}\left(u_0+{\beta }_{\mathrm{r}}{\dot{u}}_0\right)$$

(62)

$$\frac{\partial R_0}{\partial {\rho }^{\mathrm{H}}}=\frac{\partial M}{\partial {\rho }^{\mathrm{H}}}\left({\ddot{u}}_0+{\alpha }_{\mathrm{r}}{\dot{u}}_0\right)={\displaystyle \sum _{i=1}^{N^{\mathrm{mac}}}{\tilde{\xi }}_i\frac{\partial m_i^0}{\partial {\rho }^{\mathrm{H}}}}\left({\ddot{u}}_0+{\alpha }_{\mathrm{r}}{\dot{u}}_0\right)$$

(63)

for $t=1,\cdots ,\overline{N}$:

$$\begin{array}{ll}\frac{\partial R_t}{\partial D^{\mathrm{H}}}& =\frac{\partial K}{\partial D^{\mathrm{H}}}\left[\left(1-\alpha \right)\left(u_t+{\beta }_{\mathrm{r}}{\dot{u}}_t\right)+\alpha \left(u_{t-1}+{\beta }_{\mathrm{r}}{\dot{u}}_{t-1}\right)\right]\\ & ={\displaystyle \sum _{i=1}^{N^{\mathrm{mac}}}g\left({\tilde{\xi }}_i\right)\frac{\partial k_i^0}{\partial D^{\mathrm{H}}}}\left[\left(1-\alpha \right)\left(u_t+{\beta }_{\mathrm{r}}{\dot{u}}_t\right)+\alpha \left(u_{t-1}+{\beta }_{\mathrm{r}}{\dot{u}}_{t-1}\right)\right]\end{array}$$

(64)

$$\frac{\partial R_t}{\partial {\rho }^{\mathrm{H}}}=\frac{\partial M}{\partial {\rho }^{\mathrm{H}}}\left[{\ddot{u}}_t+{\alpha }_{\mathrm{r}}\left(\left(1-\alpha \right){\dot{u}}_t+\alpha {\dot{u}}_{t-1}\right)\right]={\displaystyle \sum _{i=1}^{N^{\mathrm{mac}}}{\tilde{\xi }}_i\frac{\partial m_i^0}{\partial {\rho }^{\mathrm{H}}}}\left[{\ddot{u}}_t+{\alpha }_{\mathrm{r}}\left(\left(1-\alpha \right){\dot{u}}_t+\alpha {\dot{u}}_{t-1}\right)\right]$$

(65)

4.3. Solution Procedure

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

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.

5. Adjoint Sensitivity Analysis Using Differentiate-Then-Discretize

The differentiate-then-discretize AVM constructs the adjoint equation in a semi-discretized 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 $\Phi$ in the following integral form:

$$\Phi ={\displaystyle {\int }_0^{J}c\left(u,\dot{u}\right)}\mathrm{d}\overline{t}$$

(66)

where $J$ is the duration of the dynamic event and $\overline{t}$ is the continuous time variable.

We introduce the motion Equation (14) into $\Phi$ and thereby obtain the sensitivity ${\Phi }^{\prime }$ by standard AVM:

$${\Phi }^{\prime }={\displaystyle {\int }_0^J\left(\partial c/\partial u{u}^{\prime }+\partial c/\partial \dot{u}{\dot{u}}^{\prime }\right)\mathrm{d}\overline{t}}+{\displaystyle {\int }_0^J{\lambda }^{\mathrm{T}}{\left(M\ddot{u}+C\dot{u}+Ku-f\right)}^{\prime }\mathrm{d}\overline{t}}$$

(67)

where the prime denotes differentiation regarding the design variables and $\lambda$ denotes the smooth adjoint variable. Through twice integrating-by-parts, we rearrange ${\Phi }^{\prime }$ as follows:

$$\begin{array}{ll}{\Phi }^{\prime }=& {\displaystyle {\int }_0^J{\lambda }^{\mathrm{T}}\left(M^{\prime }\ddot{u}+C^{\prime }\dot{u}+K^{\prime }u\right)\mathrm{d}\overline{t}}\\ & +{\displaystyle {\int }_0^J{u^{\prime }}^{\mathrm{T}}\left(M\ddot{\lambda }-C\dot{\lambda }+K\lambda +\partial c/\partial u-\mathrm{d}\left(\partial c/\partial \dot{u}\right)/\mathrm{d}t\right)\mathrm{d}\overline{t}}\\ & +{\left[{u^{\prime }}^{\mathrm{T}}\left(-M\dot{\lambda }+C\lambda +{\partial c/\partial \dot{u}}^{\mathrm{T}}\right)+{{\dot{u}}^{\prime }}^{\mathrm{T}}M\lambda \right]|}_{\overline{t}=J}\end{array}$$

(68)

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^{\prime }\left(J\right)$ and ${\dot{u}}^{\prime }\left(J\right)$ at the final time step, we assign the adjoint variables such that the terminal conditions are satisfied as follows:

$$\lambda \left(J\right)=0,M\dot{\lambda }\left(J\right)={{\left(\partial c/\partial u\right)}^{\mathrm{T}}|}_{\overline{t}=J}$$

(69)

To transform the adjoint problem into the initial value problem, we use a variable transformation $\overline{t}=\tau \left(s\right)=J-s$ and then construct a composite function $\Lambda$ satisfying $\Lambda \left(s\right)=\lambda \left(\tau \left(s\right)\right)$. Accordingly, we rewrite Equation (68) by transforming all the terms including $u^{\prime }$ and ${\dot{u}}^{\prime }$:

$$\begin{array}{ll}{\Phi }^{\prime }=& {\displaystyle {\int }_0^J{\lambda }^{\mathrm{T}}\left(M^{\prime }\ddot{u}+C^{\prime }\dot{u}+K^{\prime }u\right)\mathrm{d}\overline{t}+{\displaystyle {\int }_0^J{u^{\prime }}^{\mathrm{T}}\left(\tau \left(s\right)\right){\left(\partial c/\partial u-\mathrm{d}\left(\partial c/\partial \dot{u}\right)/\mathrm{d}t\right)|}_{\overline{t}=\tau \left(s\right)}\mathrm{d}s}}\\ & +{\displaystyle {\int }_0^{J}u{\prime }^{\mathrm{T}}\left(\tau \left(s\right)\right)\left(M\ddot{\Lambda }\left(\tau \left(s\right)\right)+C\dot{\Lambda }\left(\tau \left(s\right)\right)+K\Lambda \left(\tau \left(s\right)\right)\right)\mathrm{d}s}\\ & +\left[{u^{\prime }}^{\mathrm{T}}\left(J\right)\left(M\dot{\Lambda }\left(0\right)-C\Lambda \left(0\right)+{{\partial c/\partial \dot{u}}^{\mathrm{T}}|}_{\overline{t}=J}\right)+{{\dot{u}}^{\prime }}^{\mathrm{T}}\left(J\right)M\Lambda \left(0\right)\right]\end{array}$$

(70)

To annul all the terms containing $u^{\prime }$ and ${\dot{u}}^{\prime }$, we formulate the adjoint variable $\Lambda$ as follows:

$$\begin{array}{l}M\ddot{\Lambda }\left(\tau \left(s\right)\right)+C\dot{\Lambda }\left(\tau \left(s\right)\right)+K\Lambda \left(\tau \left(s\right)\right)={\left(-\partial c/\partial u+\mathrm{d}\left(\partial c/\partial \dot{u}\right)/\mathrm{d}t\right)|}_{\overline{t}=\tau \left(s\right)}\\ \Lambda \left(0\right)=0,M\dot{\Lambda }\left(0\right)=-\partial c/\partial \dot{u}{{}^{\mathrm{T}}|}_{\overline{t}=J}\end{array}$$

(71)

where the sensitivity is simplified as follows:

$${\Phi }^{\prime }={\displaystyle {\int }_0^J{\Lambda }^{\mathrm{T}}\left(J-\overline{t}\right)\left(M^{\prime }\ddot{u}\left(\overline{t}\right)+C^{\prime }\dot{u}\left(\overline{t}\right)+K^{\prime }u\left(\overline{t}\right)\right)\mathrm{d}\overline{t}}={{\Lambda }^{\mathrm{T}}\ast \left(M^{\prime }\ddot{u}+C^{\prime }\dot{u}+K^{\prime }u\right)|}_{\overline{t}=J}$$

(72)

where $\ast$ denotes the convolution operator.

Following the obtained displacement and velocity, $u_t$ and ${\dot{u}}_t$, we approximate the original objective function employing the rectangular formula:

$$\tilde{\Phi }={\displaystyle \sum _{t=0}^{\overline{N}}c\left(u_t,{\dot{u}}_t\right)}$$

(73)

Based on the discretized adjoint variables ${\mathsf{\Lambda }}_n$ solution from Equation (71), the sensitivity of objective function is approximated as follows:

$${\tilde{\Phi }}^{\prime }={\displaystyle \sum _{t=0}^{\overline{N}}{\Lambda }_{\overline{N}-t}^{\mathrm{T}}\left(M^{\prime }{\ddot{u}}_t+C^{\prime }{\dot{u}}_t+K^{\prime }{u}_t\right)}$$

(74)

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).

6. 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 $\alpha =0.05$ and determine the Newmark constants $\beta$ and $\gamma$ by employing the formulas $\beta ={\left(1+\alpha \right)}^2/4$ and $\gamma =\left(1+2\alpha \right)/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.

6.1. 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³. The Rayleigh damping parameters ${\alpha }_{\mathrm{r}}$ and ${\beta }_{\mathrm{r}}$ are assumed to be 10 s⁻¹ and 1 × 10⁻⁵ 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. Input parameters used to solve the cantilever beam problem.

Parameter	Value
Simulation time	0.05 s, 0.03 s, and 0.01 s
Number of time steps	100
Young’s modulus of base material	200 GPa
Poisson’s ratio of base material	0.3
Mass density of base material	7800 kg/m3
Rayleigh damping parameters	10 s−1 and 1 × 10−5 s
Volume fraction limit of macrostructure and unit cell	0.5 and 0.5
Filter radius in macro/micro design domain and filter exponent	[0.12, 0.002] and 3
Chosen element type	Four-nodes bilateral element
Macroscopic element thickness	0.01 m
Number of elements discretized in macroscopic design domain	5000
Number of elements discretized in microscopic design domain	2500

.

Table 2. Comparison of design sensitivity and optimum for the cantilever beam problem.

Sensitivity Analysis Method	Peak Relative Error (%)		Optimum (Nm)
	Macro Design Domain	Micro Design Domain
Discretize-then-differentiate	1.8	1.6	1.12
Differentiate-then-discretize	13.6	11.5	1.71

compares the design sensitivity between two AVM approaches, discretize-then-differentiate and differentiate-then-discretize, and FDM for this cantilever beam problem with a load application duration of $t_{\mathrm{f}}=0.05\mathrm{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.

Figure 4 shows the iterative histories during concurrent optimization and the resulting optimal designs for the load application duration of $t_{\mathrm{f}}=0.03,0.01\mathrm{s}$. It is seen that the optimized topologies at two scales for $t_{\mathrm{f}}=0.05\mathrm{s}$ (Figure 3b) and $t_{\mathrm{f}}=0.03\mathrm{s}$ (Figure 4a) are nearly identical to each other. But in the case of $t_{\mathrm{f}}=0.01\mathrm{s}$ (Figure 4b), the optimal configurations greatly distinguish from the counterparts obtained for larger load application duration. When $t_{\mathrm{f}}=0.01\mathrm{s}$, 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_{\mathrm{f}}=0.05\mathrm{s}$ resemble those for $t_{\mathrm{f}}=0.03\mathrm{s}$, while they are significantly different for $t_{\mathrm{f}}=0.01\mathrm{s}$. Particularly for $t_{\mathrm{f}}=0.01\mathrm{s}$, the optimal beam continuously deflects downward, although the external load gradually decreases following $t_{\mathrm{f}}=0.005\mathrm{s}$, which attributes to the fact that the resulting inertial force is sufficiently large to drive the beam downward.

6.2. 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³. 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 ${\alpha }_{\mathrm{r}}$ and ${\beta }_{\mathrm{r}}$ are the same as those in the first example. To solve this problem, we adopt the input parameters listed in

Table 3. Input parameters used to solve the clamped beam problem.

Parameter	Value
Simulation time	0.5 s and 0.05 s
Number of time steps	100
Young’s modulus of base material	200 GPa
Poisson’s ratio of base material	0.3
Mass density of base material	7800 kg/m3
Rayleigh damping parameters	10 s−1 and 1 × 10−5 s
Volume fraction limit of macrostructure and unit cell	0.5 and 0.5
Filter radius in macro/micro design domain and filter exponent	[0.10, 0.005] and 3
Chosen element type	Four-nodes bilateral element
Macroscopic element thickness	0.01 m
Number of elements discretized in macroscopic design domain	5000
Number of elements discretized in microscopic design domain	2500

.

Table 4. Comparison of design sensitivity and optimum for the clamped beam problem.

Sensitivity Analysis Method	Peak Relative Error (%)		Optimum (N m)
	Macro Design Domain	Micro Design Domain
Discretize-then-differentiate	2.1	1.8	0.46
Differentiate-then-discretize	18.9	16.3	0.82

compares the relative errors of the sensitivity obtained through the discretize-then-differentiate AVM with those of the sensitivity obtained with the differentiate-then-discretize AVM for this clamped beam problem for a load application duration of $t_{\mathrm{f}}=0.5\mathrm{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-then-differentiate AVM for dynamic problems in the time domain. To verify the discretize-then-differentiate 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-then-differentiate AVM is more efficient for the transient problem due to lower optimum value of the dynamic compliance.

Figure 8 shows the convergence history and the optimal design for $t_{\mathrm{f}}=0.05\mathrm{s}$. As seen from the results in Figure 7b and Figure 8, the optimal topologies are highly dependent on $t_{\mathrm{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 9 depicts the dynamic response of respective optimized design for $t_{\mathrm{f}}=0.5\mathrm{s}$ and $t_{\mathrm{f}}=0.05\mathrm{s}$, as demonstrated in Figure 7b and Figure 8. The damping effect is obviously identified from the results acquired for $t_{\mathrm{f}}=0.5\mathrm{s}$, when the amplitude of vertical displacement and transient dynamic compliance decay over the load application duration due to the energy dissipation in the damping material. In contrast to the results for $t_{\mathrm{f}}=0.5\mathrm{s}$, the dissipation effect of damping attenuates over time for $t_{\mathrm{f}}=0.05\mathrm{s}$, owing to the shorter load application duration. Therefore, the load application duration directly affects the dissipation of internal energy and the structural vibration. This explains why the optimum designs are susceptible to the load application duration.

6.3. 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³. 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 ${\alpha }_{\mathrm{r}}$ and ${\beta }_{\mathrm{r}}$ are assumed to be 50 s⁻¹ and 3 × 10⁻⁵ s, respectively. To solve this problem, we adopt the input parameters listed in

Table 5. Input parameters used to solve the support structure problem.

Parameter	Value
Simulation time	10π/ω, ω = 100π and 25π rad/s
Number of time steps	100
Young’s modulus of base material	70 GPa
Poisson’s ratio of base material	0.3
Mass density of base material	7800 kg/m3
Rayleigh damping parameters	50 s−1 and 3 × 10−5 s
Volume fraction limit of macrostructure and unit cell	0.3 and 0.5
Filter radius in macro/micro design domain and filter exponent	[0.06, 0.0015] and 3
Chosen element type	Four-nodes bilateral element
Macroscopic element thickness	0.05 m
Number of elements discretized in macroscopic design domain	5000
Number of elements discretized in microscopic design domain	2500

.

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 $\omega =100\pi \mathrm{rad}/\mathrm{s}$, as illustrated in

Table 6. Comparison of design sensitivity and optimum for the support structure problem.

Sensitivity Analysis Method	Peak Relative Error (%)		Optimum (Nm)
	Macro Design Domain	Micro Design Domain
Discretize-then-differentiate	1.4	0.9	1.96
Differentiate-then-discretize	15.6	14.2	2.91

, 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 $\omega =25\pi \mathrm{rad}/\mathrm{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.

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 Figure 11b and Figure 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 $\omega =100\pi \mathrm{rad}/\mathrm{s}$ produce the lower vibrational level than the counterpart obtained for $\omega =25\pi \mathrm{rad}/\mathrm{s}$ due to the additional lateral resistant system.

6.4. 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³, 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 ${\alpha }_{\mathrm{r}}$ and ${\beta }_{\mathrm{r}}$ are assumed to be 2 s⁻¹ and 2 × 10⁻⁶ 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\pi \mathrm{rad}/\mathrm{s}$. To solve this problem, we adopt the input parameters listed in

Table 7. Input parameters used to solve the building problem.

Parameter	Value
Simulation time	4.8 s
Number of time steps	100
Young’s modulus of base material	35 GPa
Poisson’s ratio of base material	0.25
Mass density of base material	2400 kg/m3
Rayleigh damping parameters	2 s−1 and 2 × 10−6 s
Lumped masses	0.1 × 106, 0.3 × 106 and 0.6 × 106 kg
Volume fraction limit of macrostructure and unit cell	0.3 and 0.5
Filter radius in macro/micro design domain and filter exponent	[1.0, 0.02] and 3
Chosen element type	Four-nodes bilateral element
Macroscopic element thickness	0.05 m
Number of elements discretized in macroscopic design domain	10,000
Number of elements discretized in microscopic design domain	2500

.

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\times {10}^6\mathrm{kg}$, as listed in

Table 8. Comparison of design sensitivity and optimum for the building problem.

Sensitivity Analysis Method	Peak Relative Error (%)		Optimum (Nm)
	Macro Design Domain	Micro Design Domain
Discretize-then-differentiate	2.0	1.2	23.4
Differentiate-then-discretize	23.6	15.4	31.6

, 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\pi \mathrm{rad}/\mathrm{s}$ and $m_c=0.6\times {10}^6\mathrm{kg}$. As is seen from the results in Figure 15b and Figure 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 Figure 15b and Figure 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.

6.5. 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³, 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 ${\alpha }_{\mathrm{r}}$ and ${\beta }_{\mathrm{r}}$ are assumed to be 10 s⁻¹ and 2 × 10⁻⁵ s, respectively.

Table 9. Input parameters used to solve the simply supported 3D structure problem.

Parameter	Value
Simulation time	0.05 s
Number of time steps	100
Young’s modulus of base material	200 GPa
Poisson’s ratio of base material	0.3
Mass density of base material	7800 kg/m3
Rayleigh damping parameters	10 s−1 and 2 × 10−5 s
Volume fraction limit of macrostructure and unit cell	0.5 and 0.5
Filter radius in macro/micro design domain and filter exponent	[0.15, 0.008] and 3
Chosen element type	Eight-nodes brick element
Number of elements discretized in macroscopic design domain	13,500
Number of elements discretized in microscopic design domain	8000

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].

7. 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 inertial 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.