Controlling interstory drift ratio profiles via topology optimization strategies

Abstract

An approach to control the profiles of interstory drift ratios along the height of building structures via topology optimization is proposed herein. The theoretical foundation of the proposed approach involves solving a min–max optimization problem to suppress the maximum interstory drift ratio among all stories. Two formulations are suggested: one inherits the bound formulation and the other utilizes a p-norm function to aggregate all individual interstory drift ratios. The proposed methodology can shape the interstory drift ratio profiles into inverted triangular or quadratic patterns because it realizes profile control using a group of shape weight coefficients. The proposed formulations are validated via a series of numerical examples. The disparity between the two formulations is clear. The optimization results show the optimal structural features for controlling the interstory drift ratios under different requirements.

Graphical abstract

Keywords

interstory drift ratio / aggregation function / bound formulation / min–max problem / topology optimization

1 Introduction

In the design of tall building structures, the main effort is devoted to repetitively adjusting the layouts and geometry sizes of the structural members to satisfy the specified design criteria. Among them, the earthquake-induced maximum interstory drift ratio must not exceed certain limits stipulated in global seismic design codes [1–4]. Exceeding these limits is regarded as failure in achieving the performance objectives. The interstory drift ratio is widely used as a reliable indicator of structural state owing partly to its simplicity and convenience but primarily to its capability in measuring structural deformation, which is closely associated with structural and nonstructural damage caused by different levels of earthquake motions [5,6] and wind loads [7]. In addition to its correlation with nonstructural damage to partitions and cladding, the interstory drift ratio is also correlated with the lateral stiffness distribution of building structures. A story with insufficient lateral stiffness exhibits a high interstory drift ratio and becomes a weak story, resulting in severe structural damage under earthquake excitations [5,8]. Once the weak story yields, structural deformation gradually concentrates on the softened weak story, thereby deterring the distribution of plastic hinges to other stories. This causes a collapse mechanism in ductile reinforced concrete (RC) frames, as observed in the 1999 Izmit earthquake (in Turkey) [9] and the 2017 earthquake (in Mexico) [10]. Seismic retrofitting guidelines are utilized to address this issue by enhancing the strength and stiffness of the first floor [2,5]. However, this approach does not necessarily reduce the expected total damage and building losses because the damage to the structure may be shifted to the upper floors after the weak story is strengthened [11].

If a multistory structure deforms uniformly, a weak story and concentrated damage can be avoided; hence, structural safety is improved. Researchers have focused on the effect of vertical continuous stiffness in mitigating damage concentration along the height of a structure. Lai and Mahin [12] used a mast as a strong support in concentrically braced frames to promote uniform interstory drifts along building heights. Alavi and Krawinkler [13] proposed using hinged walls to engage frames to reduce the maximum interstory drift demands, which ultimately resulted in more uniformly distributed interstory drift ratios along the height. Moghaddam et al. [14] reported that the optimized structures with uniformly distributed interstory drift ratios underwent less seismic damage than conventional structures. These studies show that a uniform interstory drift ratio distribution can effectively reduce seismic damage. However, identifying an appropriate global design encompassing all structural members to achieve the desired interstory drift ratio distribution is more difficult than designing structural member strengths at the local level. This is because any change in either the structural member size or structural layout to adjust an interstory drift ratio at one story may affect the other interstory drift ratios. Hence, controlling the interstory drift ratio distribution requires the simultaneous consideration of all elements in a structure, which is time consuming for structural engineers, who typically address this issue by performing the conventional trial-and-error and analysis–design processes. The resulting final designs may be feasible and reasonable but tend to be conservative and nonoptimal in terms of cost and performance [15].

To improve design efficiency and quality, researchers have exploited optimization techniques and their applications to create excellent structural designs for building engineering [16–20] and multiphysics problems [21–24]. In terms of optimizing the distribution of interstory drift ratios, researchers [16–20] have primarily performed the formulation, where interstory drift ratios are limited via a series of constraint functions and the structural weight or material consumption is minimized. Although this formulation enforces optimized interstory drift ratios that are lower than or equal to the prescribed limits, it does not theoretically guarantee that the resulting interstory drift ratios exhibit a certain pattern. Another shortcoming becomes evident when the design variables and floor number increase. For example, the floor number may exceed 100 in a typical skyscraper [25]. In such circumstances, the computational efficiency is severely impaired owing to the significant computational burden arising from a sensitivity analysis that involves numerous constraint functions for identifyinga feasible path that simultaneously satisfies all constraints. The situation may be worse when performing topology optimization for building structures, where the number of design variables is typically of the order of magnitude of approximately 10⁴ [26].

An alternative optimization strategy is to minimize the maximum interstory drift or interstory drift ratio. This strategy has been reported in only a few studies [27–29], where the researchers successfully demonstrated that minimizing the maximum interstory drift can transform the interstory drift profiles into a uniform shape. Therefore, the strategy is superior to one that limits interstory drifts via constraint functions. Xu et al. [27] attempted to minimize the maximum variance of interstory drift in an optimization problem targeting nonlinear structures subject to stochastic dynamic loading. The nongradient-based optimization algorithm employed in that study circumvented the crucial issue wherein the objective function, which is a max operator, is not differentiable with respect to the design variables. Wu et al. [28] suggested using the sum of all interstory drifts as an objective function, which is smooth and differentiable, to optimize story stiffness; however, the optimization is driven by a hybrid heuristic algorithm. Although minimizing the sum of interstory drifts is an indirect approach for minimizing the maximum interstory drift, it is not equivalent to minimizing the maximum interstory drift in all cases. This articulation is strongly supported by the numerical examples presented in Ref. [29]. Meanwhile, Gomez et al. [29] adopted a smooth objective function that utilizes the Kreisselmeier–Steinhauser (K–S) function [30] as an alternative to the maximum function and focused on integrating stochastic excitation into topology optimization frameworks. However, the effects of substituting a smooth objective function to an original min–max problem remain ambiguous because the details of the optimization process are not completely divulged.

Although uniformly distributed lateral drifts are preferred in multistory frame structures comprising primarily columns and beams, this characteristic is not suitable for tall or supertall buildings because of their complex structural systems. For instance, RC core tubes in supertall buildings gradually degenerate into space frames, resulting in different lateral deformation capacities along their height. For the lower stories, where the RC core tubes maintain an intact multitube layout, the allowable elastic interstory deformation is less than that for the top stories, where structural walls do not exist, and braced frames dominate the structural deformation characteristics. This indicates that a reasonable interstory drift ratio distribution for a high-rise building should vary along the height instead of being uniform at all locations. In addition, for a hybrid building structure featuring a five-story steel frame erected on top of another fivestory RC frame, enforcing a uniform distribution of interstory drift ratios is inappropriate meaningless. This is because different materials feature different lateral deformation capacities in structural members. According to the Chinese seismic design code [1], the elastic interstory drift ratio limit for steel frames is 1/250, whereas that for RC frames is 1/550. This discrepancy will widen further if the five-story steel frame is assembled as a self-centering structure [31,32].

Hence, optimization techniques must be developed to accommodate different design expectations for the distribution of interstory drift ratios. This aim of this study is to establish an optimization scheme based on continuum topology optimization [33,34] that allows interstory drift ratio profiles to be transformed into prescribed shapes, which can be uniform, inverted triangular, or quadratic. Simultaneously, the maximum interstory drift ratio will reduce gradually during optimization. Two formulations that are not affected by the nondifferential maximum function are proposed. Thus, gradient information can be readily obtained to establish a highly efficient gradient-based algorithm for topology optimization. The optimization details are presented comprehensively to provide a foundation and insights for future studies. Meanwhile, the optimization results are compared to confirm the effectiveness of the proposed optimization scheme.

2 Min–max problem, bound formulation, and aggregation function

2.1 Min–max problem

Minimizing the maximum interstory drift ratio is a min–max problem. Generally, a min–max problem can be expressed mathematically as follows:

(1)\begin{array}{r}\underset{\mathit{\rho }}{min}:max\left\{f_1\left(\mathit{\rho }\right),f_2\left(\mathit{\rho }\right),\dots ,f_i\left(\mathit{\rho }\right),\dots ,f_n\left(\mathit{\rho }\right)\right\},\\ \left(i\in 1,2,\dots ,n\right)\end{array}

where \mathit{\rho } represents the set of design variables, and f_i\left(\mathit{\rho }\right) indicates a continuous scalar function. Various min–max problems exist in mechanical and civil engineering. For instance, suppressing the maximum nominal stress in different load cases can prevent static fracture and dynamic fatigue failure. In addition, advanced applications [35,36] utilize a min–max formulation to achieve robust optimization results. The min–max problem expressed in Eq. (1) is nondifferentiable with respect to the design variables. A straightforward example is the response spectrum of a natural earthquake ground motion [37]. In this example, the spectral curve is continuous but not smooth, and the period of a single-degree-of-freedom system is regarded as a unique design variable. When the design variables vary continuously, \mathrm{m}\mathrm{a}\mathrm{x}\{f_1\left(\mathit{\rho }\right),\cdots ,f_n\left(\mathit{\rho }\right)\} yields a continuous curve enveloping all scalar functions. Nondifferentiable points on the curve should not be dismissed because one of them may be the maximal value, which should be minimized using an optimizer. This nondifferential property does not pose any difficulty for nongradient algorithms but will hinder the use of gradient-based algorithms for efficient topology optimization.

2.2 Bound formulation

The bound formulation was introduced by Bendsøe et al. [38] for circumventing the nondifferential property of a max operator, and it can be expressed as follows:

(2)\begin{array}{c}\underset{\mathit{\rho },\tau }{min}:\tau \\ \mathrm{s}.\mathrm{t}.:f_i\left(\mathit{\rho }\right)-\tau ⩽0\left(i\in 1,2,\dots ,n\right),\end{array}

where \mathit{\rho } represents the set of design variables, and f_i\left(\mathit{\rho }\right) indicates a continuous scalar function. In the bound formulation, a new design variable \tau serves as the objective function and simultaneously represents a shared upper bound over all the scalar functions. Each scalar function and its upper bound are combined to establish an independent inequality constraint. The number of constraints is n. Consequently, minimizing \tau is mathematically equivalent to minimizing the maximum value among all the scalar functions. The new objective function is linear and its sensitivity is trivial and always equal to 1.

2.3 Aggregation function

Another approach to solving the min–max problem is to replace the original maximum function with an aggregation function. This approach is widely applied in optimization problems involving multiple load cases [39] or the depression of the maximum stress in the material [40]. Typically, two types of aggregation functions exist. One is the p-norm function [40,41], and the other is the K–S function [29,39]. Without loss of generality, the p-norm function is adopted in this study to control the interstory drift ratio distribution in the optimization process. The adopted p-norm function is expressed as

(3)\begin{array}{c}{f}_{p-n}={\left\{\sum _{i=1}^n\left(f_i^p\right)\right\}}^{1/p}\\ f_i=f_i\left(\mathit{\rho }\right)⩾0\left(i\in 1,2,\dots ,n\right),\end{array}

where p is the aggregation parameter and should be positive if f_{p-n} approximates the maximum value among all f_i\left(\mathit{\rho }\right). If each f_i\left(\mathit{\rho }\right) is smooth, then f_{p-n} is a smooth function. The value of f_{p-n} is always greater than that of max{f_i\left(\mathit{\rho }\right)}, and the following inequality relationships exist:

(4)max\left\{f_i\left(\mathit{\rho }\right)\right\}⩽f_{p-n}⩽max\left\{f_i\left(\mathit{\rho }\right)\right\}\cdot n^{1/p}.

Clearly, f_{p-n} converges to max{f_i\left(\mathit{\rho }\right)} as p approaches positive infinity because n^{1/p}\to1 if p\gg n. The inequality relationship shown in Eq. (4) is derived based on the upper and lower bounds of the K–S function [42]. In principle, a greater aggregation parameter typically leads to a more compact envelope approximation; however, an extremely high value of the aggregation parameter may cause numerical instability.

3 Optimization scheme for interstory drift ratio control

3.1 Interstory drift ratio

Interstory drift is defined herein as the lateral displacement of a floor relative to the floor directly beneath it, and the interstory drift ratio can be obtained by dividing the interstory drift by the vertical distance between consecutive floors. In a finite element model of a building structure, the interstory drift ratio of the kth floor can be calculated as follows:

(5)\begin{array}{r}{\theta }_k=\frac{1}{h_k}\left(\frac{1}{D_k}\sum _{i\in {\mathrm{\Omega }}_k}\left(u_i\right)-\frac{1}{D_{k-1}}\sum _{j\in {\mathrm{\Omega }}_{k-1}}\left(u_j\right)\right)={\mathit{L}}_k^{\mathrm{T}}\mathit{U},\\ \left(k=1,2,\dots ,n_{\mathrm{f}}\right)\end{array}

where {\theta }_k is the interstory drift ratio of the kth floor, h_k the vertical distance between the kth and (k − 1)th floors, n_{\mathrm{f}} the total floor number, D_k the number of nodes in the kth floor system, {\mathrm{\Omega }}_k the set of lateral degrees of freedom in the kth floor system, and u_i the lateral displacement corresponding to the ith degree of freedom. The expression for {\theta }_k can be simplified into a vector form by introducing a constant vector {\mathit{L}}_k because D_k, u_i, and the degree of freedom are known. Thus, {\theta }_k can be calculated using the matrix product of {\theta }_k = {\mathit{L}}_k^{\mathrm{T}}\mathit{U}, which offers brevity and elegance to the subsequent formulas used in the sensitivity analysis. When k = 1, {\mathrm{\Omega }}_{k-1} becomes {\mathrm{\Omega }}_0. In this case, Eq. (5) holds because {\mathrm{\Omega }}_0 represents the set of lateral degrees of freedom at the ground level, where all the degrees of freedom are imposed with fully fixed constraints and have no displacement.

3.2 Topology optimization scheme

A density-based topology optimization approach [43,44] was employed to perform topology optimization. Therefore, the structural material in the design domain was discretized based on the following scheme: In the topology design domain discretized by finite elements, a design variable {\rho }_e was assigned to the eth element. Based on the rational approximation of the material properties (RAMP) model [45], the Young’s modulus of the eth element is correlated with the projected physical density {\overline{\rho }}_e\in \left[0,1\right] as follows:

(6)E_e=E_{\mathrm{m}\mathrm{i}\mathrm{n}}+\frac{{\overline{\rho }}_e}{1+q\left(1-{\overline{\rho }}_e\right)}\left(E_0-E_{\mathrm{m}\mathrm{i}\mathrm{n}}\right),

where q \ge 0 is the penalization parameter; E_0 is the Young’s modulus of the structural material; E_min is an extremely low value assigned to void regions to prevent the global stiffness matrix from becoming singular; and {\overline{\rho }}_e is the projected physical density [35], which is calculated as

(7){\overline{\rho }}_e=\frac{\tanh \left({\beta }_{\mathrm{f}}\eta \right)+\tanh \left({\beta }_{\mathrm{f}}\left({\tilde{\rho }}_i-\eta \right)\right)}{\tanh \left({\beta }_{\mathrm{f}}\eta \right)+\tanh \left({\beta }_{\mathrm{f}}\left(1-\eta \right)\right)},

where \eta and {\beta }_{\mathrm{f}} determine the shape of the curve expressed in Eq. (7) via the optimization process. {\overline{\rho }}_e is the filtered density of the eth element and is calculated by filtering the design variables as follows:

(8){\tilde{\rho }}_e=\frac{1}{{\displaystyle \sum _{i\in n_e}{H}_{e,i}}}\sum _{i\in n_e}\left({\rho }_i{H}_{e,i}\right),

where n_e is the set of elements i for which the center-to-center distance Δ(e, i) to the eth element is smaller than the filter radius r_min, and H_e,i is a weight factor defined as H_e,i = max{0, r_min-\left(Äe,i\right)}. Thus, the element stiffness matrix can be correlated with the design variables using the interpolated Young’s modulus E_e in Eq. (6). Hence, the global stiffness matrix K can be obtained by assembling the stiffness matrix of each element. The interpolated volume of the eth element V_e is determined based on a linear relationship with respect to the physical volume of the eth element, {\overline{V}}_e, as V_e = {\overline{\rho }}_e{\overline{V}}_e.

3.3 Formulation for controlling interstory drift ratio distribution

The two formulations proposed herein are intended to transform the interstory drift ratio distribution of an optimized structure into a prescribed profile. The first formulation based on the bound formulation is as follows:

(9)\begin{array}{rl}& \underset{\mathit{\rho },\tau }{min}:\tau \\ & \mathrm{s}.\mathrm{t}.:\{\begin{array}{l}{g}_k={\displaystyle \frac{{\theta }_k^2}{{\omega }_k^2}}-\tau ⩽0,\left(k=1,2,\dots ,n_{\mathrm{f}}\right)\\ \frac{{\displaystyle \sum _e{V}_e}}{{\displaystyle \sum _e{\overline{V}}_e}}⩽V_{\mathrm{f}\mathrm{r}\mathrm{a}},\end{array}\\ & \mathit{K}\left(\overline{\mathit{\rho }}\right)\mathit{U}={\mathit{P}}_{\mathrm{L}},\\ & \mathit{W}=\left[{\omega }_1,{\omega }_2,\dots ,{\omega }_k,\dots ,{\omega }_{n_{\mathrm{f}}}\right],\end{array}

where {\omega }_k is the kth shape weight coefficient, V_{\mathrm{f}\mathrm{r}\mathrm{a}} the volume fraction used to limit the consumption of the structural material, K the global stiffness matrix, and U the displacement vector under an external load {\mathit{P}}_{\mathrm{L}}. The inequality constraints impose a square form on the interstory drift ratios, thus ensuring that the constraints are effective regardless of whether the interstory drift ratio is positive or negative. The proportional relationship of the entries in W determines the distribution of the optimized interstory drift ratios. In principle, the optimized interstory drift ratio distribution is consistent with the proportional relationship between the shape-weight coefficients. For instance, when all entries in W are equal to 1, the formulation shown in Eq. (9) results in a uniformly distributed profile of the interstory drift ratios. If the entries in W linearly increase in the order of their floor sequence, then the optimized interstory drift ratio distribution exhibits a triangular shape. Theoretically, an arbitrary shape can be specified by adjusting the entry value of W. In this regard, a volume constraint is necessary. If the volume constraint is abandoned, then a trivial optimal solution can be a design domain fully occupied with structural materials at all locations for static structural optimization.

The second proposed formulation utilizes the p-norm function in Eq. (3) to construct the objective function, which is expressed as follows:

(10)\begin{array}{rl}& \underset{\mathit{\rho }}{min}:f_{p-n}={\left\{\sum _{k=1}^{n_{\mathrm{f}}}\left(f_k^p\right)\right\}}^{1/p},\\ & \mathrm{s}.\mathrm{t}.:\frac{{\displaystyle \sum _e{V}_e}}{{\displaystyle \sum _e{\overline{V}}_e}}⩽V_{\mathrm{f}\mathrm{r}\mathrm{a}},\\ & \mathit{K}\left(\overline{\mathit{\rho }}\right)\mathit{U}={\mathit{P}}_{\mathrm{L}},\\ & f_k={\left(\frac{{\theta }_k}{{\omega }_k}\right)}^2,\left(k=1,2,\dots ,n_{\mathrm{f}}\right)\end{array}

where p is the aggregation parameter, and the other symbols have the same definitions as previously mentioned.

Unlike Eq. (9), where multiple constraint functions must be included, Eq. (10) only requires one constraint function, i.e., the volume constraint. Hence, the computational cost for addressing the constraint functions in Eq. (9) can be avoided. The shape weight coefficients {\omega }_k \left(k=1,2,\dots ,n_{\mathrm{f}}\right) in Eq. (10) have the same effect as that used in Eq. (9). When {\omega }_k = 1 for all k, and p = 1, the objective function in Eq. (10) becomes the sum of the interstory drift ratios and is equivalent to the objective function suggested by Wu et al. [28]. The optimized interstory drift ratio profile should be consistent with the prescribed shape weight coefficients provided that the value of p is sufficiently large.

3.4 Sensitivity analysis

As mentioned in Subsection 2.2, the sensitivity of the objective function in Eq. (9) is straightforward. As for the constraint functions imposed on the interstory drift ratios, the adjoint method [46,47] is employed to avoid directly solving the displacement sensitivities. The kth constraint function in Eq. (9) can be augmented by the static equilibrium equation \mathit{K}\mathit{U} = {\mathit{P}}_{\mathrm{L}} as follows:

(11)g_k=\frac{{\theta }_k^2}{{\omega }_k^2}-\tau +{\mathit{\lambda }}_k^{\mathrm{T}}\left({\mathit{P}}_{\mathrm{L}}-\mathit{K}\mathit{U}\right),

where {\mathit{\lambda }}_k is an adjoint vector unrelated to the design variables. Taking the derivative of Eq. (11) with respect to the projected physical density {\overline{\rho }}_e yields

(12)\begin{array}{rl}\frac{\text{d}{g}_k}{\text{d}{\overline{\rho }}_e}=& \frac{2{\theta }_k}{{\omega }_k^2}{\mathit{L}}_k^{\mathrm{T}}\frac{\text{d}\mathit{U}}{\text{d}{\overline{\rho }}_e}+{\mathit{\lambda }}_k^{\mathrm{T}}\left(-\frac{\text{d}\mathit{K}}{\text{d}{\overline{\rho }}_e}\mathit{U}-\mathit{K}\frac{\text{d}\mathit{U}}{\text{d}{\overline{\rho }}_e}\right)\\ =& \left(\frac{2{\theta }_k}{{\omega }_k^2}{\mathit{L}}_k^{\mathrm{T}}-{\mathit{\lambda }}_k^{\mathrm{T}}\mathit{K}\right)\frac{\text{d}\mathit{U}}{\text{d}{\overline{\rho }}_e}-{\mathit{\lambda }}_k^{\mathrm{T}}\frac{\text{d}\mathit{K}}{\text{d}{\overline{\rho }}_e}\mathit{U},\end{array}

(13)\mathit{K}{\mathit{\lambda }}_k=\frac{2{\theta }_k}{{\omega }_k^2}{\mathit{L}}_k.

The right-hand side of Eq. (13) represents the adjoint load, whereas Eq. (13) is known as the adjoint equation. By solving the adjoint equation to yield {\mathit{\lambda }}_k, the terms in Eq. (12), which includes \mathrm{d}\mathit{U}/\mathrm{d}{\overline{\rho }}_e, disappears. Finally, \mathrm{d}{g}_k/\mathrm{d}{\rho }_e is determined using the chain rule, as follows:

(14)\frac{\text{d}{g}_k}{\text{d}{\rho }_e}=\sum _{i\in n_e}\left(\frac{\text{d}{g}_k}{\text{d}{\overline{\rho }}_i}\frac{\text{d}{\overline{\rho }}_i}{\text{d}{\tilde{\rho }}_i}\frac{\text{d}{\tilde{\rho }}_i}{\text{d}{\rho }_e}\right).

Similarly, the objective function in Eq. (10) can be augmented using the static equilibrium equation, \mathit{K}\mathit{U} = {\mathit{P}}_{\mathrm{L}}, as follows:

(15)f_{p-n}={\left\{\sum _{k=1}^{n_{\mathrm{f}}}\left(f_k^p\right)\right\}}^{1/p}+{\mathit{\lambda }}^{\mathrm{T}}\left({\mathit{P}}_{\mathrm{L}}-\mathit{K}\mathit{U}\right),

where \mathit{\lambda } is an adjoint vector unrelated to the design variables. Taking the derivative of Eq. (15) with respect to the projected physical density {\overline{\rho }}_e yields

(16)\begin{array}{rl}\frac{\text{d}{f}_{p-n}}{\text{d}{\overline{\rho }}_e}=& \sum _{k=1}^{n_{\mathrm{f}}}\left(\frac{\text{d}{f}_{p-n}}{\text{d}{f}_k}\frac{\text{d}{f}_k}{\text{d}{\overline{\rho }}_e}\right)+{\mathit{\lambda }}^{\mathrm{T}}\left(-\frac{\text{d}\mathit{K}}{\text{d}{\overline{\rho }}_e}\mathit{U}-\mathit{K}\frac{\text{d}\mathit{U}}{\text{d}{\overline{\rho }}_e}\right)\\ =& \left[\sum _{k=1}^{n_{\mathrm{f}}}\left(\frac{\text{d}{f}_{p-n}}{\text{d}{f}_k}\frac{\text{d}{f}_k}{\text{d}\mathit{U}}\right)-{\mathit{\lambda }}^{\mathrm{T}}\mathit{K}\right]\frac{\text{d}\mathit{U}}{\text{d}{\overline{\rho }}_e}-{\mathit{\lambda }}^{\mathrm{T}}\frac{\text{d}\mathit{K}}{\text{d}{\overline{\rho }}_e}\mathit{U}\end{array},

(17)\mathit{K}\mathit{\lambda }=\sum _{k=1}^{n_{\mathrm{f}}}{\left(\frac{\text{d}{f}_{p-n}}{\text{d}{f}_k}\frac{\text{d}{f}_k}{\text{d}\mathit{U}}\right)}^{\text{T}}=\sum _{k=1}^{n_{\mathrm{f}}}\left(\frac{\text{d}{f}_{p-n}}{\text{d}{f}_k}\frac{2{\theta }_k}{{\omega }_k^2}{\mathit{L}}_k\right).

The right-hand side of Eq. (17) represents the adjoint load, and Eq. (17) is the adjoint equation. By solving the adjoint equation to yield \mathit{\lambda }, the terms in Eq. (16), which includes \mathrm{d}\mathit{U}/\mathrm{d}{\overline{\rho }}_e, disappears. Finally, \mathrm{d}{f}_{p-n}/\mathrm{d}{\rho }_e is determined using the chain rule, as follows:

(18)\frac{\text{d}{f}_{p-n}}{\text{d}{\rho }_e}=\sum _{i\in n_e}\left(\frac{\text{d}{f}_{p-n}}{\text{d}{\overline{\rho }}_i}\frac{\text{d}{\overline{\rho }}_i}{\text{d}{\tilde{\rho }}_i}\frac{\text{d}{\tilde{\rho }}_i}{\text{d}{\rho }_e}\right).

Performing a sensitivity analysis of the objective function in Eq. (10) solves the adjoint problem only once; however, the adjoint problem presented in Eq. (9) must be solved n_{\mathrm{f}} times when performing the sensitivity analysis. Comparing Eqs. (13) and (17), the adjoint load represented in Eq. (17) can be regarded as the weighted sum of all adjoint loads in Eq. (13).

4 Numerical investigation

4.1 Structural model

In this study, a numerical structural model was developed to investigate the performance of the proposed formulations. The model exhibits a plane-frame structure with a topological design domain, as illustrated in Fig.1(a). Similar testing models have been reported in the literature [29,48,49]. The plane frame subjected to lateral loads comprised three bays and nine stories. Each bay spanned 9 m in the lateral direction, with a floor height of 4.5 m. The entire design domain was modeled by a 27 m × 40.5 m rectangular geometry with a thickness of 0.04 m. The geometric information of the beam and column sections is provided in Tab.1 and Tab.2. The structural material properties of the plane frame and design domain reflect those of steel: Young’s modulus E_0 = 200 GPa and Poisson’s ratio \nu = 0.23. Two-node Euler beam elements were employed to model the columns and beams in the frame. Each node possessed two transitional degrees of freedom and one rotational degree of freedom. The design domain was uniformly discretized using square -node quadrilateral (Q4) in-plane stress elements with a 240 × 360 mesh grid. The lateral and vertical degrees of freedom at the beam-to-column nodes in the frame system were fully coupled to those (in the Q4 elements) at the same position in the design domain. Consistent with the results of previous studies [14,50], the earthquake-induced load pattern, as shown in Fig.1(b), demonstrated a linear incremental behavior along the height. Lateral forces were imposed on the beam-to-column nodes, and each component was calculated as follows: Q = \overline{Q}{h}_{\mathrm{Q}}, where \overline{Q} = 50 kN/m, and h_{\mathrm{Q}} is the height of the node position. Because the vertical loads did not significantly affect the lateral drifts in the optimized results, they are not considered herein. The column bases in the frame and the bottom boundary of the design domain were completely restrained. Evidently, both the structure and loads shown in Fig.1 are symmetric.

By controlling the axial stiffness of the beam elements in the frame, we implemented two diaphragm models in this study, that is, a rigid and a flexible diaphragm model. In the rigid diaphragm model, the beam-section areas are amplified by a sufficiently large factor such that the beams engage all columns at the same level, thus causing those columns to yield similar lateral drifts. In the flexible diaphragm model, the values listed in Tab.2 are adopted in the beam-section areas to achieve the actual axial stiffness of the beams.

In all the presented optimization cases, the penalization parameter q and E_min in Eq. (6) are 3 and 10⁻⁷E_0, respectively, and the filter radius r_min for Eq. (8) spans eight Q4 elements; V_{\mathrm{f}\mathrm{r}\mathrm{a}} in Eqs. (9) and (10) is 0.25, \eta in Eq. (7) is 0.5, and {\beta }_{\mathrm{f}} gradually increases to a maximum value of 50. Each optimization case begins from the same initial design, in which all the design variables are equal to 0.25. The method of moving asymptotes (MMA) proposed by Svanberg [51,52] was used as an optimizer to drive the optimization process. The convergence criterion specifies that the iteration continues until the maximum change among the updated design variables is less than 0.03.

4.2 Effect of uniform distribution control on interstory drift ratios

4.2.1 Effect of aggregation parameter

To investigate the effect of the aggregation parameter p in Eq. (10) on the optimization results, all the shape weight coefficients were set to 1 to obtain a uniform profile of interstory drift ratios, and the value of p was varied from 1 to 24 at a sampling interval of 2. The projected physical density ({\overline{\rho }}_e) fields of the optimization results with flexible diaphragms are shown in Fig.2, whereas those with rigid diaphragms are presented in Fig.3.

In the optimized results, the structural layouts and members, which are represented by solid elements in the design domain, exhibit a common feature. Amega chevron brace appeared in the lower section from the first to fifth stories, but the upper section showed different structural layouts composed of slender tree branch-like members owing to the variation in the aggregation parameter p. When the value of p increased from 1 to 14, the corresponding optimized results, as shown in Fig.2(a)–Fig.2(h), indicated distinct changes. However, the optimized results exhibited similar layouts with negligible differences in the details at p values ranging from 16 to 24, as shown in Fig.2(i)–Fig.2(m), and presented a convergent tendency. The transition features shown in Fig.3 concerning the variation in the aggregation parameter p, exhibits a similar pattern. These optimization results indicate that excessively increasing the value of p is unnecessary, partly because the variation tendency of the optimized layouts has already converged when p reaches 24, and partly because an excessively large value of p would cause numerical instability.

4.2.2 Characteristics of optimization history

The convergence histories of the interstory drift ratio profiles are shown in Fig.4: each case exhibited an identical unevenly distributed interstory drift ratio profile at the beginning of the optimization process. The initial interstory drift ratio profile confirmed the characteristics of the frame structures, with higher interstory drift ratios located at the lower stories and smaller ones clustered at the upper stories. Because the initial element stiffness for the Q4 elements in the design domain was discounted by approximately 69% based on the RAMP interpolation scheme, the frame contributed more significantly to the global lateral stiffness than the lateral load-resisting system embedded in the design domain; therefore, it dominated the structural deformation. As the optimization proceeded, the unevenly distributed profile gradually varied. In the first variation phase, the upper interstory drift ratios increased, whereas the lower ones decreased. Subsequently, the upper interstory drift ratios stopped increasing; in fact, they decreased simultaneously with the lower ones in the second variation phase until the optimization was terminated.

Because the proposed formulations involved a minimization process, the optimization histories of the maximum interstory drift ratio, as illustrated in Fig.5, declined as the iterations proceeded in all cases, except for the case in which the aggregation parameter was 1. In this particular case, the maximum interstory drift ratio increased before the optimization was terminated, as shown in the inset of Fig.5(b). This clearly shows that minimizing the sum of the interstory drift ratios is not equivalent to minimizing the maximum interstory drift ratio. In fact, performing optimization to minimize the sum of the interstory drift ratios can magnify the maximum interstory drift ratio in some cases.

4.2.3 Comparison between two formulations

A comparison of the interstory drift ratio profiles obtained using the aggregation and bound formulations is shown in Fig.6. The results show that the aggregation formulation could not ensure that the optimized interstory drift ratio profiles converged to a perfectly uniform shape when the aggregation parameter p was equal to 1. This implies that minimizing the sum of the interstory drift ratios cannot yield a perfectly uniform distribution. The higher the value of p, the more uniform is the distribution of the interstory drift ratios. In addition, the converged maximum interstory drift ratios in each case are presented along with their aggregation parameters, as shown in Fig.7. As the value of p increased, the maximum interstory drift ratio decreased more readily, although this tendency was locally disturbed by the fluctuations at p = 12 and 14 (in Fig.7(a)) and p = 16 (in Fig.7(b)). These results indicate that the proposed bound formulation shown in Eq. (10) yielded the best performance, which featured an almost vertical profile and the same values in the first four digits of the interstory drift ratios, as listed in Tab.3. The bottom blue line in Fig.7 represents the maximum interstory drift ratio yielded by the bound formulation. Similarly, the bound formulation exhibited better performance.

Furthermore, a comparison between Fig.2(m) and Fig.2(n) shows that the aggregation formulation failed to achieve the same or similar optimization results as the bound formulation, even when the value of p reached 24. This proves that the two proposed formulations are not equivalent from the perspective of numerical implementation results.

4.3 Effect of nonuniform distribution control on interstory drift ratios

In the proposed formulations shown in Eqs. (9) and (10), the shape weight coefficients determine the optimized interstory drift ratio distribution. The effectiveness of those formulations in achieving a uniform interstory drift ratio distribution is validated in Subsection 4.2.1. To further validate the suggested control strategy for realizing a nonuniform distribution of interstory drift ratios, two other groups of {\omega }_k were tested, as listed in Tab.4. The purpose of the first coefficient group was to realize an inverted triangular profile for interstory drift ratios along the height, whereas that of the second coefficient group was to transform a curved interstory drift ratio distribution based on a quadratic polynomial. Nonuniform distribution control was performed to relax the deformation restriction on the upper stories. Because the rigid rotation affects the interstory drift ratios of the upper stories but is rarely associated with structural damage, the upper interstory drift ratios can be set slightly higher than the lower ones. Based on the effect of the aggregation parameter p mentioned in Subsection 4.2.1, the value of p was set to 24 to simultaneously achieve a firm control effect and numerical stability. For comparison, the bound formulation shown in Eq. (9) was used under the same conditions. The projected physical density ({\overline{\rho }}_e) fields of the optimization results are shown in Fig.8. Owing to the nonuniform distribution control, the structural layouts shown in Fig.8 differed significantly from those under uniform distribution control shown in Fig.3. The upper structural layouts showed the most significant difference, whereas the lower sections appeared similar.

The convergence histories of the interstory drift ratio profiles corresponding to the designs shown in Fig.8 are presented in Fig.9. Their variation patterns were the same, as illustrated in Fig.4, and two phases were indicated. A comparison between the converged interstory drift ratio profiles and the reference lines is shown in Fig.10. These reference lines strictly reflect the proportional relationship exhibited by the shape-weight coefficients listed in Tab.4. The profiles generated by the proposed bound formulation fully overlapped their reference lines, which reflects the excellent capability of the nonuniform distribution control on the interstory drift ratios. Unlike the bound formulation, the proposed aggregation formulation presents a discounted control effect because of the divergence (see Fig.10(a) and Fig.10(c)) between the optimized results and their reference lines, although the divergence is almost invisible in Fig.10(b) and Fig.10(d). In general, the converged profiles eliminated the initial-stage characteristics of the frame structures and successfully transformed into an inverted triangular or quadratic shape instead of a disordered one.

5 Discussion regarding effects of diaphragm models on optimization

The effects of diaphragm models have been widely discussed in previous studies [29,53,54]. This study revisits this topic by comparing the optimized results with those of flexible and rigid diaphragm models. For convenience, the comparative pairs and their information are listed in Tab.5. Generally, using different diaphragm models yields different results. This observation is consistent with those of previous studies [29,53,54]. Because the axial stiffness of the beams was unrealistically amplified, the finite element model with rigid diaphragms exhibited greater lateral stiffness than its counterpart with flexible diaphragms. This can be justified by comparing the interstory drift ratios shown in Fig.10(a) and Fig.10(c), as well as by comparing the interstory drift ratios shown in Fig.10(b) and Fig.10(d).

Another concern, which has not yet been reported in the literature, is related to the deviation in the interstory drift ratios of the optimized results when they are re-evaluated using a different diaphragm model. Notably, a realistic floor system comprising beams and slabs provides an axial stiffness between that of a flexible and a rigid diaphragm model.

To investigate the deviation, the typical optimized results were re-evaluated using both the flexible and rigid diaphragm models for comparison. For example, Fig.11(a) and Fig.11(c) show a comparison of the interstory drift ratios of the optimized results using the flexible diaphragm model with those re-evaluated using the rigid diaphragm model. Meanwhile, Fig.11(b) and Fig.11(d) show a comparison of the interstory drift ratios of the optimized results using the rigid diaphragm model with those re-evaluated using the flexible diaphragm model. In these examples, the deviation effect was more evident in the optimized results yielded by the rigid diaphragm model. Hence, the flexible diaphragm model is more suitable for optimization when the floor system is to be simplified to reduce the computational cost.

6 Conclusions

In this study, the interstory drift ratio profile was controlled using topology optimization strategies. The two formulations proposed herein were based on bound and aggregation formulations. By setting the shape weight coefficients in the proposed formulations, the interstory drift ratio profiles of the optimized structures can be transformed into any prescribed shape. The numerical examples validated the efficacy of the control in generating structural designs with uniform, inverted triangular, or quadratic distributions of interstory drift ratios.

The effect of the aggregation parameter on the p-norm function was comprehensively investigated. The optimized structural layouts varied significantly as the aggregation parameter increased from 1 to 14 and then gradually converged as the aggregation parameter increased to 24. When p was equal to 1, the aggregation formulation degenerated into the sum of all interstory drift ratios, which is not equivalent to the original min–max problem, that is, minimizing the maximum interstory drift ratio, and did not achieve a uniform distribution of interstory drift ratios in the numerical examples presented. A higher value of the aggregation parameter rendered the resulting interstory drift ratio profile more similar to the prescribed shape, which is beneficial for lower maximum interstory drift ratios.

In addition, the structural layouts based on the aggregation formulation differed from those generated by the bound formulation, even when the value of the aggregation parameter was sufficiently high. Although the two proposed formulations successfully controlled the optimized interstory drift ratios, the bound formulation yielded optimized results that satisfied the prescribed shape more closely than the aggregation formulation. However, the proposed aggregation formulation is more pragmatic for general tall buildings in terms of reducing computational costs.

The effects of the rigid and flexible diaphragm models on the optimized results were investigated via cross-validation. The results showed that the interstory drift ratios from the optimized results based on the flexible diaphragms deviated slightly from those based on the rigid diaphragm model. Therefore, using the flexible diaphragm model is more conservative and reliable than using the rigid diaphragm model in such an optimization.

We believe that the reliability and effectiveness of the formulations proposed herein are not confined to topology optimization problems. In fact, the proposed formulations can be used to optimize the size or shape of building structures that require the distribution control of interstory drift ratios.

References

[1] GB50011-2010. Code for Seismic Design of Buildings. Beijing: Ministry of Construction of China, 2010 (in Chinese)

[2] ATC. Guidelines for SEISMIC PERFORMANCE ASSessment of Buildings. Redwood City (CA): Applied Technology Council, 2007

[3] SNZ. Concrete Structures Standard. Wellington: Standards New Zealand, 2004

[4] ASCE/SEI7-02. Minimum Design Loads for Buildings and Other Structures. Reston, VA: American Society of Civil Engineers, 2016

[5] Moehle J P, Mahin S A. Observations on the behavior of reinforced concrete buildings during earthquakes. American Concrete Institute Special Publication, Earthquake-Resistant Concrete Structures—Inelastic Response and Design, 1991, 127: 67–90

[6] Mayes R L. Interstory drift design and damage control issues. Structural Design of Tall Buildings, 1995, 4(1): 15–25

[7] Griffis L G. Serviceability limit states under wind load. Engineering Journal AISC, 1993, 30(1): 1–16

[8] PaulayTPriestleyM J N. Seismic Design of reinforced Concrete and Masonry Buildings. New York: John Wiley and Sons, 1992

[9] Kirac N, Dogan M, Ozbasaran H. Failure of weak-storey during earthquakes. Engineering Failure Analysis, 2011, 18(2): 572–581

[10] Jara J M, Hernández E J, Olmos B A, Martínez G. Building damages during the September 19, 2017 earthquake in Mexico City and seismic retrofitting of existing first soft-story buildings. Engineering Structures, 2020, 209: 109977

[11] Agha BeigiHSullivanT JCalviG MChristopoulosC. Controlled soft storey mechanism as a seismic protection system. In: The 10th International Conference on Urban Earthquake Engineering. Tokyo: Tokyo Institute of Technology, 2013

[12] Lai J W, Mahin S A. Strongback system: A way to reduce damage concentration in steel-braced frames. Journal of Structural Engineering, 2015, 141(9): 04014223

[13] Alavi B, Krawinkler H. Strengthening of moment-resisting frame structures against near-fault ground motion effects. Earthquake Engineering & Structural Dynamics, 2004, 33(6): 707–720

[14] Moghaddam H, Hajirasouliha I, Doostan A. Optimum seismic design of concentrically braced steel frames: Concepts and design procedures. Journal of Constructional Steel Research, 2005, 61(2): 151–166

[15] Lagaros N D, Papadrakakis M. Seismic design of RC structures: A critical assessment in the framework of multi-objective optimization. Earthquake Engineering & Structural Dynamics, 2007, 36(12): 1623–1639

[16] Farahmand-Tabar S, Ashtari P. Simultaneous size and topology optimization of 3D outrigger-braced tall buildings with inclined belt truss using genetic algorithm. Structural Design of Tall and Special Buildings, 2020, 29(13): e1776

[17] Kim C K, Kim H S, Hwang J S, Hong S M. Stiffness-based optimal design of tall steel frameworks subject to lateral loading. Structural Optimization, 1998, 15(3−4): 180–186

[18] Chan C M, Zou X K. Elastic and inelastic drift performance optimization for reinforced concrete buildings under earthquake loads. Earthquake Engineering & Structural Dynamics, 2004, 33(8): 929–950

[19] Zou X K, Chan C M. An optimal resizing technique for seismic drift design of concrete buildings subjected to response spectrum and time history loadings. Computers & Structures, 2005, 83(19−20): 1689–1704

[20] Tomei V, Imbimbo M, Mele E. Optimization of structural patterns for tall buildings: the case of diagrid. Engineering Structures, 2018, 171: 280–197

[21] Vu-Huu T, Phung-Van P, Nguyen-Xuan H, Abdel Wahab M. A polytree-based adaptive polygonal finite element method for topology optimization of fluid-submerged breakwater interaction. Computers & Mathematics with Applications (Oxford, England), 2018, 76(5): 1198–1218

[22] Ghasemi H, Park H S, Rabczuk T. A multi-material level-set based topology optimization of flexoelectric composites. Computer Methods in Applied Mechanics and Engineering, 2018, 332: 47–62

[23] Ghasemi H, Park H S, Alajlan N, Rabczuk T. A computational framework for design and optimization of flexoelectric materials. International Journal of Computational Methods, 2020, 17(1): 1850097

[24] Hamdia K M, Ghasemi H, Zhuang X Y, Rabczuk T. Multilevel Monte Carlo method for topology optimization of flexoelectric composites with uncertain material properties. Engineering Analysis with Boundary Elements, 2022, 134: 412–418

[25] Zhang J, Li Q. Identification of modal parameters of a 600-m-high skyscraper from field vibration tests. Earthquake Engineering & Structural Dynamics, 2019, 48(15): 1678–1698

[26] Beghini L L, Beghini A, Katz N, Baker W F, Paulino G H. Connecting architecture and engineering through structural topology optimization. Engineering Structures, 2014, 59: 716–726

[27] Xu J, Spencer B F Jr, Lu X. Performance-based optimization of nonlinear structures subject to stochastic dynamic loading. Engineering Structures, 2017, 134: 334–345

[28] Wu S, He H, Cheng S, Chen Y. Story stiffness optimization of frame subjected to earthquake under uniform displacement criterion. Structural and Multidisciplinary Optimization, 2021, 63(3): 1533–1546

[29] Gomez F, Spencer B F Jr, Carrion J. Topology optimization of buildings subjected to stochastic base excitation. Engineering Structures, 2020, 223: 111111

[30] KreisselmeierGSteinhauserR. Systematic control design by optimizing a vector performance index. In: International Federation of Active Controls Symposium on Computer-Aided Design of Control Systems. Zurich: Pergamon Press Ltd., 1979

[31] Lu X, Cui Y, Liu J, Gao W. Shaking table test and numerical simulation of a 1/2-scale self-centering reinforced concrete frame. Earthquake Engineering & Structural Dynamics, 2015, 44(12): 1899–1917

[32] Gao W, Lu X. Modelling unbonded prestressing tendons in self-centering connections through improved sliding cable elements. Engineering Structures, 2019, 180: 809–828

[33] Bendsøe M P, Kikuchi N. Generating optimal topologies in structural design using a homogenization method. Computer Methods in Applied Mechanics and Engineering, 1988, 71(2): 197–224

[34] BendsøeM PSigmundO. Topology Optimization––Theory, Methods and Applications. Berlin: Springer, 2003

[35] Wang F, Lazarov B S, Sigmund O. On projection methods, convergence and robust formulations in topology optimization. Structural and Multidisciplinary Optimization, 2011, 43(6): 767–784

[36] Gao W, Wang F, Sigmund O. Systematic design of high-Q prestressed micro membrane resonators. Computer Methods in Applied Mechanics and Engineering, 2020, 361: 112692

[37] ChopraA K. Dynamics of Structures: Theory and Applications to Earthquake Engineering. New Jersey: Prentice-Hall, 1995

[38] Bendsøe M P, Olhoff N, Taylor J E. A variational formulation for multicriteria structural optimization. Journal of Structural Mechanics, 1983, 11(4): 523–544

[39] James K A, Hansen J S, Martins J R R A. Structural topology optimization for multiple load cases using a dynamic aggregation technique. Engineering Optimization, 2009, 41(12): 1103–1118

[40] Le C, Norato J, Bruns T, Ha C, Tortorelli D. Stress-based topology optimization for continua. Structural and Multidisciplinary Optimization, 2010, 41(4): 605–620

[41] Gao W, Lu X, Wang S. Seismic topology optimization based on spectral approaches. Journal of Building Engineering, 2022, 47: 103781

[42] Poon N M K, Martins J R R A. An adaptive approach to constraint aggregation using adjoint sensitivity analysis. Structural and Multidisciplinary Optimization, 2007, 34(1): 61–73

[43] Sigmund O A. 99 line topology optimization code written in Matlab. Structural and Multidisciplinary Optimization, 2001, 21(2): 120–127

[44] Andreassen E, Clausen A, Schevenels M, Lazarov B S, Sigmund O. Efficient topology optimization in Matlab using 88 lines of code. Structural and Multidisciplinary Optimization, 2011, 43(1): 1–16

[45] Stolpe M, Svanberg K. An alternative interpolation scheme for minimum compliance topology optimization. Structural and Multidisciplinary Optimization, 2001, 22(2): 116–124

[46] Feng T T, Arora J S, Haug E J. Optimal structural design under dynamic loads. International Journal for Numerical Methods in Engineering, 1977, 11(1): 39–52

[47] Arora J S, Haug E J. Methods of design sensitivity analysis in structural optimization. AIAA Journal, 1979, 17(9): 970–974

[48] Mijar A R, Swan C C, Arora J S, Kosaka I. Continuum topology optimization for concept design of frame bracing systems. Journal of Structural Engineering, 1998, 124(5): 541–550

[49] Stromberg L L, Beghini A, Baker W F, Paulino G H. Topology optimization for braced frames: combining continuum and beam/column elements. Engineering Structures, 2012, 37: 106–124

[50] Zhou Y, Zhang C, Lu X. An inter-story drift-based parameter analysis of the optimal location of outriggers in tall buildings. Structural Design of Tall and Special Buildings, 2016, 25(5): 215–231

[51] Svanberg K. The method of moving asymptotes—A new method for structural optimization. International Journal for Numerical Methods in Engineering, 1987, 24(2): 359–373

[52] Svanberg K. A class of globally convergent optimization methods based on conservative convex separable approximations. SIAM Journal on Optimization, 2002, 12(2): 555–573

[53] Allahdadian S, Boroomand B, Barekatein A R. Towards optimal design of bracing system of multi-story structures under harmonic base excitation through a topology optimization scheme. Finite Elements in Analysis and Design, 2012, 61: 60–74

[54] Allahdadian S, Boroomand B. Topology optimization of planar frames under seismic loads induced by actual and artificial earthquake records. Engineering Structures, 2016, 115: 140–154

RIGHTS & PERMISSIONS

The Author(s). This article is published with open access at link.springer.com and journal.hep.com.cn