Topology optimization of multi-story buildings under fully non-stationary stochastic seismic ground motion

Abstract

Topology optimization has been mainly addressed for structures under static loads using a deterministic setting. Nonetheless, many structural systems are subjected to uncertain dynamic loads, and thus efficient approaches are required to evaluate the optimal topology in such kind of applications. Within this framework, the present paper deals with the topology optimization of multi-story buildings subjected to seismic ground motion. Because of the inherent randomness of the earthquakes, the uncertain system response is determined through a random vibration-based approach in which the seismic ground motion is described as filtered white Gaussian noise with time-varying amplitude and frequency content (i.e., fully non-stationary seismic ground motion). The paper is especially concerned with the assessment of the dynamic response sensitivity for the gradient-based numerical solution of the optimization problem. To this end, an approximated construction of the gradient is proposed in which explicit, exact derivatives with respect to the design variables are computed analytically through direct differentiation for a sub-assembly of elements (up to a single element) resulting from the discretization of the optimizable domain. The proposed strategy is first validated for the simpler case of stationary base excitation by comparing the results with those obtained using an exact approach based on the adjoint method, and its correctness is ultimately verified for the more general case of non-stationary seismic ground motion. Overall, this validation demonstrates that the proposed approach leads to accurate results at low computational effort. Further numerical investigations are finally presented to highlight to what extent the features of the non-stationary seismic ground motion influence the optimal topology.

1 Introduction

Automated structural optimization procedures have proven to be a valuable support in real-life structural design in order to enhance the performance of the built environment while reducing its environmental and economic impacts (Lagaros 2018; Aldwaik and Adeli 2014). In this context, topology optimization is nowadays widely recognized as an effective way for finding the optimal material distribution (i.e., location and connectivity of structural members) within the prescribed design space. It can be implemented conveniently for identifying the optimal load-resisting system that ensures the highest structural performance at reduced material consumption. Several topology optimization applications have been presented until now considering two-dimensional (Sigmund 2001; Stromberg et al. 2012) or three-dimensional (Liu and Tovar 2014; Angelucci et al. 2021c) domains under static loads. On the other hand, the search for the optimal topology under dynamic loading has been primarily addressed using a deterministic setting by means of system frequencies optimization (Ma et al. 1994; Du and Olhoff 2007; Maeda et al. 2006). Alternative approaches have been also explored to take into account explicitly the system dynamic response, either via modal decomposition (Zhao and Wang 2016; Angelucci et al. 2020; Martin and Deierlein 2020) or direct integration (Kaveh et al. 2012; Chahande and Arora 1994; Min et al. 1999; Balling et al. 2009; Allahdadian and Boroomand 2016).

It is worth noting that a deterministic structural design optimization has some limitations when dealing with dynamic loads affected by inherent uncertainties, such as in civil engineering structures subjected to wind loads and/or earthquake. Accounting for uncertain dynamic loads in topology optimization, however, poses additional significant computational challenges that ultimately makes a crude use of Monte Carlo simulations unfeasible. For such applications, the random vibration theory provides a more efficient approach to describe the dynamic response of structural systems under stochastic excitation. Along this line, Zhang et al. (2015) developed a computational strategy that merges the pseudo-excitation method and the mode acceleration method in order to optimize the topology of large-scale structures subjected to stationary random excitation. Chun et al. (2016) formulated a reliability-based design optimization problem using a discrete representation of the stochastic excitation. The focus of such work is the topology optimization of linear structures subjected to a stationary random process, which allows the use of a closed-form solution for the estimation of the probability of failure that quantifies the system reliability (in detail, the probability of failure is taken equal to the instantaneous failure probability at the point-in-time corresponding to the onset of the stationarity response). Time domain, frequency domain, and state-space methods have been compared by Zhu et al. (2017) when looking for the optimal lateral load-resisting system that minimizes the variance of frame systems under stationary stochastic ground motion. Gomez and Spencer (2019) addressed the topology optimization of structural systems under stationary stochastic excitation via random vibration theory. They developed an exact, yet very efficient, approach based on the adjoint method to facilitate the gradient-based numerical solution of the corresponding optimization problem, see also (Gomez et al. 2020, 2021). An energy-based design criterion has been considered by Angelucci et al. (2021a) for the topology optimization of linear elastic multi-story buildings subjected to earthquake, under the simplifying assumption of stationary stochastic seismic ground motion. Instead of exploiting the random vibration theory, Suksuwan and Spence (2018) explored a different way for optimizing the topology of structural systems under wind and seismic loads, in which the original probabilistic and dynamic problem is converted into an equivalent static and deterministic one.

Most of the researches that employ the random vibration theory in topology optimization problems presume a stationary excitation, but such simplification does not allow to reflect the time-dependent variation of amplitude and frequency content in some dynamic loads. In civil engineering applications, for instance, ordinary wind loading can be often represented as stationary process, but seismic loading is essentially a non-stationary phenomenon. Actually, the seismic design of the structures should take into proper account that both intensity and frequency content of the earthquake are time-dependent (Basone et al. 2017). The time-varying intensity of the seismic ground motion (i.e., non-stationarity in amplitude) typically manifests with an initial built-up phase, a gradually decaying tail portion at the end, and a fairly uniform intense strong-motion phase in-between. On the other hand, the time-varying frequency content (i.e., non-stationarity in frequency content) is related to the phase difference of various seismic waves (Ohsaki 1979). The role of a non-stationary seismic ground motion on the structural response has been extensively studied in previous researches. For instance, Jangid (2004) demonstrated that the root mean square of the displacements undergone by a linear elastic single-degree-of-freedom system subjected to stochastic seismic ground motion depends significantly on the shape of the amplitude modulation of the base excitation. The influence of the temporal non-stationarity in the frequency content of the seismic ground motion on the response of linear elastic single- and multi-degree-of freedom systems was later on demonstrated by Li et al. (2016), who concluded that a fully non-stationary earthquake model (i.e., a seismic ground motion with time-varying amplitude and frequency content) is needed for rational and reliable earthquake-resistant design. Hence, non-stationary loading conditions to which structural systems can be exposed during their lifetime, such as seismic loading, have to be considered properly for the optimization of their topology.

Although handling with non-stationary excitation is of utmost importance in some applications, such dynamic load case is seldom considered in the current literature about topology optimization. A somewhat pertinent attempt about that has been made by Xu et al. (2017), who addressed the optimal stiffness distribution in multi-story buildings under seismic ground motion. As a matter of fact, however, the study by Xu et al. (2017) resembles a size optimization problem, which is addressed by simulating the earthquake as amplitude-modulated and filtered white Gaussian noise (the frequency content of the seismic ground motion is instead assumed time-independent). After having pointed out that the stochastic structural response differs depending on whether a stationary or amplitude-modulated excitation is considered, Xu et al. (2017) performed the design optimization by means of a gradient-free method. To the best of the authors’ knowledge, the only topology optimization approaches proposed so far accounting for a non-stationary excitation have been developed by Hu et al. (2018) and Li et al. (2020), in which an explicit time domain approach (Su and Xu 2014; Hu et al. 2016) has been implemented in order to find the optimal lateral bracing system of frame structures and the optimal bar elements for truss structures, respectively. It is too evident that research on topology optimization under non-stationary excitation is still lacking, and thus the present work is meant at providing a contribution in this field.

Specifically, the present paper illustrates a computational framework for the topology optimization of the lateral load-resisting systems in multi-story buildings subjected to earthquake. The optimized topology is found by discretizing the structural system with the standard finite element method under the assumption of linear behavior and implementing the solid isotropic material with penalization approach. The seismic ground motion is herein modeled as fully non-stationary random process (i.e., filtered white Gaussian noise whose amplitude and frequency vary in time), unlike most studies in the field of topology optimization which instead consider a stationary dynamic loading (e.g., Gomez and Spencer 2019; Gomez et al. 2020). A gradient-based technique is adopted to solve the design optimization problem, which aims at minimizing the stochastic compliance of the structural system under a prescribed maximum volume. This, in turn, also calls for the sensitivity analysis of the objective function with respect to the topological design variables. Since accomplishing this task at global level requires prohibitive computational time and memory consumption, an approximated approach is herein introduced. Indeed, while the assumption of stationary excitation allows to calculate the exact gradient efficiently, for example by means of the adjoint method (e.g., Gomez and Spencer 2019; Gomez et al. 2020), an approximated way to compute the gradient is here proposed to deal with a non-stationary excitation. This relies on the concept of sub-assembly, that is a subset of the optimizable domain consisting of the finite element with respect to which the derivative is performed and a few finite elements within its neighborhood to improve the accuracy of the approximation if needed. Once the sub-assembly is defined, explicit, exact derivatives are computed analytically through direct differentiation of the related Lyapunov equation for non-stationary conditions considering a finite set of time steps. An extensive preliminary validation demonstrates accuracy, stability and efficiency of the proposed computational framework. The final numerical investigation highlights the practical application of the proposed approach as well as the role of the non-stationary seismic ground motion on the optimized topology of the lateral load-resisting system in multi-story buildings.

2 Formulation of the structural topology optimization problem

2.1 System description and design domain

The topology optimization of the lateral load-resisting system (LRS) of a multi-story building subjected to earthquake starts from the definition of the design domain. To this end, it can be noted that the region where the optimal material distribution has to be determined is naturally represented by the external envelope of the building. The continuous design domain $\Upomega$ so defined is here discretized using 4-node quadrilateral (Q4) elements as illustrated in Fig. 1.

Fig. 1

Structural system discretization and identification of the design domain

Together with the bracing system resulting from the optimization procedure, the LRS of multi-story buildings also consists of a secondary frame system $\Upgamma$ that withstands vertical and lateral loads, which is here discretized using Euler–Bernoulli beam elements. Past researches have highlighted the importance of including such a secondary system for bounding the continuous design domain (Stromberg et al. 2012; Bobby et al. 2014). On one hand, this frame allows to model properly the overall stiffness of the structural system. On the other hand, it is useful to prevent the direct application of external loads to the design domain, which may jeopardize the objectivity of the final results. The introduction of the secondary system is especially beneficial in case of topology optimization problems under assigned volume constraint. Indeed, since the presence of vertical members helps reducing the concentration of structural material along the boundaries of the domain, bracing layouts of the LRS are better defined while the identification of the working points is simplified. According to Stromberg et al. (2012), the boundary members are discretized into smaller elements with translational and rotational degrees of freedom, such that the end nodes of the frame elements are attached continuously at every node of the optimizable domain. The complete LRS is thus fully defined by combining the (optimizable) continuous design domain $\Upomega$ with the (not optimizable) secondary frame system $\Upgamma$ as shown in Fig. 1.

Each floor accommodates structural (i.e., slabs, floor framing) and non-structural (i.e., finishes, partition walls) elements that—although not part of the LRS—provide supplementary structural mass. This additional mass is estimated by considering the floor tributary areas. It generally represents the most of the mass in typical multi-story buildings while the contribution attributable to the LRS is almost negligible. As a consequence, it is reasonable to assume that the mass is primarily located at the floors level.

In the design of multi-story buildings, it is customary to consider that the floor system supporting the gravity loads behaves like a diaphragm rigid in its own plane. This legitimizes the reduction of the mass of each floor to a single point (i.e., the center of mass of the floor). However, in topology optimization problems, the introduction of the displacement constraints needed to simulate such a behavior may produce locations with unrealistically high axial stiffness that may undermine the objectivity of the final solution. Since the rigidity/flexibility of the floor diaphragms has a significant influence on the building behavior, the effects of the floor deformability cannot be neglected, and a reasonable stiffness must thus be assigned to each floor. Furthermore, the mass is distributed at a discrete number of points of the floor system, namely the master nodes, placed at the intersections with the secondary system as shown in Fig. 1. For further details about the influence of the mass modeling on the optimization of the LRS for multi-story buildings, the reader is referred to the work by Angelucci et al. (2021b). As already observed by Bobby et al. (2017), since all loads are applied at the master nodes only, the design domain topology does not depend on their particular arrangement as far as the secondary system (to which the master nodes belong to by definition) is not included in the optimizable domain. That is to say, so doing the system topology is insensitive with respect to the position of the nodes where loads apply, which is a desirable feature of the optimization process because their definition is arbitrary and should not influence the design.

2.2 Optimum structural design problem

The multi-story building is modeled as linear elastic system with linear viscous damping. This is a reasonable assumption for the optimization of buildings subjected to frequent and low intensity earthquakes (e.g., Elettore et al. 2020), especially when the design is targeted at ensuring the operativity and mitigating the damage. Additionally, the hypothesis of linear elastic behavior is fully compliant with most existing studies about topology optimization of multi-story buildings subjected to seismic stochastic excitation (e.g., Hu et al. 2018; Chun et al. 2016).

Let $N_f$ be the total number of degrees of freedom resulting from the finite element-based system discretization. The equation of motion under a horizontal seismic ground motion $\ddot{x}_g$ acting over the time window [0, T] reads:

$$\begin{aligned} \ddot{\mathbf{u}} + \mathbf{M}^{-1}\mathbf{C}\dot{\mathbf{u}}+\mathbf{M}^{-1}\mathbf{K}\mathbf{u}=-\mathbf{r}\ddot{x}_g , \end{aligned}$$

(1)

where $\mathbf{M}$, $\mathbf{C}$ and $\mathbf{K}$ are mass matrix, linear viscous damping matrix and linear elastic stiffness matrix, respectively. Moreover, $\mathbf{u}$ are the displacements of the structure whereas $\mathbf{r}$ is the incidence vector. Here and henceforth, the dependency from the time variable t is omitted for the sake of conciseness. The linear viscous damping matrix $\mathbf{C}$ is formulated in agreement with the Rayleigh’s model, and thus:

$$\begin{aligned} \mathbf{C}= a_0 \mathbf{M}+ a_1 \mathbf{K}, \end{aligned}$$

(2)

where $a_0$ and $a_1$ are two positive constants calculated by imposing that the viscous damping ratio $\xi _s$ is the same for the first two modes of the structure. Therefore, it is obtained:

$$\begin{aligned}&a_0 = \xi _s \dfrac{2\omega _1\omega _2}{\omega _1+\omega _2} , \end{aligned}$$

(3a)

$$\begin{aligned}&a_1 = \xi _s \dfrac{2}{\omega _1+\omega _2} , \end{aligned}$$

(3b)

where $\omega _1$ and $\omega _2$ are the first two circular natural frequencies of the structure.

The topology optimization problem for such multi-story building aims at minimizing the compliance under the seismic ground motion acting over the time window [0, T] while fulfilling the prescribed maximum volume $V_{\text {max}}$. Formally, the optimum design problem reads:

$$\begin{aligned} \begin{aligned} \min _{\varvec{\rho }}&\left\{ \int _0^{T}{ \mathbf{u}^\top \mathbf{K}\mathbf{u}\text{ d }t } \right\} \\&\text{ s.t. } \\&\sum _{e \in \Upomega } V_{\Upomega _e} = V_{\text {max}}\\&\varvec{\rho }_{\text {min}} \le \varvec{\rho }\le \varvec{\rho }_{\text {max}},\\ \end{aligned} \end{aligned}$$

(4)

where $\varvec{\rho }$ is vector collecting the design variables that define the topology of the structural system, with $\varvec{\rho }_{\text {min}}$ and $\varvec{\rho }_{\text {max}}$ being its lower and upper bound, respectively. Moreover, $V_{\Upomega _e}$ is the volume of the eth finite element within the optimizable domain $\Upomega$.

It is well-known that an inherent randomness exists in the occurrence of earthquakes, which reflects into significant uncertainty in predicting the seismic response of structural systems. In order to cope with the randomness of the seismic ground motion in structural design, the earthquake ground accelerations can be simulated conveniently as a stochastic process. This approach was developed through the decades starting from the pioneering work by Housner (1947), and it is herein employed because it allows to calculate efficiently the statistical dynamic response of structural systems under earthquake. Since the seismic ground motion is here modeled as random process, the compliance needs to be redefined in a stochastic sense. The optimization problem dealt with in the present work is thus obtained by revising Eq. (4) accordingly as follows:

$$\begin{aligned} \begin{aligned} \min _{\varvec{\rho }}&\left\{ \text {E} \left[ \int _0^{T}{ \mathbf{u}^\top \mathbf{K}\mathbf{u}\text{ d }t } \right] \right\} \\&\text{ s.t. } \\&\sum _{e \in \Upomega } V_{\Upomega _e} = V_{\text {max}}\\&\varvec{\rho }_{\text {min}} \le \varvec{\rho }\le \varvec{\rho }_{\text {max}},\\ \end{aligned} \end{aligned}$$

(5)

where $\text {E}[\cdot ]$ is the mean value operator. The objective function f of the optimization problem in Eq. (5) represents the expected compliance and it can be rewritten as follows:

$$\begin{aligned} f = \text {E} \left[ \int _0^{T}{ {\mathbf{u}}^\top \mathbf{K}{\mathbf{u}}\text{ d }t } \right] = \int _0^{T}{ \mathbf{r}^\top \left( \mathbf{K}\otimes \mathbf{R}_{ {\mathbf{u}} {\mathbf{u}}} \right) \mathbf{r}\text{ d }t } , \end{aligned}$$

(6)

where $\mathbf{R}_{ {\mathbf{u}} {\mathbf{u}}}$ is the covariance matrix of the system response and $\otimes$ is the element-wise multiplication operator.

3 Structural response under non-stationary stochastic seismic ground motion

3.1 Seismic ground motion modeling

The horizontal component of the seismic ground motion $\ddot{x}_g$ is modeled as fully non-stationary filtered white Gaussian noise. In the context of earthquake engineering, the Kanai–Tajimi filter and the Clough–Penzien filter are two commonly adopted linear filtering techniques. The Kanai–Tajimi model is the most simple one since it employs a single linear second-order filter. The Clough–Penzien filter is more appropriate for modeling an earthquake because an additional linear filter with respect to the Kanai–Tajimi model removes the embedded low-frequency components. Therefore, the Clough–Penzien filter model is here adopted. The assumption of fully non-stationary seismic ground motion implies that both amplitude and frequency content of the seismic ground motion change throughout the earthquake duration. With these premises, the random seismic ground motion is here defined as follows:

$$\begin{aligned}&\ddot{x}_g = \mathbf{a}_p^\top \mathbf{y}_p , \end{aligned}$$

(7a)

$$\begin{aligned}&\dot{\mathbf{y}}_p = \mathbf{D}_p \mathbf{y}_p + \mathbf{v}_p W , \end{aligned}$$

(7b)

where

$$\begin{aligned}&\mathbf{a}_p = \begin{Bmatrix} -\omega _p^2 \\ -2\xi _p\omega _p \\ \omega _k^2 \\ 2\xi _k\omega _k \end{Bmatrix} , \end{aligned}$$

(8a)

$$\begin{aligned}&\mathbf{y}_p = \begin{Bmatrix} x_p \\ \dot{x}_p \\ x_k \\ \dot{x}_k \end{Bmatrix} , \end{aligned}$$

(8b)

$$\begin{aligned}&\mathbf{D}_p = \begin{bmatrix} 0 &{} 1 &{} 0 &{} 0 \\ -\omega _p^2 &{} -2\xi _p\omega _p &{} \omega _k^2 &{} 2\xi _k\omega _k \\ 0 &{} 0 &{} 0 &{} 1 \\ 0 &{} 0 &{} -\omega _k^2 &{} -2\xi _k\omega _k \end{bmatrix} , \end{aligned}$$

(8c)

$$\begin{aligned}&\mathbf{v}_p = \begin{Bmatrix} 0 \\ 0 \\ 0 \\ -\varphi \end{Bmatrix} . \end{aligned}$$

(8d)

Herein, $\omega _k$ and $\xi _k$ are the dominant circular frequency and dominant damping of the soil, respectively. On the other hand, $\omega _p$ and $\xi _p$ denote the parameters of the filter hindering the low-frequency components of the dynamic excitation. Moreover, W is a zero-mean white Gaussian noise having constant power spectral density $S_0$ whereas $\varphi$ rules the amplitude modulation. Under the hypothesis of stationary seismic ground motion, all filter parameters (i.e., $\omega _k$, $\xi _k$, $\omega _p$ and $\xi _p$) are constants, and the parameter ruling the amplitude modulation (i.e., $\varphi$) is equal to one. Conversely, the assumption of fully non-stationary seismic ground motion here adopted implies that all these parameters, in general, can vary in time. Specifically, both $\omega _k$ and $\varphi$ are commonly assumed as time-dependent function in order to account for the temporal non-stationarity in frequency and amplitude of the seismic ground motion, respectively.

Several formulations that relate the peak ground acceleration (PGA) $\ddot{x}_g^{\text {max}}$ to $S_0$ have been proposed. Following Liu et al. (2016), the relationship used in this study is:

$$\begin{aligned} S_0 = \dfrac{\left( \ddot{x}_g^{\text {max}}\right) ^2}{\delta ^2 \left[ \pi {\omega }_k \left( 2{\xi }_k+\dfrac{1}{2{\xi }_k}\right) \right] } , \end{aligned}$$

(9)

where the peak factor $\delta$ is equal to 2.80. In case of fully non-stationary seismic ground motion, the constant value of $S_0$ is obtained from Eq. (9) by assuming a constant reference value for both dominant circular frequency ${\omega }_k$ and dominant damping of the soil ${\xi }_k$.

3.2 Non-stationary stochastic structural response analysis

By introducing the $N \times 1$ state-space vector (with $N=2N_f$):

$$\begin{aligned} \mathbf{y}_s = \begin{Bmatrix} \mathbf{u}\\ \dot{\mathbf{u}} \end{Bmatrix} , \end{aligned}$$

(10)

the equation of motion in Eq. (1) is rewritten as follows:

$$\begin{aligned} \dot{\mathbf{y}}_s = \mathbf{A}_s \mathbf{y}_s + \mathbf{H}_p \mathbf{y}_p, \end{aligned}$$

(11)

where

$$\begin{aligned} \mathbf{A}_s= & {} \begin{bmatrix} \varvec{0}_{N_f} &{} \mathbf{I}_{N_f} \\ -\mathbf{M}^{-1}\mathbf{K}&{} -\mathbf{M}^{-1}\mathbf{C}\end{bmatrix} , \end{aligned}$$

(12)

$${\mathbf{H}}_{p} = \left[ {\begin{array}{ll} {0_{{N_{f} \times 4}} } \\ { - {\mathbf{ra}}_{p}^{{ \top }} } \\ \end{array} } \right]$$

(13)

Herein, $\varvec{0}$ is a (square or rectangular) zero matrix or vector and $\mathbf{I}$ is the (square) identity matrix, respectively (their size is given in the corresponding subscripts, whereby a single or a double entry are used for a square matrix or rectangular matrix, respectively).

The state-space representation of the overall dynamics is obtained by assembling the filter equation in Eq. (7a, 7b) and the equation of motion in Eq. (11) as follows:

$$\begin{aligned} \underbrace{\begin{Bmatrix} \dot{\mathbf{y}}_s \\ \dot{\mathbf{y}}_p \end{Bmatrix}}_{\dot{\mathbf{y}}}= \underbrace{\begin{bmatrix} \mathbf{A}_s &{} \mathbf{H}_p \\ \varvec{0}_{4\times N} &{} \mathbf{D}_p \end{bmatrix}}_{\mathbf{A}} \underbrace{\begin{Bmatrix} \mathbf{y}_s \\ \mathbf{y}_p \end{Bmatrix}}_{\mathbf{y}} + \underbrace{\begin{Bmatrix} \varvec{0}_{N\times 1} \\ \mathbf{v}_pW \end{Bmatrix}}_{\mathbf{f}} . \end{aligned}$$

(14)

The corresponding covariance matrix is:

$$\begin{aligned} \mathbf{R}= \text {E}\left[ \mathbf{y}\mathbf{y}^\top \right] =\text {E} \left[ \begin{Bmatrix} \mathbf{y}_s \\ \mathbf{y}_p \end{Bmatrix} \begin{Bmatrix} \mathbf{y}_s&\mathbf{y}_p \end{Bmatrix} \right] = \begin{bmatrix} \mathbf{R}_{\mathbf{y}_s \mathbf{y}_s} &{} \mathbf{R}_{\mathbf{y}_s \mathbf{y}_p} \\ \mathbf{R}_{\mathbf{y}_p \mathbf{y}_s} &{} \mathbf{R}_{\mathbf{y}_p \mathbf{y}_p}\end{bmatrix} , \end{aligned}$$

(15)

in which

$$\begin{aligned} \mathbf{R}_{\mathbf{y}_s \mathbf{y}_s} = \begin{bmatrix} \mathbf{R}_{\mathbf{u}\mathbf{u}} &{} \mathbf{R}_{\mathbf{u}\dot{\mathbf{u}}} \\ \mathbf{R}_{\dot{\mathbf{u}} \mathbf{u}} &{} \mathbf{R}_{\dot{\mathbf{u}} \dot{\mathbf{u}}}\end{bmatrix} . \end{aligned}$$

(16)

The matrix $\mathbf{R}$, in turn, is the solution of the Lyapunov equation in non-stationary conditions, which reads:

$$\begin{aligned} \mathbf{A}\mathbf{R}+ \mathbf{R}\mathbf{A}^\top + \mathbf{B}= \dot{\mathbf{R}} , \end{aligned}$$

(17)

where $\mathbf{B}$ is a matrix whose elements are equal to zero expect that the element whose index is $(N+4,N+4)$, which is equal to $2\pi S_0 \varphi ^2$.

A discrete set of time instants labeled as $1,\ldots ,i,i+1,\ldots$ is considered within the time window [0, T]. Evenly spaced time instants with a constant time step $\Updelta t$ are assumed. If a linear variation of $\dot{\mathbf{R}}$ within each time interval $\Updelta t$ is supposed, then:

$$\begin{aligned} {^{i+1}}\mathbf{R}= {^{i}}\mathbf{R}+ \dfrac{1}{2}\Updelta t \left( {^{i+1}}\dot{\mathbf{R}} + {^{i}}\dot{\mathbf{R}} \right) , \end{aligned}$$

(18)

under the assumption that ${^{1}}\mathbf{R}= \varvec{0}_{N+4}$.

After some straightforward manipulations, the time discretization of the non-stationary version of the Lyapunov equation given by Eq. (17) based on Eq. (18) leads to the following stationary version of the Lyapunov equation at each time step:

$$\begin{aligned} {^{i+1}}\bar{\mathbf{A}} {^{i+1}}\mathbf{R}+ {^{i+1}}\mathbf{R}{^{i+1}}\bar{\mathbf{A}}^\top + {^{i+1}}\bar{\mathbf{B}} = \varvec{0}_{N+4} , \end{aligned}$$

(19)

where

$$\begin{aligned}&{^{i+1}}\bar{\mathbf{A}} = \dfrac{1}{2} \left( -\mathbf{I}_{N+4} + \Updelta t {^{i+1}}\mathbf{A}\right) , \end{aligned}$$

(20)

$$\begin{aligned}&{^{i+1}}\bar{\mathbf{B}} = \left( \mathbf{I}_{N+4} + \dfrac{1}{2}\Updelta t {^{i}}\mathbf{A}\right) {^{i}}\mathbf{R}+ \dfrac{1}{2}\Updelta t {^{i}}\mathbf{R}{^{i}}\mathbf{A}^\top \nonumber \\&\qquad\quad + \dfrac{1}{2}\Updelta t \left( {^{i+1}}\mathbf{B}+ {^{i}}\mathbf{B}\right) . \end{aligned}$$

(21)

This time discretization of the non-stationary Lyapunov equation given by Eqs. (19)–(21) is detailed in Appendix 1.

Once Eq. (19) is solved at each time instant, the objective function given by Eq. (6) can be calculated through a numerical approximation of the involved integral (i.e., via trapezoidal numerical integration).

4 Solution of the structural topology optimization problem

4.1 Optimal material distribution

The global mass matrix $\mathbf{M}$ can be expressed as the superposition of two contributions, i.e., the mass matrix referring to the design domain $\mathbf{M}_\Upomega$ and that of the secondary system including any carried mass $\mathbf{M}_\Upgamma$. Similarly, the global stiffness matrix $\mathbf{K}$ can be expressed as the superposition of two contributions, i.e., the stiffness matrix referring to the design domain $\mathbf{K}_\Upomega$ and that of the secondary system $\mathbf{K}_\Upgamma$. The solid isotropic material with penalization (SIMP) approach (Bendsøe and Sigmund 1999) is adopted in the present study to look for the optimal material distribution in the design domain. It assumes a continuous material distribution, by which mass and stiffness for the design domain are interpolated using a smooth convex function defined through a heuristic power-law as follows:

$$\begin{aligned}\mathbf{M}(\varvec{\rho })&= \mathbf{M}_\Upomega (\varvec{\rho }) + \mathbf{M}_\Upgamma =\bigwedge _{e \in \Upomega }{\mathbf{M}_{\Upomega _e}(\varvec{\rho })} + \mathbf{M}_\Upgamma \nonumber \\ &=\bigwedge _{e \in \Upomega }{\rho _e^q \mathbf{M}_{\Upomega _e}^0} + \mathbf{M}_\Upgamma , \end{aligned}$$

(22)

$$\begin{aligned}\mathbf{K}(\varvec{\rho })&= \mathbf{K}_\Upomega (\varvec{\rho }) + \mathbf{K}_\Upgamma =\bigwedge _{e \in \Upomega }{\mathbf{K}_{\Upomega _e}(\varvec{\rho })} + \mathbf{K}_\Upgamma \nonumber \\& =\bigwedge _{e \in \Upomega }{\rho _e^p \mathbf{K}_{\Upomega _e}^0} + \mathbf{K}_\Upgamma , \end{aligned}$$

(23)

where $\bigwedge$ stands for the standard finite element assembly operator, $\rho _e$ is the topological design variable of the eth finite element, whereas $\mathbf{M}_{\Upomega _e}$ and $\mathbf{K}_{\Upomega _e}$ are mass and stiffness matrix of the eth finite element within the optimizable domain $\Upomega$, respectively (e.g., mass and stiffness matrix of the Q4 finite element). The value of $\rho _e$ is constrained to vary between a small (non-zero) positive value for ensuring the positive definiteness of the global stiffness matrix and one. The penalty factors $q \ge 1$ and $p \ge 1$ are introduced to avoid undesirable fictitious intermediate densities in the relaxed setting and to steer the solution toward binary values (i.e., black-and-white image). Taking into account these penalty factors, it follows that $\mathbf{M}_{\Upomega _e}(\varvec{\rho })=\rho _e^q \mathbf{M}_{\Upomega _e}^0$ and $\mathbf{K}_{\Upomega _e}(\varvec{\rho })=\rho _e^p \mathbf{K}_{\Upomega _e}^0$ as in Eqs. (22) and (23), where $\mathbf{M}_{\Upomega _e}^0$ and $\mathbf{K}_{\Upomega _e}^0$ are mass and stiffness matrix of the eth finite element of the optimizable domain $\Upomega$ with solid material, respectively. It is here assumed $q=1$ and $p=3$, since it has been shown that these values generally provide superior convergence properties toward binary solutions. In order to cope with the non-convexity of the density-based topology optimization problem formulation and to stabilize the convergence, a continuation method is implemented as suggested by Bendsøe and Sigmund (2003). In detail, the penalization parameter p is slowly increased starting from $p=1$ up to $p=3$ for cautiously shifting from convexity to non-convexity. Finally, in order to avoid numerical instabilities (i.e., check-board pattern and mesh-dependency) as well as to impose minimum length control, a linear density filter is also employed using the H-filter operator (Bourdin 2001).

4.2 Sensitivity analysis

The solution of the minimum compliance problem given by Eq. (5) can be accomplished via a gradient-based numerical procedure. This calls for a sensitivity analysis of constraint and objective function with respect to the design variables. Since the constraint function is an explicit linear combination of design variables, its sensitivity analysis is straightforward. As regards the objective function, it is obtained from Eq. (6):

$$\begin{aligned} \dfrac{\partial f }{\partial \rho _e} = \int _0^{T}{ \mathbf{r}^\top \left( \dfrac{\partial \mathbf{K}}{\partial \rho _e} \otimes \mathbf{R}_{{\mathbf{u}}{\mathbf{u}}} \right) \mathbf{r}\text{ d }t } + \int _0^{T}{ \mathbf{r}^\top \left( \mathbf{K}\otimes \dfrac{\partial \mathbf{R}_{{\mathbf{u}}{\mathbf{u}}}}{\partial \rho _e} \right) \mathbf{r}\text{ d }t } . \end{aligned}$$

(24)

The derivative of $\mathbf{R}$ with respect to $\rho _e$—from which $\partial \mathbf{R}_{{\mathbf{u}}{\mathbf{u}}}/\partial \rho _e$ in Eq. (24) is extracted according to Eqs. (15) and (16)—can be determined from another stationary Lyapunov equation once Eq. (19) is solved. In fact, taking into account the time discretization, it is obtained:

$$\begin{aligned} {^{i+1}}\bar{\mathbf{A}} \dfrac{\partial {^{i+1}}\mathbf{R}}{\partial \rho _e} + \dfrac{\partial {^{i+1}}\mathbf{R}}{\partial \rho _e} {^{i+1}}\bar{\mathbf{A}}^\top + {^{i+1}}\bar{\bar{\mathbf{B}}} = \varvec{0}_{N+4} , \end{aligned}$$

(25)

where

$$\begin{aligned} {^{i+1}}\bar{\bar{\mathbf{B}}} = \dfrac{1}{2}\Updelta t \dfrac{\partial {^{i+1}}\mathbf{A}}{\partial \rho _e} {^{i+1}}\mathbf{R}+ \dfrac{1}{2}\Updelta t {^{i+1}}\mathbf{R}\dfrac{\partial {^{i+1}}\mathbf{A}^\top }{\partial \rho _e} + \frac{\partial {^{i+1}}\bar{\mathbf{B}}}{\partial \rho _e} . \end{aligned}$$

(26)

Herein, the term $\partial {^{i+1}}\bar{\mathbf{B}}/\partial \rho _e$ is obtained from Eq. (21) as follows:

$$\begin{aligned} \begin{aligned} \dfrac{\partial {^{i+1}}\bar{\mathbf{B}}}{\partial \rho _e} =&\dfrac{\partial {^{i}}\mathbf{R}}{\partial \rho _e} + \dfrac{1}{2}\Updelta t \dfrac{\partial {^{i}}\mathbf{A}}{\partial \rho _e} {^{i}}\mathbf{R}+ \dfrac{1}{2}\Updelta t {^{i}}\mathbf{A}\dfrac{\partial {^{i}}\mathbf{R}}{\partial \rho _e} \\&+ \dfrac{1}{2}\Updelta t \dfrac{\partial {^{i}}\mathbf{R}}{\partial \rho _e} {^{i}}\mathbf{A}^\top + \dfrac{1}{2}\Updelta t {^{i}}\mathbf{R}\dfrac{\partial {^{i}}\mathbf{A}^\top }{\partial \rho _e} , \end{aligned} \end{aligned}$$

(27)

under the assumption that $\partial {^{1}}\mathbf{R}/\partial \rho _e = \varvec{0}_{N+4}$ $\forall e$. The sensitivity analysis given by Eqs. (25) and (26) is detailed in Appendix 2.

It is noted that:

$$\begin{aligned} \dfrac{\partial {^{i+1}}\mathbf{A}}{\partial \rho _e} = \begin{bmatrix} \dfrac{\partial \mathbf{A}_s}{\partial \rho _e} &{} \varvec{0}_{N\times 4} \\ \varvec{0}_{4\times N} &{} \varvec{0}_{4} \end{bmatrix} , \end{aligned}$$

(28)

and thus it is constant $\forall i$ regardless the time-variation of amplitude and frequency content of the seismic ground motion, with

$$\begin{aligned} \dfrac{\partial \mathbf{A}_s}{\partial \rho _e} = \begin{bmatrix} \varvec{0}_{N_f} &{} \varvec{0}_{N_f} \\ \mathbf{M}^{-1}\dfrac{\partial \mathbf{M}}{\partial \rho _e} \mathbf{M}^{-1} \mathbf{K}- \mathbf{M}^{-1}\dfrac{\partial \mathbf{K}}{\partial \rho _e} &{} \mathbf{M}^{-1}\dfrac{\partial \mathbf{M}}{\partial \rho _e}\mathbf{M}^{-1}-\mathbf{M}^{-1}\dfrac{\partial \mathbf{C}}{\partial \rho _e} \end{bmatrix} . \end{aligned}$$

(29)

The derivatives $\partial \mathbf{M}/\partial \rho _e$, $\partial \mathbf{K}/\partial \rho _e$ and $\partial \mathbf{C}/\partial \rho _e$ are determined in agreement with the SIMP approach as detailed in the Appendix 3. Once Eqs. (19) and (25) are solved at each time instant, the gradient of the objective function given by Eq. (24) can be obtained by means of a numerical approximation of the involved integrals (i.e., via trapezoidal numerical integration).

4.3 Approximated gradient evaluation

Although Eq. (24) provides the exact way to estimate the gradient of the objective function, its use for topology optimization problems is not viable because the calculation of the explicit derivatives of the covariance matrix for the whole design domain through Eq. (25) requires huge elaboration time and memory consumption. The computational effort is further aggravated in case of non-stationary random vibrations because of the time discretization. In case of stationary random vibrations, Gomez and Spencer (2019) proposed a very attractive way to reduce the elaboration cost by means of the adjoint method. Unfortunately, that approach cannot be implemented directly herein to deal with non-stationary random vibrations. Because of this, an approximated construction of the objective function gradient is here proposed for the topology optimization of multi-story buildings subjected to earthquake.

First, it is opportune to recall that each finite element within the optimizable domain brings one design variable only according to the SIMP approach (i.e., the density value within the system topology), which implies that ${\partial \mathbf{K}}/{\partial \rho _e}={\partial \mathbf{K}_{\Upomega _e}}/{\partial \rho _e}$, i.e., ${\partial \mathbf{K}}/{\partial \rho _e}$ is null except for the degrees of freedom of the eth finite element within the optimizable domain $\Upomega$ for which the derivative is calculated. Therefore, the first addend in Eq. (24) simplifies as follows:

$$\begin{aligned} \mathbf{r}^\top \left( \dfrac{\partial \mathbf{K}}{\partial \rho _e} \otimes \mathbf{R}_{{\mathbf{u}}{\mathbf{u}}} \right) \mathbf{r}= \mathbf{r}_{\Upomega _e}^\top \left( \dfrac{\partial \mathbf{K}_{\Upomega _e}}{\partial \rho _e} \otimes \mathbf{R}_{{\mathbf{u}}{\mathbf{u}},{\Upomega _e}} \right) \mathbf{r}_{\Upomega _e}, \end{aligned}$$

(30)

where $\mathbf{r}_{\Upomega _e}$ and $\mathbf{R}_{{\mathbf{u}}{\mathbf{u}},{\Upomega _e}}$ are the incidence vector and the covariance matrix of the structural response that refers to the degrees of freedom of the eth finite element within the optimizable domain $\Upomega$.

A similar argument does not apply rigorously to the second addend in Eq. (24), but it is here conjectured that it holds true asymptotically for a sub-assembly including the eth finite element with respect to which the derivative is calculated and some finite elements within its neighborhood. Formally, it means that the following approximation is here accepted:

$$\begin{aligned} \dfrac{\partial \mathbf{R}_{{\mathbf{u}}{\mathbf{u}}} }{\partial \rho _e } \approx \bigwedge _{e \in \Uppi _e}{ \dfrac{\partial \mathbf{R}_{{\mathbf{u}}{\mathbf{u}},{\Upomega _e}}}{\partial \rho _e} }, \end{aligned}$$

(31)

where $\Uppi _e \subseteq \Upomega$ is the sub-assembly consisting of the eth finite element with respect to which the derivative is calculated and the finite elements within its neighborhood. Because of the time discretization, ${\partial \mathbf{R}_{{\mathbf{u}}{\mathbf{u}},{\Upomega _e}} }/{\partial \rho _e }$ is defined through a set of time instants as ${\partial {^{i+1}}\mathbf{R}_{{\mathbf{u}}{\mathbf{u}},{\Upomega _e}} }/{\partial \rho _e }$, which are obtained by solving Eq. (25) once it is specialized for the degrees of freedom involved in the sub-assembly $\Uppi _e$ only. Figure 2 illustrates some alternative definitions of the sub-assembly $\Uppi _e$.

Fig. 2

Some definitions of the sub-assembly for the approximated evaluation of the gradient: sub-assembly coinciding with the whole design domain (a), sub-assembly with 8 neighborhood finite elements (b), sub-assembly with 4 neighborhood finite elements (c), sub-assembly including the single finite element only (d)

The underlying idea in Eq. (31) rests on a simple mechanics-based reasoning: the elements of the covariance matrix of the structural response that refers to degrees of freedom far enough from those of the eth finite element should be almost insensitive to variations of the eth design variable value. The error associated with this approximation tends to zero as $\Uppi _e \rightarrow \Upomega$, but this causes an increment of the computational effort. Conversely, the error is maximum when $\Uppi _e$ degenerates into the eth finite element only, but this requires the lowest computational effort.

Taking into account Eqs. (30) and (31), the least computationally demanding approximation of Eq. (24) is thus based on a sub-assembly consisting of a single finite element, and it reads:

$$\begin{aligned} \begin{aligned} \dfrac{\partial f }{\partial \rho _e} \approx&\int _0^{T}{ \ \mathbf{r}_{\Upomega _e}^\top \left( \dfrac{\partial \mathbf{K}_{\Upomega _e}}{\partial \rho _e} \otimes \mathbf{R}_{{\mathbf{u}}{\mathbf{u}},{\Upomega _e}} \right) \mathbf{r}_{\Upomega _e} \text{ d }t }\\&+ \int _0^{T}{ \mathbf{r}_{\Upomega _e}^\top \left( \mathbf{K}_e \otimes \dfrac{\partial \mathbf{R}_{{\mathbf{u}}{\mathbf{u}},{\Upomega _e}}}{\partial \rho _e} \right) \mathbf{r}_{\Upomega _e} \text{ d }t } \end{aligned} , \end{aligned}$$

(32)

where $\mathbf{K}_e$ includes the stiffness matrix of the eth element within the optimizable domain $\mathbf{K}_{\Upomega _e}$ and the pertinent contribution of the stiffness attributable to the (non-optimizable) secondary frame system $\mathbf{K}_{\Upgamma _e}$.

It is evident that computational time and memory consumption for calculating the approximated gradient according to Eq. (32) are much less than the ones required for the exact estimate in Eq. (24). This is mostly due to the huge reduction of the size of the Lyapunov equation given by Eq. (25) required to compute the second addend of the approximation in Eq. (32). Because of the time discretization, the approximated gradient of the objective function given by Eq. (32) is once again estimated through a numerical approximation of the involved integrals (i.e., via trapezoidal numerical integration). To further reduce the overall elaboration time, efficient algorithms for solving the involved Lyapunov equations should be preferred in place of the standard ones (e.g., Bartel-Stewards algorithm), see about this matter the recent survey by (Benner and Saak 2013).

If needed, the accuracy of the approximated gradient evaluation given by Eq. (32) can be improved by expanding the size of the sub-assembly $\Uppi _e$ in such a way to include other finite elements within the neighborhood apart from the eth finite element for which the derivative is calculated (this will reduce the error attributable to the second term in Eq. (32), whereas the first term remains unchanged regardless the number of elements included in the sub-assembly).

It is highlighted that the proposed approach for topology optimization deals with damage control under frequent and low intensity earthquakes, for which the structure must be designed to remain elastic and undergo slightly damage only (Elettore et al. 2020). Conversely, the design for collapse prevention in case of rare and high intensity seismic events requires to account for nonlinear behavior and extensive damage (Elettore et al. 2020). In such a case, the proposed approach can be still implemented, provided that the nonlinear behavior (e.g., a hysteretic response, possibly together with stiffness and/or strength degradation) is taken into account through a suitable linearization procedure (Roberts and Spanos 2003).

For the sake of completeness, it is finally pointed out that the conversion of the optimal topological solution into a manufacturable structural layout consisting of an assembly of classical members (e.g., beam-type elements) requires the application of image processing techniques. To this end, several post-processing procedures have been proposed to translate the element densities of the optimized topology (i.e., low-level abstraction) into a skeleton structure (i.e., high-level abstraction) (e.g., Chou and Lin 2010; Gamache et al. 2018; Kazakis et al. 2017; Nana et al. 2017). Once topology results are interpreted coherently into parametric computer-aided design models, the resulting layout is ready for size optimization and/or final manufacturing. In the case of concrete structures, to which the following numerical applications refer to, the final realization can be carried out using either cast-in-place or precast elements, or both. A viable process for concrete bracing frames has been presented by Qiao et al. (2017).

5 Numerical investigations

5.1 Validation

The proposed gradient approximation based on Eq. (32) needs a preliminary validation in such a way to assess its accuracy and computational efficiency. To this end, it is considered the topology optimization of a simple ideal structural system, so as to facilitate the reproducibility of the following results. The optimizable domain of such system has a square shape whose side length is equal to 5 m, and it is discretized with Q4 finite elements whose thickness is equal to 0.1 m. The structure is idealized with fixed supports at the base. Elastic modulus and density are taken equal to 21 GPa and 2400 kg/m$^3$, respectively. The viscous damping matrix follows the Rayleigh’s model, where the viscous damping ratio is assumed equal to 5% for the first two modes of vibration. The mass is calculated considering dead and live loads equal to 7.0 kN/m$^2$ and 2.0 kN/m$^2$, respectively, and is concentrated in both corners of the top. The topology design starts with a uniform density distribution within the optimizable domain equal to 0.30, and symmetry constraints about the centerline are enforced. A linear density-based regularization filter with radius equal to 0.3 m is employed. The approximated gradient evaluation is performed according to the proposed procedure by considering a single finite element only as well as 8 or 4 additional neighborhood finite elements as shown in Fig. 2b–d.

The validation is conducted initially for the simpler case of stationary base excitation with $\ddot{u}_g^{\text {max}}=0.3\text {g}$, $\xi _k=\xi _p=0.60$, $\omega _k=15.00 \text{ rad/s }$, $\omega _p=1.50 \text{ rad/s }$. The reference optimized solution is obtained using the adjoint method as described by Gomez and Spencer (2019), which is a convenient, yet exact, way for topology optimization in case of stationary dynamic loads. The optimum design is first performed for a structured mesh consisting of squares with size $0.50 \text{ m } \times 0.50 \text{ m }$, corresponding to a total of 100 Q4 finite elements within the optimizable domain. Such coarse mesh is useful to facilitate the element-by-element check of the approximation error corresponding to the proposed approach as compared to an exact computation of the problem gradient. Figure 3 shows an excellent agreement between the optimized topology calculated by means of the adjoint method according to Gomez and Spencer (2019) (Fig. 3a) and the results obtained by means of the proposed approach based on the approximated gradient evaluation (Fig. 3b–d). Notably, the optimized topologies calculated through the proposed approach are basically equal to each other (Fig. 3b–d) regardless the number of finite elements in the sub-assembly. The optimal values of the design variables are also made explicit in Fig. 3 for better appreciating the differences among the solution strategies. This series of plots confirms that only a few design variables calculated by means of the proposed approach differ from the corresponding reference solution obtained according to (Gomez and Spencer 2019). As expected, the larger is the number of finite elements included in the sub-assembly employed for the approximated evaluation of the objective function gradient, the smaller is the difference (i.e., the optimal solution obtained through the approximated gradient evaluation tends toward the exact one as the number of finite elements involved in the sub-assembly grows up). Anyway, the difference in terms of design variable values is negligible from a practical standpoint for all sizes of the sub-assembly.

Fig. 3

Validation for an optimizable domain consisting of Q4 finite elements $0.50 \text{ m } \times 0.50 \text{ m }$ under stationary base motion: optimized topologies (first row) and density values (second row) calculated via adjoint method (Gomez and Spencer 2019) (a) and the proposed approximation considering sub-assemblies with 8 or 4 neighborhood finite elements or a single finite element only (b–d). The underlined numbers identify the density values obtained by means of the proposed approximation that differ from the solution calculated according to (Gomez and Spencer 2019)

Figure 4 is helpful for confirming the agreement among the solution strategies in term of objective function sensitivity as calculated at the last iteration of the optimization procedure. In this regard, it is recalled that observed differences are attributable to the evaluation of the derivatives of the covariance matrix of the structural response according to Eq. (31). These plots confirm that the proposed approximation leads to an objective function gradient that very well agrees with the reference solution carried out following Gomez and Spencer (2019) regardless the assumed size of the sub-assembly. As a matter of fact, the differences can be inferred only by looking at the colormap interval, and do not produce practical effects on the final density values as already pointed out in Fig. 3.

Fig. 4

Validation for an optimizable domain consisting of Q4 finite elements $0.50 \text{ m } \times 0.50 \text{ m }$ under stationary base motion: sensitivity of the objective function calculated via adjoint method (Gomez and Spencer 2019) (a) and the proposed approximation considering sub-assemblies with 8 or 4 neighborhood finite elements or a single finite element only (b–d)

In order to evaluate the influence of the mesh refinement, a structured mesh consisting of squares with size $0.25 \text{ m } \times 0.25 \text{ m }$ is also considered in Fig. 5, which corresponds to a total of 400 Q4 finite elements within the optimizable domain. These results confirm that, even for a finer mesh, there exist an almost perfect agreement between the results carried out by means of the proposed procedure and the reference solution obtained through the adjoint method-based approach developed by Gomez and Spencer (2019).

Fig. 5

Validation for an optimizable domain consisting of Q4 finite elements $0.25 \text{ m } \times 0.25 \text{ m }$ under stationary base motion: optimized topology (first row) and sensitivity of the objective function (second row) calculated via adjoint method (Gomez and Spencer 2019) (a) and the proposed approximation considering sub-assemblies with 8 or 4 neighborhood finite elements or a single finite element only (b–d)

Altogether, Figs. 3, 4 and 5 demonstrate that the least computationally demanding approximation of the objective function gradient given by Eq. (32) is accurate enough since it replicates almost exactly density and sensitivity values of the reference results calculated according to Gomez and Spencer (2019), whereas increasing the size of the sub-assembly (i.e., including some finite elements within the neighborhood) leads to slight numerical improvements with no practical effects on the final topology. The implemented density-based filtering technique also plays an important role in this regard. In fact, it smears and alleviates the approximation error in the estimated density values by computing the density of each finite element as the weighted average of the corresponding values within its neighborhood before the finite element-based solver is invoked, and the sensitivities are modified next in a consistent way (Bendsøe and Sigmund 1995). In terms of computational efficiency, the proposed approach also well compares with the exact methodology proposed by Gomez and Spencer (2019). The topology optimization based on the approximated gradient evaluation according to Eq. (32) (i.e., sub-assembly including a single finite element) was found $+6$% slower on average. This slight increment of the elaboration time is considered acceptable in light of the fact that the proposed approach can be applied to deal with the more general load case of fully non-stationary stochastic base motion without special restrictions about the time-variation of amplitude and frequency content of the excitation.

The general validity of the proposed framework has been confirmed numerically in case of non-stationary base motion as well. Since the approach by Gomez and Spencer (2019) cannot be applied in this scenario, the reference solution in such a case is obtained by direct differentiation over the whole optimizable domain according to Eq. (24). Note that the elaboration cost is even more important under non-stationary excitation because of the dramatic increment of the computational efforts due to the discretization of the time window. In this sense, Fig. 6 highlights the elaboration times for a non-stationary base motion. Considering the first 10 iterations of the optimization procedure, the use of the approximation given by Eq. (32) (that refers to a single finite element in the sub-assembly) resulted about 100 and 1000 times faster than the exact formulation in Eq. (24) (which also corresponds to a sub-assembly coinciding with the whole design domain) for meshes with 100 or 400 finite elements, respectively, without producing any evident effect on the final topology.

Fig. 6

Elaboration times for an optimizable domain consisting of Q4 finite elements $0.50 \text{ m } \times 0.50 \text{ m }$ and $0.25 \text{ m } \times 0.25 \text{ m }$ under non-stationary base motion (reported values are normalized with respect to the elaboration time required for the sub-assembly including the single finite element only)

It is important to remark that the feasibility of the present approximated construction of the gradient has been thoroughly verified for the class of topology optimization problems of interest in the present study. A preliminary investigation is needed to verify its feasibility for applications other than the ones addressed in the present work.

5.2 Applications

After having demonstrated the correctness of the proposed approach, it is now applied to determine the LRS of multi-story buildings subjected to non-stationary seismic ground motion. So doing, the proposed approximation of the gradient is calculated for a sub-assembly including a single finite element only according to Eq. (32), since the numerical validation has confirmed that it practically reproduces the exact optimized topology. Moreover, the associated minimum computational effort makes this level of approximation especially attractive for the optimization of large structures. The main aim of the following applications is to assess the role of the seismic ground motion characteristics on topology optimization.

The non-stationarity of the frequency content of the stochastic seismic ground motion is taken into account by assuming that the dominant circular frequency $\omega _k$ is time-dependent. To this end, the following expressions due to Fan and Ahmadi (1990) are considered:

$$\begin{aligned}&\omega _k = 9.425+59.722(\exp {(-0.0625t)}-\exp {(-0.15t)}), \end{aligned}$$

(33a)

$$\begin{aligned}&\omega _k = 3.456+2.827\sin {\left( 0.17(t-2)\right) }, \end{aligned}$$

(33b)

which are representative of a firm soil and a soft soil, respectively. Particularly, these expressions were calibrated by Fan and Ahmadi (1990) for El Centro 1940 earthquake and Mexico City 1985 earthquake, respectively. The other filter parameters are assumed as constant values according to Yeh and Wen (1990). Specifically, $\xi _k=0.65$, $\omega _p=2 \text { rad/s}$ and $\xi _p=0.6$ for firm soil and $\xi _k=0.10$, $\omega _p=2.3 \text { rad/s}$ and $\xi _p=0.1$ for soft soil. In order to calculate the constant value of $S_0$, it is assumed $\omega _k=19 \text { rad/s}$ and $\omega _k=4.2 \text { rad/s}$ for firm and soft soil, respectively, as proposed by Yeh and Wen (1990). These constant values of the filter parameters will be also considered when performing the topology optimization under stationary stochastic seismic ground motion. It is pointed out that all these constant filter parameters were derived by Yeh and Wen (1990) for El Centro and Mexico City earthquake, respectively.

The non-stationary in amplitude of the stochastic seismic ground motion is simulated by means of the modulation function due to Jennings et al. (1969), which reads:

$$\begin{aligned} \varphi = \left\{ \begin{array}{lll} \left( t/t_a\right) ^2 &{}t<t_a\\ 1 &{}t_a \le t \le t_b\\ \exp {(-a(t-t_b))} &{}t>t_b\\ \end{array} \right. , \end{aligned}$$

(34)

where $t_a$ and $t_b$ denotes starting and ending time instant of the strong-motion phase while a is the shape parameter of the decaying phase. Taking into account the study by Fan and Ahmadi (1990), it is here assumed $T=20 \text { s}$ with $t_b-t_a=5\text { s}$ for firm soil, whereas $T=50 \text { s}$ with $t_b-t_a= 30\text { s}$ is considered for soft soil. For all soil types, it is adopted $t_a=1 \text { s}$ and $a=0.5$. The PGA value is $\ddot{u}_g^{\text {max}}=0.2\text {g}$.

Multi-story buildings are hereafter analyzed by splitting the envelope into four panels such that each façade can be considered as an independent planar continuous optimizable domain, which is discretized by using Q4 finite elements with uniform thickness equal to 0.1 m and mesh size $0.1 \text { m} \times 0.1 \text { m}$. The structure is idealized with fixed supports at the base, and symmetry constraints with respect to the vertical centerline are enforced in the optimization procedure. Two lateral columns bounding the continuous domain are included in the non-optimizable domain. Normal-weight concrete is adopted for both domains, with elastic modulus and density equal to 21 GPa and 2400 kg/m$^3$, respectively. The optimization is performed using a linear density-based regularization filter with radius equal to 0.2 m, and initial volume fraction over the optimizable domain equal to 0.25. An additional floor mass equal to 4000 kg is lumped at the master nodes of each story level. The Rayleigh’s model of the viscous damping matrix is adopted with a viscous damping ratio equal to 5%.

The first case-study deals with an ideal 3-story building, a well-known benchmark for topology optimization algorithms (Chun et al. 2016; Bobby et al. 2017). Height and width of the façade are equal to 15 m and 5 m, respectively, and floor slabs are placed every 5 m. The optimizable domain is discretized with a structure mesh using 150 and 50 Q4 finite elements along the vertical and the horizontal direction, respectively, for a total of 7500 finite elements. The resulting mesh size is similar to that adopted in previous studies for this benchmark (e.g., Chun et al. 2016). Lateral columns of the secondary system are modeled using beam finite elements whose length is equal to the side of the single Q4 finite element (i.e., 0.1 m), thereby leading to a total of 150 beam elements into each column. Two different cross-sections of the lateral columns are considered, namely $0.5\text { m} \times 0.5\text { m}$ and $0.6\text { m} \times 0.6\text { m}$. Two master nodes with lumped mass are considered at each floor level. General layout and discretization of this 3-story building are shown in Fig. 7a. The optimized topologies in case of non-stationary stochastic seismic ground motion for different soil types and cross-section size of the lateral columns are given in Fig. 7b–e.

Fig. 7

Topology optimization of 3-story building under non-stationary stochastic seismic ground motion: layout and discretization of the building with master nodes position marked by red dots (a), optimum design for firm soil condition with cross-sections of the lateral columns equal to $0.5\text { m} \times 0.5\text { m}$ (b), optimum design for firm soil condition with cross-sections of the lateral columns equal to $0.6\text { m} \times 0.6\text { m}$ (c), optimum design for soft soil condition with cross-sections of the lateral columns equal to $0.5\text { m} \times 0.5\text { m}$ (d), optimum design for soft soil condition with cross-sections of the lateral columns equal to $0.6\text { m} \times 0.6\text { m}$ (e)

Final results indicate that the seismic ground motion parameters can have a significant impact on the geometry of the final layouts. This manifests in terms of brace-to-brace and brace-to-column working points as well as on the bracings shape. The different arrangement of the structural material is especially evident in the lower parts of the topologies, whereas it maintains a 45$^\circ$ angle in the upper levels. This suggests that the most effective way to meet the stiffness demand in seismically excited structures is to strengthen the lower levels of the structural system first by adjusting the bracing geometry. Figure 7b and d highlights that both bracing-column points and intersection points of the bracings at the lower floor slightly move up vertically when passing from firm to soft soil conditions even if the same column cross-section is considered (i.e., $0.5\text { m} \times 0.5\text { m}$). On the other hand, the increment of the column cross-section size leads to the overall uplift of the first bracing module, as shown in Fig. 7c and e. Both anchoring nodes at the column base and bracing-column points at the first level moves upward, thus resulting in steeper geometry for the two upper braces of the module. The diagonal arrangement results in a stiffer configuration because of the larger inertial loads applied to the master nodes along the column height due to the increment of the mass associated with the larger column cross-section size. Meanwhile, the larger contribution of the columns of the LRS in reducing the stiffness demand of the optimizable domain alleviates the material concentration at the edges of the base module. Additionally, the increment of the column size under soft soil condition leads to a more traditional X-bracing topology, for which stiffening members are no longer required. Corresponding convergence histories for volume and objective function are plotted into Fig. 8.

Fig. 8

Topology optimization of 3-story building under non-stationary stochastic seismic ground motion: convergence history of constraint (a) and objective function (b) for different soil types and cross-sections of the lateral columns

The curves in Fig. 8 show that convergence has been achieved successfully. Particularly, a structural material distribution satisfying the volume constraint is rapidly identified, whereas further iterations are required to attain a stable, minimum objective function value. The optimized topologies exhibit improved dynamic performance though only about 25% of the initial volume is used. For an equal volume cut, a reduction in the objective function value equal to 55% and 80% is observed for firm and soft soil conditions, respectively, with respect to the initial design. The same case-study is also examined under a stationary stochastic seismic ground motion to infer whether this leads to different results. The optimized topologies under stationary seismic ground motion are reported in Fig. 9.

Fig. 9

Topology optimization of 3-story building under stationary stochastic seismic ground motion: layout and discretization of the building with master nodes position marked by red dots (a), optimum design for firm soil condition with cross-sections of the lateral columns equal to $0.5\text { m} \times 0.5\text { m}$ (b), optimum design for firm soil condition with cross-sections of the lateral columns equal to $0.6\text { m} \times 0.6\text { m}$ (c), optimum design for soft soil condition with cross-sections of the lateral columns equal to $0.5\text { m} \times 0.5\text { m}$ (d), optimum design for soft soil condition with cross-sections of the lateral columns equal to $0.6\text { m} \times 0.6\text { m}$ (e)

The comparison of Figs. 7 and 9 highlights that final topologies are similar, except in the case of soft soil condition and cross-sections of the lateral columns equal to $0.5\text { m} \times 0.5\text { m}$. As it is evident by comparing Figs. 7d and 9d, the additional stiffening branches at the first module that were found in case of non-stationary stochastic seismic ground motion disappear in case of stationary excitation.

The last case-study deals with a 5-story building. Height and width of the façade are equal to 25 m and 5 m, respectively, and floor slabs are placed every 5 m. The optimizable domain is discretized using 250 and 50 Q4 finite elements along the vertical and the horizontal direction, respectively, for a total of 12,500 finite elements. Lateral columns of the secondary system have a cross-section size equal to $0.7\text { m} \times 0.7\text { m}$. Two master nodes with lumped mass are considered at each floor location. The building is analyzed under firm soil condition, and it was initially designed in such a way that its initial fundamental frequency (i.e., fundamental frequency corresponding to the initial volume fraction over the optimizable domain equal to 0.25) was as close as possible to the constant value of the dominant frequency of the soil (i.e., $\omega _k=19 \text { rad/s}$). Specifically, the initial fundamental frequency is $\omega _1=19.5 \text { rad/s}$. General layout and discretization of this 5-story building are shown in Fig. 10a. The optimized topologies in case of non-stationary and stationary stochastic seismic ground motion for firm soil condition are given in Fig. 10b–c, respectively.

Fig. 10

Topology optimization of 5-story building for firm soil condition: layout and discretization of the building with master nodes position marked by red dots (a), optimum design for non-stationary stochastic seismic ground motion (b), optimum design for stationary stochastic seismic ground motion (c), optimum design for equivalent static load (d)

At the end of the optimization procedure, it is found $\omega _1=26.6 \text { rad/s}$ and $\omega _1=22.6 \text { rad/s}$ for a non-stationary and a stationary stochastic seismic ground motion. Therefore, the optimization procedure distributes the structural material in order to increase the fundamental natural frequency of the building with respect to that of the soil beneath. Additionally, the optimum design under a non-stationary stochastic seismic ground motion leads to a LRS stiffer than that obtained in case of stationary excitation (i.e., the stiffness demand is higher for non-stationary base excitation rather than for stationary ground motion). The topology optimization procedure fulfills the volume constraint in both seismic scenarios by strengthening the lower levels while producing almost the same layout in the upper part. However, the optimized topologies in Fig. 10b, c are qualitatively different in terms of bracing arrangement at the first three modules. For non-stationary stochastic seismic ground motion, the stiffness demand is mainly fulfilled by modifying the geometry of the braces through the generation of additional strengthening elements. Conversely, for a stationary stochastic seismic ground motion, the lateral bracing system minimizes the compliance by placing more structural material at the second module and adding stiffening arms at the third module. The topology optimization for this last case-study is also performed by means of an equivalent static load in order to better understand the role of the loading conditions on final results. Specifically, following a well-established procedure in seismic analysis of structures, it is assumed an inverted triangle-type distribution of forces proportional to the floor masses. The spectral acceleration of the El Centro earthquake is determined after PGA scaling, in agreement with the data employed for modeling the stochastic seismic ground motion. Figure 10d highlights that the use of an equivalent static load leads to brace-to-brace and column-to-brace nodes having different arrangements as compared to a dynamic excitation (as shown in Fig. 10b and c for non-stationary and stationary seismic ground motion, respectively). This is especially evident in the lower half of the building. Overall, these results indicates that the optimized stiffness demand can strongly depend on the representation of the seismic loading.

6 Conclusions

The present work has addressed the topology optimization of the lateral resisting system in multi-story buildings under earthquake. To this end, a random vibration-based approach has been adopted to cope with the inherent randomness of the earthquakes, where the seismic ground motion is simulated as filtered white Gaussian noise with time-varying amplitude and frequency content (i.e., fully non-stationary seismic ground motion). Because of the large computational effort associated with the resolution of this optimization problem, an approximated construction of the gradient has been proposed in which explicit, exact derivatives with respect to the design variables are computed analytically through direct differentiation for a sub-assembly of finite elements resulting from the discretization of the optimizable domain. Within the framework of the proposed approximation, the sub-assembly consists of the finite element with respect to which the derivative is calculated and some finite elements within its neighborhood. The minimum sub-assembly consisting of a single element was found to be good enough for obtaining a satisfactory accuracy at lowest computational effort.

Numerical investigations have substantiated accuracy, efficiency and objectivity of the proposed approach based on the approximated evaluation of the gradient. The results have also demonstrated that the characteristics of the input seismic ground motion can drastically affect the final topologies. Particularly, modeling the earthquake in a more effective way as a fully non-stationary stochastic process can lead to optimized layouts very different from that carried out under the simplified assumption of stationary excitation. It is thus concluded that a reliable topology optimization under non-stationary dynamic loads should take into proper account the time-variation of amplitude and frequency content of the excitation.

Appendices

Appendix 1

This Appendix details the time discretization of the non-stationary Lyapunov equation given by Eqs. (19)–(21). For each time instant $i+1$ and i, the time discretization of Eq. (17) implies that:

$$\begin{aligned}&{^{i+1}}\mathbf{A}{^{i+1}}\mathbf{R}+ {^{i+1}}\mathbf{R}\mathbf{A}^\top + {^{i+1}}\mathbf{B}= {^{i+1}}\dot{\mathbf{R}}, \end{aligned}$$

(41)

$$\begin{aligned}&{^{i}}\mathbf{A}{^{i}}\mathbf{R}+ {^{i}}\mathbf{R}{^{i}}\mathbf{A}^\top + {^{i}}\mathbf{B}= {^{i}}\dot{\mathbf{R}}. \end{aligned}$$

(42)

Once that Eqs. (41)–(42) are added and the result is multiplied by $\frac{1}{2}\Updelta t$, by taking into account Eq. (18) is obtained:

$$\begin{aligned} \begin{aligned}&{^{i+1}}\mathbf{R}- {^{i}}\mathbf{R}= \\&\dfrac{1}{2}\Updelta t \left( {^{i+1}}\mathbf{A}{^{i+1}}\mathbf{R}+ {^{i+1}}\mathbf{R}{^{i+1}}\mathbf{A}^\top + {^{i+1}}\mathbf{B}+ {^{i}}\mathbf{A}{^{i}}\mathbf{R}+ {^{i}}\mathbf{R}{^{i}}\mathbf{A}^\top + {^{i}}\mathbf{B}\right) . \end{aligned} \end{aligned}$$

(43)

Equation (43) is first rewritten as follows:

$$\begin{aligned} \begin{aligned}&\left( \dfrac{1}{2}\Updelta t {^{i+1}}\mathbf{A}{^{i+1}}\mathbf{R}- {^{i+1}}\mathbf{R}\right) + \dfrac{1}{2}\Updelta t {^{i+1}}\mathbf{R}{^{i+1}}\mathbf{A}^\top + \\&\left( {^{i}}\mathbf{R}+ \dfrac{1}{2}\Updelta t {^{i}}\mathbf{A}{^{i}}\mathbf{R}+ \dfrac{1}{2}\Updelta t {^{i}}\mathbf{R}{^{i}}\mathbf{A}^\top + \dfrac{1}{2}\Updelta t {^{i+1}}\mathbf{B}+ \dfrac{1}{2}\Updelta t {^{i}}\mathbf{B}\right) = \varvec{0}_{N+4}, \end{aligned} \end{aligned}$$

(44)

from which:

$$\begin{aligned} \begin{aligned}&\left( - \dfrac{1}{2}{^{i+1}}\mathbf{R}+ \dfrac{1}{2}\Updelta t {^{i+1}}\mathbf{A}{^{i+1}}\mathbf{R}\right) + \left( - \dfrac{1}{2}{^{i+1}}\mathbf{R}\right. \\&\left. \quad + \dfrac{1}{2}\Updelta t {^{i+1}}\mathbf{R}{^{i+1}}\mathbf{A}^\top \right) \\&\left( {^{i}}\mathbf{R}+ \dfrac{1}{2}\Updelta t {^{i}}\mathbf{A}{^{i}}\mathbf{R}+ \dfrac{1}{2}\Updelta t {^{i}}\mathbf{R}{^{i}}\mathbf{A}^\top + \dfrac{1}{2}\Updelta t {^{i+1}}\mathbf{B}+ \dfrac{1}{2}\Updelta t {^{i}}\mathbf{B}\right) \\&\quad = \varvec{0}_{N+4}. \end{aligned} \end{aligned}$$

(45)

Equation (45) coincides with Eq. (19) by introducing the quantities in Eqs. (20)–(21).

Appendix 2

More details about the sensitivity analysis given by Eqs. (25)–(26) are provided in this Appendix. Initially, Eq. (19) is rewritten as follows:

$$\begin{aligned}&-\dfrac{1}{2}{^{i+1}}\mathbf{R}+ \dfrac{1}{2}\Updelta t {^{i+1}}\mathbf{A}{^{i+1}}\mathbf{R}- \dfrac{1}{2}{^{i+1}}\mathbf{R}\nonumber \\& + \dfrac{1}{2}\Updelta t{^{i+1}}\mathbf{R}{^{i+1}}\mathbf{A}^\top + {^{i+1}}\bar{\mathbf{B}}=\varvec{0}_{N+4}. \end{aligned}$$

(46)

By applying the derivative operator $\partial /\partial \rho _e$ to Eq. (46) is obtained:

$$\begin{aligned} \begin{aligned}&-\dfrac{1}{2}\dfrac{\partial {^{i+1}}\mathbf{R}}{\partial \rho _e} + \dfrac{1}{2}\Updelta t \dfrac{\partial {^{i+1}}\mathbf{A}}{\partial \rho _e} {^{i+1}}\mathbf{R}+ \dfrac{1}{2}\Updelta t {^{i+1}}\mathbf{A}\dfrac{\partial {^{i+1}}\mathbf{R}}{\partial \rho _e} - \\&\dfrac{1}{2}\dfrac{\partial {^{i+1}}\mathbf{R}}{\partial \rho _e} + \dfrac{1}{2}\Updelta t \dfrac{\partial {^{i+1}}\mathbf{R}}{\partial \rho _e}{^{i+1}}\mathbf{A}^\top + \dfrac{1}{2}\Updelta t{^{i+1}}\mathbf{R}\dfrac{\partial {^{i+1}}\mathbf{A}^\top }{\partial \rho _e} \\&+ \frac{\partial {^{i+1}}\bar{\mathbf{B}}}{\partial \rho _e}=\varvec{0}_{N+4}. \end{aligned} \end{aligned}$$

(47)

Equation (47) is finally rewritten as follows:

$$\begin{aligned} \begin{aligned}&\left( -\dfrac{1}{2}\mathbf{I}_{N+4} + \dfrac{1}{2}\Updelta t {^{i+1}}\mathbf{A}\right) \dfrac{\partial {^{i+1}}\mathbf{R}}{\partial \rho _e} \\&\quad + \dfrac{\partial {^{i+1}}\mathbf{R}}{\partial \rho _e}\left( -\dfrac{1}{2}\mathbf{I}_{N+4} + \dfrac{1}{2}\Updelta t {^{i+1}}\mathbf{A}^\top \right) \\&\quad +\left( \dfrac{1}{2}\Updelta t \dfrac{\partial {^{i+1}}\mathbf{A}}{\partial \rho _e} {^{i+1}}\mathbf{R}+ \dfrac{1}{2}\Updelta t {^{i+1}}\mathbf{R}\dfrac{\partial {^{i+1}}\mathbf{A}^\top }{\partial \rho _e} + \frac{\partial {^{i+1}}\bar{\mathbf{B}}}{\partial \rho _e} \right) =\varvec{0}_{N+4}. \end{aligned} \end{aligned}$$

(48)

Equation (48) coincides with Eq. (25) by introducing the quantities in Eqs. (20) and (26).

Appendix 3

This Appendix details the calculation of $\partial \mathbf{M}/\partial \rho _e$, $\partial \mathbf{K}/\partial \rho _e$ and $\partial \mathbf{C}/\partial \rho _e$ in Eq. (29). To this end, by virtue of the SIMP approach, it is obtained that:

$$\begin{aligned}&\dfrac{\partial \mathbf{M}}{\partial \rho _e} = \bigwedge _{e \in \Upomega }{\dfrac{\mathbf{M}_{\Upomega _e}}{\partial \rho _e}}= \bigwedge _{e \in \Upomega }{q\rho _e^{q-1}\mathbf{M}_{\Upomega _e}^0}. \end{aligned}$$

(49)

$$\begin{aligned}&\dfrac{\partial \mathbf{K}}{\partial \rho _e} = \bigwedge _{e \in \Upomega }{\dfrac{\mathbf{K}_{\Upomega _e}}{\partial \rho _e}}= \bigwedge _{e \in \Upomega }{p\rho _e^{p-1}\mathbf{K}_{\Upomega _e}^0}, \end{aligned}$$

(50)

Because of the Rayleigh’s model of the linear viscous damping matrix $\mathbf{C}$ in Eq. (2), it is also obtained:

$$\begin{aligned} \dfrac{\partial \mathbf{C}}{\partial \rho _e} = \dfrac{\partial a_0}{\partial \rho _e} \mathbf{M}+ a_0 \dfrac{\partial \mathbf{M}}{\partial \rho _e} + \dfrac{\partial a_1}{\partial \rho _e} \mathbf{K}+ a_1 \dfrac{\partial \mathbf{K}}{\partial \rho _e} , \end{aligned}$$

(51)

where

$$\begin{aligned}&\dfrac{\partial a_0}{\partial \rho _e} = \dfrac{\partial a_0}{\partial \omega _1}\dfrac{\partial \omega _1}{\partial \rho _e} + \dfrac{\partial a_0}{\partial \omega _2}\dfrac{\partial \omega _2}{\partial \rho _e}, \end{aligned}$$

(52)

$$\begin{aligned}&\dfrac{\partial a_1}{\partial \rho _e} = \dfrac{\partial a_1}{\partial \omega _1}\dfrac{\partial \omega _1}{\partial \rho _e} + \dfrac{\partial a_1}{\partial \omega _2}\dfrac{\partial \omega _2}{\partial \rho _e}. \end{aligned}$$

(53)

It is pointed out that the search for $\partial \mathbf{C}/\partial \rho _e$ as performed in the present study differs from that presented by Gomez and Spencer (2019). Both studies implement the Rayleigh’s model but they differ in the way by which $a_0$ and $a_1$ are considered. In fact, Gomez and Spencer (2019) assume that $a_0$ and $a_1$ are constants that do not depend on the design variables. According to the Rayleigh’s model, this implies that the viscous damping ratio $\xi _s$ is not an invariant system design property. However, given the building material and the construction system, the viscous damping ratio value $\xi _s$ is usually an assigned constant design data at least for the lowest modes of vibration. That’s why the two parameters $a_0$ and $a_1$ are here defined in such a way that $\xi _s$ is constant for the first two modes of vibration as shown in Eqs. (3a)–(3b). This, in turn, implies that $a_0$ and $a_1$ must depend on the design variables. In detail, for the calculation of $\partial a_0/\partial \rho _e$ and $\partial a_1/\partial \rho _e$, it is first observed that:

$$\begin{aligned}&\dfrac{\partial a_0}{\partial \omega _1} = \dfrac{2\xi _s\omega _2^2}{\left( \omega _1+\omega _2\right) ^2}, \ \ \dfrac{\partial a_0}{\partial \omega _2} = \dfrac{2\xi _s\omega _1^2}{\left( \omega _1+\omega _2\right) ^2}, \end{aligned}$$

(54)

$$\begin{aligned}&\dfrac{\partial a_1}{\partial \omega _1} = \dfrac{-2\xi _s}{\left( \omega _1+\omega _2\right) ^2}, \ \ \dfrac{\partial a_1}{\partial \omega _2} = \dfrac{-2\xi _s}{\left( \omega _1+\omega _2\right) ^2}, \end{aligned}$$

(55)

whereas

$$\begin{aligned} \dfrac{\partial \omega _1}{\partial \rho _e} = \dfrac{1}{2\omega _1}\dfrac{\partial \lambda _1}{\partial \rho _e}, \ \ \dfrac{\partial \omega _2}{\partial \rho _e} = \dfrac{1}{2\omega _2}\dfrac{\partial \lambda _2}{\partial \rho _e} \end{aligned}$$

(56)

are the derivatives of the first two circular natural frequencies of the structure. The terms $\lambda _1$ and $\lambda _2$ are the first two eigenvalues of the following eigenproblem:

$$\begin{aligned} \varvec{\phi }_j^\top \mathbf{K}\phi _k = \lambda _k\varvec{\phi }_k^\top \mathbf{M}\varvec{\phi }_k, \end{aligned}$$

(57)

where $\varvec{\phi }_k$ are the eigenvectors that constitute the mode shapes matrix of the structural system. The derivatives of $\partial \lambda _1 / \partial \rho _e$ and $\partial \lambda _2 / \partial \rho _e$ are evaluated as follows (Fox and Kapoor 1968):

$$\begin{aligned} \dfrac{\partial \lambda _1}{\partial \rho _e}=\varvec{\phi }_1^\top \left[ \dfrac{\partial \mathbf{K}}{\partial \rho _e}-\lambda _1\dfrac{\partial \mathbf{M}}{\partial \rho _e} \right] \varvec{\phi }_1, \ \ \dfrac{\partial \lambda _2}{\partial \rho _e}=\varvec{\phi }_2^\top \left[ \dfrac{\partial \mathbf{K}}{\partial \rho _e}-\lambda _2\dfrac{\partial \mathbf{M}}{\partial \rho _e} \right] \varvec{\phi }_2, \end{aligned}$$

(58)

under the assumption of distinct real eigenvalues and mass orthonormalization condition for the mode shapes matrix.