Shape optimization of structures under earthquake loadings

The optimum design of structures under static loads is well-known in the design world; however, structural optimization under dynamic loading faces many challenges in real applications. Issues such as the time-dependent behavior of constraints, changing the design space in the time domain, and the cost of sensitivities could be mentioned. Therefore, optimum design under dynamic loadings is a challenging task. In order to perform efficient structural shape optimization under earthquake loadings, the finite element-based approximation method for the transformation of earthquake loading into the equivalent static loads (ESLs) is proposed. The loads calculated using this method are more accurate and reliable than those obtained using the building regulations. The shape optimization of the structures is then carried out using the obtained ESLs. The proposed methodologies are transformed into user-friendly computer code, and their capabilities are demonstrated using numerical examples.


Introduction
Most of the forces in the real world are dynamic in nature; in addition, their magnitudes are variable in the time domain. A dynamic analysis of structures requires the large amounts of information processing and data interpretation. Especially, when the optimization of structures under earthquake loading is required, the data processing is a cumbersome task. Accordingly, static loads could be utilized as a substitute for earthquake loadings, if they produce the same responses as the dynamic loads at the arbitrary time. Because of the facilities of the static analyses, the users of building codes have an interest to apply them in their design purposes. Application of dynamic coefficients or factors is a common way for transformation of dynamic loadings into static ones. However, the dynamic factors are not based on mathematical logic. They are mostly determined by engineering judgments and experience. Usually, these coefficients produce over-estimate loads and render the designs uneconomical (Humar and Maghoub 2003).
In competitive world, optimum design has great importance in economic design of structures. Therefore, optimum design of structures under dynamic loading is active in the fields of structural design. Because of the time-dependent behavior of the constraints and sensitivity analysis difficulties, optimization for earthquake loadings is a heavy duty. Accordingly, researchers have worked for many years in this field to determine the simple methods for optimum design of structures under such loads (Kang et al. 2006). Some researchers have focused on the methods that directly deal with dynamic loadings (Cassis and Schmit 1976;Feng et al. 1977;Yamakawa 1981;Mills-Curran and Schmit 1985;Haftka 1989, 1991;Chahande and Arora 1994).
Mathematics-based optimization procedures involve sensitivity calculations. Several methods have been proposed for sensitivity analysis in the structural optimization field under transient loadings (Yamakawa 1981;Mills-Curran and Schmit 1985;Haftka 1989, 1991). Dynamic response optimization is still difficult due to the large amounts of computational time required for analyses and gradient calculations. For that reason, the methods have been limited to small-scale structures with few degrees of freedom. For large-scale structures such as buildings and dams, the optimum design under dynamic loadings seems to be impossible, because difficulties arise in treating time-dependent behavior of constraints and objective functions. When structural optimization problems under dynamic loads are substituted with static response optimizations, two aspects should be considered. The first is the reliable transformation of dynamic loadings into static loads and the second is appropriate application of the resulting static loads into the optimization procedure. Cheng and Truman (1983) studied the optimization of the structures by modal response spectrum method. Truman and Petruska (1991) applied optimality criterion techniques for optimization of two-dimensional structures using time history analysis of seismic excitations. Although, the above researches are limited to seismic loads, the transformation is based on experimental codes, and it lacks generality. Choi et al. (2005) and Choi and Park (1999a, b) proposed several methods to find reliable equivalent static loads from dynamic loading. Two kinds of ESLs are considered in the literature: exact and approximate. Although, the location of the approximated ESLs should be pre-assumed in these methods, and the ESLs are calculated at some peaks of important locations (Kang et al. 2001). The pre-assumptions for load locations create different ESLs. It is difficult to choose the locations and to calculate the sensitivities of the ESLs with respect to the design variables.
In order to overcome these difficulties, two methods have been proposed for dynamic response optimization (Choi and Park 2002). First, exact ESLs are calculated at the time intervals during the dynamic loading. Second, an ESL set is defined as a static load set to generate the same response field that occurs under the dynamic load at a certain time. As a result, ESLs are generated for all nodes, and the pre-assumed locations are not required. In addition, the ESLs are calculated for all time intervals (Choi and Park 2002;Park et al. 2005;Hong et al. 2010;Kim and Park 2010). However, the proposed methods are complicated for optimization procedure. For engineering applications, performing optimization for all time intervals is difficult and requires a large amount of calculations. In addition, the dynamic response optimization studies using the above methods are limited to small cases with few degrees of freedoms and for very simple loadings such as impacts. In this paper, the ESLs are calculated by an approximation method, and the optimization procedure is very easy to apply.
Recently, structural optimization under earthquake loadings has been investigated in the civil engineering field. To our knowledge, comprehensive studies for shape optimization under direct earthquake excitation have not been reported in the literature. Accordingly, the present paper, implements a displacement-based finite element approach for transformation of earthquake loading into equivalent static loads. In order to calculate the equivalent static loads from earthquake loading, a mathematics-based optimization technique is used. In addition, shape optimization is applied to the combination of gravity, hydrostatic and seismic loadings due to its importance in the structural design criteria. Here, two-dimensional continuous structures are optimized directly under earthquake loading without any assumptions or complicated formulations. Accordingly, the optimum design of large-scale structures under earthquake loadings can be obtained by the proposed method.

Transformation of earthquake loading into ESLs
In spite of static loads, the values and the directions of earthquake loadings are time variants. It means that the static load has a different effect as the dynamic one. Since the static loads are easy to utilize in the optimization procedure, earthquake loadings are transformed into equivalent static loads (ESLs). The ESLs are static loads that produce the same displacement field as the earthquake loadings.
The process of ESLs calculation is presented in this section. The dynamic equation of motion for structures using the finite element method is explained as (1) Where M, C, K are the mass, damping and stiffness matrices, respectively. Vectors of displacements, velocities and accelerations of a structure at the arbitrary time t a are defined by U,U,Ü, respectively.Ü g is the ground acceleration, and r is the influence vector that shows the direction of the applied loading on the structure. In static analysis using finite element methods, the equilibrium equation is defined as follows: Here, D, P show the vectors of static displacements and static loads, respectively. Equation (3) shows the relation between earthquake loading and exact equivalent static loads. According to this equation, the ESLs generate the same displacement field as that of the earthquake loading at the arbitrary time t a .
KU (t a ) = P In (3), U (t a ) = u 1 , u 2 , . . . u q , . . . u N is the vector of nodal displacements in dynamic situation, P = p 1 , p 2 , . . . p j , . . . p N is a vector of exact equivalent static loads. This equation shows that the ESLs exist for earthquake loading at the arbitrary time. Substituting (3) into (1), relation between exact ESLs and earthquake loading can be defined as (4).
Equation (4) explains that for obtaining the values of ESLs (P), the transient analysis of structure is required, and this analysis also needs extensive calculations and information processing. Especially for optimization purposes, the computations are not efficient. Hence, in order to perform shape optimization, the approximate ESLs are calculated and applied on the structure. Equation (1) or (4) are coupled equations, and are not easy to solve for large-scale systems. Direct solution of these equations is easily possible using high-speed computers and commercial softwares. However, in linear seismic analysis and design of structures, eigenvalues (ω) and eigenvectors ( ) are available data, and using these data for efficient solution of (1) is a logical strategy. Here, the modal approach (U = Y) is applied for transformation of (1) into uncoupled second-order equation as (5).
y n + 2ξ n ω nẏ n + ω 2 n y n = T n M n MrÜ g Where, = [ 1 , 2 , .... n ] shows the modal matrix, T n = [φ 11 , φ 12 , . . . φ 1n ] is the vector of normalized components for n-th mode shape, and Y = y 1 y 2 , . . . y nm refers to the vector of modal displacements. nm is the number of applied modes in the calculation of dynamic responses. y n ,ẏ n ,ÿ n , refer to the modal displacement, velocity and acceleration at n-th mode, respectively. ω n is the frequency of n-th mode. Damping ratio and the modal mass of the n-th mode are defined by ξ n , M n .
Because of the complex and irregular forming of earthquake loading, a numerical method should be applied for solving (5). For this purpose, the implicit β-Newmark average acceleration algorithm is utilized (Newmark 1959). At a desired time (t a ), the dynamic displacement of the q-th degree of freedom is determined as (6).
Here, N refers to the total number of degrees of freedom for structure, and φ qk is the normalized component of a mode shape matrix. Equation (6) shows the modal concept in the dynamic analysis of structures. The following conversion is used for efficient solution of this equation.
For the normalized values of a mode shape matrix, the k-th component of Y could be calculated as (8) Dynamic displacement at q-th degree of freedom can be computed as (9) by substituting (8) into (6).
Equation (9) can be defined as a system of simultaneous linear equations with N variables for P j . Direct solution of this equation to produce the same displacement field as that of earthquake loading is a tiresome task. As a result, this equation is modified to inequality equation as follows: In engineering designs, the approximate methods are utilized in calculation of the ESLs. There are two approximate approaches for ESLs computations: displacement-based and stress-based techniques. In the standard finite element analysis, stresses in the elements or Gauss points are obtained from differentiation of the displacements. For this reason, the accuracies of stresses are less than those of displacements. Accordingly, the displacement-based method is utilized in the ESLs calculations. In dynamic response of structures, higher frequencies have fewer effects on their responses. Thus, only first nm natural frequencies are considered in practical applications, therefore, (10) can be approximated to (11) Where nm refers to the number of the first applied modes, and ml indicates to the number of the locations of the ESLs in structure. In order to find the numerical values of the ESLs, the optimization problem should be solved as (12) (Choi et al. 2005).
Where, d q is the dynamic displacement of structure at the critical time for the q-th degree of freedom. Selecting the objective function as a square sum of P j is the best way to find the smallest loads. The obtained ESLs in optimization formulation create a displacement field equalling the dynamic loading at the critical times (Grandhi et al. 1986).
The critical times are those in which the selected degrees of freedom of a structure have a maximum response during earthquake loading. As the geometry of structure changes during the optimization iterations, the shape of the structure modifies. Subsequently, the extreme responses maybe occur at other critical times. Equation (12) shows the produced approximated ESLs are usually greater than corresponding dynamic loads at critical times. The obtained static loads are directly applied to the structures to perform static optimization.

Finite element model
Gravity, hydrostatic, and earthquake loadings have been included in the finite element model. In civil engineering, linear design of structures under earthquake loading is a common methodology. For that reason, linear behavior of the structure is taken into account here. Four-node quadrilateral solid elements are used for modelling the plane strain and plane stress type of structures (see Fig. 1). For numerical calculation of responses, the required vectors and matrices for the elements are calculated using standard finite element analysis (Cook et al. 2002). In order to compute the seismic responses of structures, Rayleigh damping is used as (13) (Chopra 2001).
Where, α, β are the constants. According to the (13), these coefficients can be determined from specified damping ratios ξ i , ξ j for the i-th and j-th modes, respectively (Chopra 2001).
If both modes are assumed to have the same damping ratio ξ , which is reasonable based on experimental data, the constants are obtained by solving the two algebraic equations. Dynamic loading is a single component of earthquake acceleration here. The first ten modes of vibration have been included in the modal response analysis of structures. The implementation of the above methodologies in the developed computer program has been verified based-on the Choi benchmarks (Choi et al. 2005).

Shape optimization under seismic excitation
Mathematically, it is impossible to optimize a large-scale structure under earthquake loadings because the involved functions are defined over the time domain, and sensitivity analysis is very difficult. As mentioned earlier, only smallscale problems have been solved or the dynamic loads have been transformed into static loads using dynamic factors. The optimization process under earthquake excitation can be divided into two parts: the analysis and the design domains. In the analysis domain, the earthquake loading is transformed into the equivalent static loads. The transformed ESLs are used in the design domain as the external forces. Optimum values are found in the design domain. If the design is changed, the characteristics of the structure are changed and the transformation process should be conducted again. Because of this, the entire optimization process circulates between the two domains. The iterative procedure between the two domains is defined as a design cycle. The design cycle is performed iteratively until the design variables converge. The flow diagram of the shape optimization of structures under earthquake loading is illustrated in Fig. 2.
Since a structure vibrates under an earthquake loading, tension and compression can be observed in the same node (see Fig. 3). The static load makes only tension or compression at the arbitrary node. It should be mentioned that at least two static loads are needed to simulate states (B) and (C), respectively. Consequently, multiple static loads should be utilized to represent the effect of an earthquake loading. Then, the critical load should be applied to perform the static optimization of the structures. The approximated ESL can reduce considerable calculations. Structural optimization under dynamic load seems to be quite difficult. By using this method, dynamic response optimization can be accomplished for structures. A dynamic load can be transformed to the multiple static loads. The multiple loads can be handled as a multiple loading condition in optimization procedure. As be mentioned, the earthquake loading is transformed into the equivalent static loads in the appropriate degrees of freedom. Thus, the dynamic response optimization problem is converted to the static optimization formulation. Each optimization problem can be written in the general form as (15). minimize f(X) subject to: Here, f(X) is the objective function, X L , and X U are the lower and the upper bounds of the vector of the design variables, and nc is the total number of constraints. Since the costs of structural responses and sensitivity analyses during optimization iterations are expensive, direct solution of (15) is not easy. Using the finite element method for structural response, explicit relations between design variables and behavioral constraints are not available. Hence, the constraints and the objective function are approximated to the form of Taylor series expansion. For instance, for stress vector (σ e ) at the Gauss or nodal points, implicit and explicit relations can be defined as (16).
Where E, B, d e are the constitutive matrix, straindisplacement matrix and vector of nodal displacement, respectively. X = {x 1 , x 2 , x 3 , .
. . x n } shows the vector of design variables in the optimization process, and X 0 is the vector of design variables in which the Taylor series are expanded around it. Consequently, each cycle of the static shape optimization problem is as the following explicit form.
Minimize weight subject to: (σ t ) i and (σ c ) i refer to the tension, and compression stresses, respectively. The allowable values of the tension and compression stresses are shown by σ t allw , σ c allw . Forward finite difference method is applied for gradient calculations of stresses.

Sensitivities calculations
To calculate the sensitivities using forward finite difference method, the finite element analysis should be performed for each design variable. In forward finite difference method, the derivatives are approximated from the exact response at the original design point X and at the perturbed design point X + δX as (18).
Here, R(X) refers to the response of the structure, and δx i indicates the predetermined step-size of each design variable. Finite difference methods are the easiest to apply, and they are attractive for many applications. When R(X) is known, application of (18) involves only one additional calculation of the responses at X + δx i . For a problem with n design variables, finite-difference derivative calculations require repetition of the analysis for n + 1 (using (18)) different design points. Finite-difference approximations might have accuracy problems. Two sources of errors should be considered whenever these approximations are used. First, the truncation error, which is a result of neglecting terms in the Taylor series expansion of the perturbed response. The second is the condition error, which is the difference between the numerical evaluation of the function and its exact value. These are two conflicting considerations. That is, a small step size δx i will reduce the truncation error, but may increase the condition error. To eliminate round-off errors due to approximations it is recommended to increase the step-size. Here, the step size is assigned to the δx i = 0.005x i .

Numerical studies
Shape optimization of two-dimensional structures under earthquake loading is studied here. The first case is a plane stress cantilever beam, and the second one is a plane strain concrete gravity dam. Vertical component of Kobe earthquake is used for the beam excitation, and horizontal component of El-Centro earthquake is applied for the dam excitation. Design variables in both cases are X = {h c , h m , h b }. Where h c , h m , h b are the widths of the structures at the reference points (two ends and middle), respectively.

Case 1
This case study is a weight minimization of a cantilever beam. According to Fig. 4, during the shape optimization process, length and thickness of the beam are fixed to 6.0 m Fig. 4 The geometry and finite element mesh of the cantilever beam for initial design and 0.30 m, respectively. The material properties for the beam are also presented in Table 1.
Here, f cc is the uniaxial compressive strength of concrete at the age of 28 days, E c is the Young's modulus of concrete, γ c weight specific of concrete, and υ is Poisson ratio. The beam widths (h c , h m , h b ) should not be less than 30 cm anywhere. Initial values of design variables are set to 1.50 m. The beam's support is excited in y direction by TAZ090 component of Kobe earthquake loading (Fig. 5). Damping ratios for all modes have been assigned to 5 percent.
Free end vertical displacements of the beam are calculated using dynamic analysis, and their values are illustrated in Fig. 6. According to this figure, the critical time at initial design is at 5.32 (sec). At this time, maximum negative displacement of the free end of the beam is -0.68 mm. At the next critical time (6.15 (sec)), maximum positive displacement of the beam is 0.71 mm. In the optimization procedure, these critical times produce the multiple ESLs. As a result, for the first cycle of the optimization process, the critical time is 6.15 (sec).
The following optimization problem is utilized to compute the ESLs of the beam under Kobe earthquake.
Finding P j to minimize P 2 j subject to Here i is the index for the numbers of ESLs in structure. The locations of d i are shown in Fig. 7. In order to show the details of calculations, linear constraints of (19), are As seen in Fig. 7, the values of ESLs (P j ) are calculated using the SQP algorithm as a solution algorithm of (19). In addition, nodal displacements of the beam under original loading and ESLs, and errors are summarized in Table 2. Node numbers in this table indicate to the locations of the ESLs. e.g node number 1 refers to the node under P 1 . Table 2 shows that the displacements under ESLs are greater than those in the dynamic status. The last column of the table shows the errors of the selected degrees of freedom. The average error for this case is 0.15. The values of error in compare with Choi et al. (2005) are relatively reasonable. As the number of nodes where the ESLs are applied increases, the errors decrease. According to (17), to evaluate the values of constraints in the shape optimization procedure, stress analysis is required. Comparison of principal stresses for selected nodes as seen in Fig. 8, is presented in Table 3. Whereas in the ESL loading for nodes As seen from the table, the value of error for node 8 is not acceptable becuase the characteristics of the earthquake loadings and the ESLs are not same. Earthquake loadings are imposed at the supports of structures, and the stiffness, damping, and mass components of structures have an effect on produced displacements. However, in the ELSs, the displacements are affected only by the stiffness part of the system, and deformations are usually made under concentrated loads. Therefore, to our knowledge, the stresses are not comparable in these cases.
According to (12), the values of displacements under ESLs are usually greater than those from original dynamic loading. In spite of displacements, for stresses this criterion is not satisfied. As seen from Table 3, at some nodes such as 9, 12, 13 and 14, values of principal stresses under ESLs are less than those under earthquake loading.
Applying the calculated ESLs in the desired nodes, as shown in Fig. 7, shape optimization of the beam under Kobe earthquake is transformed into a static optimization problem. For initial design variables, the shape optimization formulation is as (21).
For finding the new values of design variables, a numerical optimization algorithm is used to solve (21) (Arora 2004). Accordingly, the shape of the beam is modified as the design variables are changed. The optimization process converged to the optimum design after three cycles, as shown in Fig. 9.
In initial design, weight of the beam is 211.50 kN and in final iteration, it reaches to 79.90 kN.

Case 2
This case is a weight minimization of a concrete gravity dam. The shape of the dam, as shown in Fig. 10, is similar to that of the Koyna dam. Free-bord of the Koyna dam, the part higher than EL66.5 m, is not modeled here. The upper face of the dam is subjected to hydrostatic pressure. The water level is at EL66.5 m. Moreover, the foundation of the dam is rigid and has not been included in the finite element model. The mesh dimensions in the finite element analysis have been selected based on sensitivity analysis studies. In addition to the earthquake loading, the weight of the dam and hydrostatic pressure loading are included in the finite element model. The dam is subjected to the horizontal component of El-Centro earthquake as seen in Fig. 11. In this case, damping ratios for all modes have been designated to 4.5 percent. Initial values of design variables are set to 16.30, 43.25 and 70.20 m, respectively. Values of design variables (h c , h m , h b ) should not be less than 15 m anywhere. In Table 4, material properties for this case are presented, in which γ w is weight specific of water.
In Fig. 12, dynamic displacements of the crest of the dam under El-Centro earthquake are illustrated. Two critical times should be considered here. The first time is at 5.20 (sec) and the maximum positive displacement of the crest is 5.70 mm. The second one is at 6.12 (sec) and the maximum negative displacement of crest is −8.80 mm. For that reason, the calculation of the ESLs should be carried out at 6.12 (sec).
Values of the ESLs, displacements of the selected nodes under earthquake loading and ESLs, and also errors are presented in Table 5.
Solving (23)       loadings are not same. Accordingly, calculation of ESLs using stress-based finite element approach could be an alternative method. This method could be used as a new research methodology for the optimum design of structures. 5. In the previous researches, shape optimization has been done at each time intervals. In comparison with the past researches, using obtained ESLs by this method to perform the optimum design in the critical times is more efficient and cost effective. 6. The proposed method needs transient dynamic analysis of the structure in calculation of the ESLs in each optimization cycle. However, the optimum design of structure is independent of performing the dynamic analysis. Similarly, the building codes and static optimizations do not require performing a dynamic analysis. 7. Shape optimizations are achieved with only a few iterations, and the computations are efficient in comparison with past researches.
Open Access This article is distributed under the terms of the Creative Commons Attribution License which permits any use, distribution, and reproduction in any medium, provided the original author(s) and the source are credited.