Shape optimization of tall buildings cross-section: Balancing profit and aeroelastic performance

Summary Shape optimization is an effective tool to improve the aerodynamic performance of tall buildings by introducing minor modifications to the original project. Nevertheless, economic criteria demand efficient cross sections aiming at maximizing the building's profitability. These two contradictory criteria are commonly handled by adopting multi-objective optimization approaches seeking the definition of Pareto fronts. However, the aerodynamic nonlinear features of low-aspect-ratio cross sections typically adopted in architectural practice can cause wind-induced acceleration response surfaces over the considered design domain with multiple local minima that eventually lead to discontinuous Pareto fronts with non-convex regions. This study delves into this problem and proposes a design framework that effectively combines the reduced basis method with multi-objective optimization techniques to carry out the aerodynamic shape optimization using surrogates trained with CFD simulations. The ability of the optimization strategy to properly define the non-convex regions of discontinuous Pareto fronts is successfully leveraged by adopting the weighted min – max method.

The mutual interdependency between cross-section shape, building height, acceleration response, and profitability is certainly complex.In Kwok, [6] this design problem is explained in the following terms: "… an aerodynamically efficient plan shape may come at the cost of usable floor area, in turn, may require additional compensatory storeys, which obviously increase the wind loads and construction cost again.On one hand, a less-efficient planform … usually attracts lower rental or sales returns; on the other hand, additional upper level storeys have the potential to increase revenue.Therefore, it is not straightforward for engineers to decide on aerodynamic modifications."This problem can be tackled efficiently by adopting multi-objective optimum design techniques, which can help the designer to identify the solution that offers an adequate structural response for the intended uses of the building while maximizing profitability. [7]cent contributions where optimum design techniques have been incorporated in the shape design of tall buildings including their aerodynamic performance are Kareem et al., [8] who combined CFD simulations with optimization algorithms to carry out the aerodynamic optimization of tall buildings cross-sections; Bernardini et al., [9] where a Kriging surrogate model based upon two-dimensional (2D) URANS simulations was adopted to speed up the optimization of the corners of a high-rise building cross-section using evolutionary algorithms; Elshaer et al., [1] in which the geometry of the corners of a tall building is optimized using genetic algorithms and an artificial neural network is adopted for the evaluation of the objective functions based on three-dimensional (3D) large eddy simulations (LES) considering turbulent flow; and Ding and Kareem, [10] who adopted a multi-fidelity surrogate modeling approach based on CFD simulations for addressing the shape optimization of the cross section of tall buildings.An interesting review paper on the subject is Asghari Mooneghi and Kargarmoakhar, [11] where the main challenges in the optimum design of buildings are discussed.The nonlinear behavior of the aerodynamic responses when the tall building cross section is modified is reflected in the Pareto fronts obtained when solving multi-objective optimization problems.In the studies reported by Ding and Kareem [10] and Elshaer and Butsuamlak, [12] discontinuities in the Pareto fronts can be seen as a consequence of the aerodynamic nonlinear features.Addressing this interesting and complex problem requires the adoption of the most adequate formulation from the optimization perspective.
It must be also noted that, in the aforementioned references, the aerodynamic behavior of candidate designs or surrogate training samples is assessed using CFD simulations and the modifications of the building's cross-sectional geometry are ruled by a set of shape design variables that define the position of a number of control points along the perimeter of the building's cross section.However, the effectiveness of this approach can be advanced by adopting a more efficient definition of the shape design variables, such as the reduced basis method, [13] which is capable of managing the changes in the cross-section geometry with a small number of design variables.The reformulation of the shape variables into basis geometries allows the definition of any candidate geometry in the design space as a combination of those bases.This effective method permits a more efficient exploration of the design domain, increasing the accuracy of the aerodynamic response emulation, and lowering the computational burden related to the numerical evaluation of the aerodynamic responses.
Discussions on the computational advantages of the reduced basis method can be found in Rozza et al. [14] and Bachmayr and Cohen, [15] while applications in airfoils or bypass designs are reported in Lassila and Rozza, [16] Skinner et al., [17] and Manzoni et al., [18] among others.
Given the clear advantages of this method with regard to other formulations, its application is gaining momentum in the aerospace engineering field (see, for instance, Skinner et al. [17] ), and it has a great potential to be applied in wind engineering applications, particularly in civil and architectural structures.
In tall buildings design, this approach was adopted in the author's previous work Cid Montoya et al., [19] where a multi-objective surrogatebased optimization framework was developed to identify the Pareto front of optimal cross sections maximizing profitability and minimizing windinduced accelerations.In the aforementioned reference, design constraints were imposed to keep the accelerations caused by some specific wind scenarios, characterized by specific wind directions, below a set of thresholds based on the ISO10137 2007 standard. [20]While some nonlinear features were clearly noted in the response surfaces, the uncommon characteristics on the Pareto front found in that study were mainly driven by the nonlinear design constraints caused by the wind-induced accelerations related to the considered wind scenarios.In the present study, the problem is formulated as a bi-objective optimization problem balancing profitability with the building wind-induced accelerations.The wind scenario selected highlights the nonlinear features of the wind-induced response as a function of the tall building cross section, as well as their effect on the Pareto front properties.
The formulation of the multi-objective surrogate-based shape optimization is reviewed in Section 2, where the formulation of the reduced basis method and the weighted sum and min-max method are discussed.Then, the specific formulation applied in this study to the cross-section optimization of tall buildings is introduced considering the objectives related to along-wind peak acceleration and profitability after a period of several years.Changes in the cross-sectional geometry may decrease the wind drag force, and consequently the along-wind acceleration, but at the expense of reducing the floor area of the building, and consequently its profitability.The optimization problem is solved by adopting a surrogate model whose input is the set of combination factors defining any cross-sectional geometry inside the shape design space, and the output is its aerodynamic response in terms, for instance, of the drag coefficient that is used to evaluate the along wind acceleration.The objectives considered in the optimization problem are described in Section 4. Section 5 provides the details of the application case used in this study, including the parameters of the profit and load models, the basis adopted, the design domain considered, the CFD simulations, and the Kriging surrogate model.Finally, in Section 6, the results obtained for the Pareto optimal set are reported and discussed, focusing on the aerodynamic nonlinear features that drive the complex characteristics of the Pareto front.FORMULATION

| Reduced basis concept and synergies in the surrogate model definition
The framework for shape optimum design of the cross-section of tall buildings presented herein is based upon the application of the so-called reduced basis concept. [13]Initially, a set of different geometries, generally accepted as cross sections for high-rise buildings, but that are not optimal in a mathematical sense, are adopted as the design bases.It is assumed that a certain combination of those generally accepted geometries can perform even better for a certain response.Hence, any modified design is defined as a linear combination of G basis designs: where G is the number of basis vectors and a i is a coefficient establishing the contribution of each basis geometry Y i .Y is the resultant geometry of the cross section of the high-rise building.
Given the multiple symmetry axes usually found in the cross sections of tall buildings, the designs may be advantageously defined in polar coordinates (ρ, θ).In this manner, every point of each cross-sectional design can be obtained following the expression: When the perimeter of the basis geometries is defined by straight lines, the following equation in polar coordinates may be applied to define each basis' segments, given the polar coordinates of two points Figure 1 illustrates with a general example how this method generates new geometries for the cross section of a tall building modifying the corners using two basic geometries as bases: a square and a rhombus.The figure depicts the geometry resulting from applying as combination factors a = (0.5, 0.5), according to Equation (4): Definition of building cross-sectional design from two generic bases A (square) and B (rhombus) using the reduced basis approach and polar coordinates The combination of bases in polar coordinates may lead to a wide range of different geometries.For instance, in the former case, the combination of two bases, whose geometries are defined by straight lines, leads to curved convex segments (see Figure 1 for the differences between the design, in blue, and the reference line, in red).The presence of curved and linear geometries may allow the definition of a wider range of different designs and explore innovative alternatives.
The main advantage of this formulation is that it enables the introduction of relevant shape modifications by adopting a short number of shape design variables.This is of great interest in defining surrogate models for the responses over the design domain as the size of the samples set required to train a surrogate depends heavily in the number of design variables (dimensions of the sampling plan).It must be remarked that the computational burden of the numerical simulations required to obtain the aerodynamic response of each sample is relatively high, which explains the interest in formulating the problem in terms of a small number of design variables capable of setting a relatively broad design domain.
On the contrary, the standard definition of candidate geometries based on a set of control points along the perimeter may require a larger number of shape design variables, with the consequent penalization in the number of samples required in an eventual sampling plan to build a surrogate model.
2.2 | Multi-objective optimization: General definition and methods

| General definition
Multi-objective optimization aims at identifying the values of the design variables that produce the best values of a set of objective functions.The general formulation is given in Equation ( 5), where a, f, g are the vectors of the design variables, objective functions, and design constraints and a i L and a i U are the lower and upper limits of the design variables [21][22][23] : Usually, the objective functions conflict among them, so there are no values of the design variables that produce simultaneously the optimum of each one.Thus, the goal is to obtain a set of vectors a P known as Pareto points that represent designs for which any variation of the variables aimed to enhance the performance of an objective function will worsen, at least, another one.If the aim is to minimize the objective functions, they have the property that it does not exist any point a that accomplishes and The Pareto points can be obtained as the solution of a single objective optimization problem in which the objective function contains all the components of the f(a).Each combination provides an a P set, and the so-called Pareto front is the curve defined by the values of f (a P ).

| Weighted sum method for the identification of the Pareto optimal set
The identification of the optimal Pareto set may be conducted by applying the weighted sum method, presented in expression 7, which is the most common method adopted in this type of problem. [21]It can be seen that the multi-objective formulation is transformed into a single-objective optimization and the values of the weight coefficients represent the relative importance of each objective function considered.By solving the problem for different sets of weights, the Pareto optimal points are obtained.
The main drawback of the weighted sum method is that its effectiveness may be compromised due to two main limitations [24,25] : First, the selection of the weights is conditioned by the problem to be solved and requires expertise to find the most adequate selection of weights to achieve an accurate representation of the complete set of Pareto points.Otherwise, only a part of it is obtained.And the second issue is related to the performance of the method in non-convex regions, [26] since the weighted sum method only provides points in convex regions of the Pareto front.In these cases, more efficient methods are needed, and among the ones available, the min-max method has been used in this research.
2.2.3 | Weighted min-max method for the identification of the Pareto optimal set The weighted min-max method, also known as weighted Tchebycheff method, is formulated as This problem is commonly solved by reformulating Equation (8) as where β is a new design variable that is included in the original vector of design variables a and f j max , f j 0 are the worst and best values of the objec- The main advantage of this method with respect to the alternative weighted sum method in Section 2.2.2 is that it can provide all optimal points in the Pareto set, even in non-convex regions.Nevertheless, it is computationally more expensive.For the application case considered herein, this is not a relevant issue since the problem is formulated as a surrogate-based optimization, decreasing the computational demands of the optimization algorithm.On the contrary, the good performance of this method in non-convex regions is of utmost interest due to the high nonlinear behavior of the aerodynamic responses, as it will be shown in the next sections.

| Formulation of the multi-objective optimization for high-rise buildings
In the application case that is later reported, it is considered of interest the identification of the Pareto optimal set that results from the trade-off between the profit of the planned development and the along-wind response (peak acceleration of the top floor).As will be further discussed in Section 4, the gross internal volume of the building is kept constant, and the height h depends on the area of the cross section of the building.A higher building may provide higher revenues (number of stories M); however, its along-wind response may also be higher, due to its increased flexibility, and the acceleration would also depend strongly on the geometry of the building cross section, giving place to multiple Pareto optima.
It is highlighted that the analytical definition of the basis geometries is made using the line equations in polar coordinates of the segments defining the bases, as discussed in Section 2.1.The modified designs are obtained by applying the pertinent combination factors a i to the analytical line equations of the three basis geometries adopted.Therefore, in the modified geometries, not only the position of the corner is modified, but also the geometry of the curves connecting the corners, obtaining slightly curved segments.This introduces noticeable changes in the aerodynamic responses (see Section 5.4.6)without significantly impacting the floor area of the building, which enables the optimization algorithm to explore efficient design alternatives.
The proposed problem formulation is presented in Equation (10).The objective functions have been normalized introducing the maximum and minimum values of the along-wind peak acceleration ( b € x max , b € x min ) and profit at year Q (P Q max , P Q min ) over the shape design space to guarantee that they are in the interval [0,1]: A being the cross-sectional area of the building for a given value of a, h the height of the building, and V the constant gross volume of the building.
The weighted sum and the weighted min-max methods described above are applied to obtain the Pareto optimal set.In both cases, the optimization problem has been solved using the SQP algorithm. [21,23] Figure 2, a scheme summarizing the process outlined above to conduct the multi-objective shape optimization of the cross section of tall buildings, based upon the application of the reduced basis concept, is provided.

| Surrogate model formulation
Once the geometry of each basis is selected, a design domain is generated where the design variables are the combination factors of those bases.
When three bases are considered, the design domain can be depicted using a ternary plot representation because of Equation ( 1).In the multiobjective optimization problem defined in Equation ( 10), the along-wind response depends on the drag coefficient C D in Equation ( 17), as well as other parameters related to the wind at the particular building location and its structural characteristics.Hence, the drag coefficient over the entire design domain would be required by the optimization algorithm.Similarly, revenues and construction costs over the entire design domain would be required to evaluate the second objective function considered in the optimization problem.It must be remarked that the aerodynamic performance of each candidate design cannot be obtained "on the run" during the optimization due to the intractability of conducting wind tunnel tests or completing CFD simulations for candidate geometries.Because of this, engineers must resource surrogate models [27][28][29] to emulate the aerodynamic responses.
In this work, a Kriging surrogate model is adopted [30,31] to approximate the surface response of the force coefficients over the design space considered for the cross-sectional geometries of the building.The aerodynamic response of the set of samples used to train the surrogate is obtained from 2D URANS simulations.The model may be formulated as Flowchart for the multiobjective shape optimization of the cross section of tall buildings where κ a ð Þ T γ is the trend function and ε a ð Þ is the stationary Gaussian process error model, which permits that the model comprises all the samples used in the training process by correcting the trend function.In this work, a quadratic trend has been adopted, using the implementation included in the DAKOTA framework. [32]

| DESIGN CRITERIA IN TALL-BUILDINGS SHAPE DESIGN
In the following subsections the specific formulation of the two contradictory objectives considered in Equation ( 10) and the methods used for their evaluation are described.

| Profit model for tall buildings
In high-rise developments, the expected costs and returns over time must be carefully considered in order to assess their profitability.In this work, both construction costs and returns models, which are dependent on the building's height and floor area, are considered.In this study, the lot's acquisition cost is not included in the costs model, since it is highly variable and it is not dependent on the building's characteristics.Neither the maintenance costs over the years have been considered.Height, shape, and geometry are the most sensitive cost drivers in tall buildings.In the cost model adopted herein, the so-called Cost Model One based on height in Strelitz [33] is combined with the Cost Model Two based on the plan size, taken from the same reference.The model provides the costs per square meter of gross internal area (GIA) relative to a baseline 15-story building, whose costs are obtained from Watts and Langdon, [34] considering concepts such as the substructure, superstructure, façades, mechanical, electrical, and plumbing services, lifts and escalators, internal walls and finishes, and preliminary works.
For the returns model, the formulation proposed in Tse et al. [35] has been adopted, where rental rates are quantified using a continuous exponential function: In the above equation, RT n k is the annual rent per square meter of GIA for a future year n, at floor k, f is the fixed rate of inflation, r h is a rate adjustment coefficient representing the rent increment with the number of story, and VR n is the assumed vacancy rates and bad debt losses for year n.The total revenue at year Q of a M-story building of uniform floor area A, considering a floor area efficiency E, that is, the ratio of net to gross area, is The profit at year Q of a planned real estate development can be assessed using the next equation, where only construction costs CC are considered: It must be borne in mind that both profit and construction cost are dependent upon the floor area and the height of the building.Therefore, the height of the building will be treated as a non-categorical variable, since it will be assumed a continuous behavior during the optimization process, but its discrete nature will arise in the interpretation of the final results and the final design.In the optimization problem addressed herein, the gross internal volume of the building is imposed as constant.Assuming the same constant inter-story height for all the candidate designs, this means that the total gross floor area is the same for all the considered buildings.In this manner, the Pareto optimal designs would correspond to buildings with the same total GIA, providing a fair comparison among different designs.

| Along-wind service response of tall buildings
The structural performance of the candidate tall buildings inside the considered design domain is assessed based on the acceleration responses for the 5-year return period wind speed, following the formulation in the Eurocode EN 1991-1-4. [36]Further information may be found in the extensive comparative study of various international wind codes in Kwon and Kareem. [37]ccording to Clause (B.4) in Annex B, the standard deviation of the along-wind acceleration at height z for a vertical structure can be evaluated as follows: where C D is the drag coefficient (see Section 4.3.),ρ is the air density, D is the building's cross section width, I v z s ð Þ is the turbulent intensity at the reference height for obtaining the structural factor c s c d , v m z s ð Þ is the mean wind speed at the reference height, R is the resonant response factor, K x is a nondimensional coefficient, m 1,x is the equivalent mass in the along-wind direction, and Φ 1,x z ð Þ is the fundamental mode shape in the along-wind direction (defined according to Clause F.3).This formulation enables estimating the along-wind response from 2D aerodynamic information of the tall building cross section.
The peak acceleration b x is calculated by multiplying the standard deviation of the acceleration by a peak factor defined in equation (B.4) of Clause B.2: ν being the natural frequency of the structure, and T = 600 s.
It must be noted that Eurocode EN 1991-1-4 [36] does not provide expressions to estimate across-wind accelerations for serviceability design.

| Force coefficients definition
The drag, lift, and moment coefficients, as well as the Strouhal number, are defined according to the following formulae: In the former expressions, F D is the mean drag force per span length, positive in the along-wind direction, F L is the mean across-wind force per span length, positive upward, and M is the mean twist moment per unit of span length, positive in the clockwise direction.D is the width of the body, U is the flow reference speed, ρ is the density of the air, and the dominant frequency in the lift force time history is f.The standard deviation of the force coefficients is represented in the following as C D 0 , C L 0 , and C M 0 .

| Fluid mechanics numerical formulation
The unsteady flow around the different geometries adopted as samples for defining the required surrogate model is studied solving the timeaveraged equations for conservation of mass and momentum [38] : where U i is the mean velocity vector, x i is the position vector, t is the time, ρ is the fluid density, u , i is the fluctuating velocity and the over-bar represents the time average, P is the mean pressure, μ is the fluid viscosity, and S ij is the mean strain-rate tensor.From the former equation, the specific Reynolds stress tensor is defined as which is an additional unknown to be modeled based on the Boussinesq assumption for one-and two-equation turbulence models.
where ν T is the kinematic eddy viscosity, S ij is the mean strain-rate tensor, and k is the kinetic energy per unit mass of the turbulent fluctuation.
In this study, the closure problem is solved applying Menter's k-ω SST model for incompressible flows. [39]| APPLICATION CASE: HIGH-RISE BUILDING MULTI-OBJECTIVE SHAPE OPTIMIZATION 5.1 | Definition of the shape design domain for the cross-section of the high-rise building In the frame of the reduced basis shape optimization application case presented herein, three different geometries are adopted as bases: square cylinder (Basis 1), a geometry where the corners of the square cylinder are recessed 0.0535B (Basis 2) and a chamfered cylinder (Basis 3), which are presented in Figure 3a, along with the corresponding combination factors for each basis (a 1 , a 2 , a 3 ).For all the candidate designs, the 50% central length of the lateral sides is not modified.
These three geometries define the shape design space, where modified candidate geometries can be obtained based on a ternary diagram, being the vertices of the triangle the basis geometries.Any point inside the triangle represents the fraction (coefficient a i in Equation 1) of those basis geometries that define a particular candidate geometry.It must be remarked that the modified geometries are a combination of these three bases, whose sides are defined analytically by means of line equations in Cartesian coordinates.Hence, the modified geometries, while showing recessing corners, also feature curved sides as discussed in Section 2.1., which play an important role in the distinctive aerodynamic performance of the candidate designs.This will be further discussed in Section 5.4.
The gross internal volume of the building is assumed to be constant, which results in tall building designs with different heights depending on their cross-sectional area.The reference floor area considered to establish the gross internal volume is a square of 30 (Basis 1, in Figure 3), while the reference height is h ref ¼ 160 m. Figure 3b shows the corresponding designs for the three bases adopted in this work, satisfying the constant gross internal volume condition.

| Definition of the parameters in the profit model
Since the considered application case does not correspond to a particular project, the data used to evaluate returns and construction costs are taken from several references in the literature. [35]The building height can range from 160 to 182.861 m, and the floor area efficiency is assumed to be E = 0.65.The annual rent per square meter on the first floor RT 1 is 680 €/m 2 GIA, the rent increment with the number Basis definition: (a) basis geometries and corresponding contribution factors (equation 1) and (b) constant gross internal volume tall building designs for the three basis cross-sectional geometries of stories (r h ) is set to 0.5%, and the vacancy rate and bad debt losses (VR n ) are expected to be 40% in the first year, 20% in the second year, and 10% in the following years.Ten years (Q = 10) and a fixed rate of inflation f = 1% are considered for the calculations of the returns.The application of the proposed approach for a specific project in a precise location would require a careful assessment of the expected costs and returns and the inclusion of the lot acquisition and maintenance costs.Evidently, adopting different values in the parameters considered in the profit model would cause changes in the Pareto optimal set obtained by solving the multi-objective optimization problem.
The construction costs model applied herein has been tailored based on the data provided for low-medium rise offices in London (Watts and Langdon [34] ) and the trends in cost with height (Cost Model One) and plan size (Cost Model Two) reported in Strelitz. [33]In Table 1, the cost values adopted for buildings with 15, 30, 45, and 60 stories assuming constant floor area are provided.In Figure 4, the curves representing the trends in the construction cost increments for the principal cost drivers are reported for a number of stories in the range (15, 60).The functions modeling the costs per m 2 of GIA depending on the number of stories are provided in Equations ( 21)- (27).It can be appreciated the important impact that height has on the superstructure's cost, which, according to the data in Table 1, was already the main cost driver for low-medium rise office buildings.
The assessment of construction costs may be further refined by considering the changes in the floor area for a given height since smaller floor areas result in a higher proportion of costs per GIA being attributed to the envelope and to the stronger structural system required by slender tall buildings.Based on the data provided in Strelitz [33] for the Cost Model Two, the following expression is used to approximate the increment in construction costs per square meter of GIA with smaller floor areas: where A ref ¼ 900 m 2 m 2 was introduced in Section 5.1.
For any design included inside the shape design domain, its profitability at year Q may be assessed using Equation ( 14), and this information can be used in the multi-objective shape optimization problem in Equation (10).

| Definition of the parameters in the along-wind response model
The along-wind response in serviceability conditions has been evaluated following the procedure in Annex B of EN 1991-1-4, [33] according to Equation (15).In Table 2, the fundamental parameters considered to evaluate the along-wind acceleration at the top of the building for any design inside the shape design domain are summarized.
The along-wind response model enables the evaluation of the peak along-wind acceleration at building height h, which is a function of the drag coefficient, which is itself dependent on the cross-sectional geometry.In the proposed approach, a surrogate model must be defined, whose aim is to provide the aerodynamic data of any candidate design, as it is demanded by the optimization algorithm.

| Design of experiments
The design of experiments has been done in a deterministic way, adopting as criterion a uniform distribution over the shape design space.The three bases along with the nine samples are plotted in Figure 5a over the design domain using a ternary graph in order to show the evenness of the samples' distribution in the shape design space.
A detailed image of the corner region is provided in Figure 5b, highlighting that the corner sides are not straight but curved, as discussed in Section 2.1.Consequently, the perimeters of the considered geometries may cross each other, due to differences in the convexity of the produced designs.In other words, the combination of the three bases produces geometries that are more complex than just the recessing of the corners between Basis 1 and Basis 3.

| Flow modeling and computational approach
The aerodynamic response of the basis geometries and the samples required for the definition of the surrogate model are obtained by means of 2D URANS simulations, adopting Menter's k-ω SST turbulence model.The dimensions of the flow domain are 135 Â 90D.The center of the studied geometries is located at 45D from the inlet and 45D from the upper and lower bounds of the fluid domain.This provides a negligible blockage ratio and prevents the interaction of the boundary conditions with the flow field.All the simulations that are going to be reported were conducted at a Reynolds number, based on the width of the body, of Re ¼ 5 Â 10 4 , which is of the same order of magnitude as the ones in standard aerodynamic wind tunnel tests that are used later for validation.In Figure 6, a scheme of the fluid domain, including boundary conditions, is provided.It must be noticed that the square cylinder is represented out of scale.
The open-source CFD finite volume solver used in this application is OpenFOAM.In the simulations, the time discretization is conducted by means of the two-step second-order backward scheme.Second-order central differencing schemes were also adopted for the divergence and the gradient terms.The PIMPLE algorithm has been chosen for conducting the pressure-velocity coupling in the transient cases considered herein.
The spatial discretization of the computational domain has been done by means of a hexahedral structured mesh, with the exception of the region in the vicinity of the considered geometries where an unstructured mesh, using triangular prisms, has been adopted (buffer zone mesh).
Attached to the body, a high-density structured boundary layer (BL) mesh has been defined (BL mesh).In the BL mesh, the maximum value of the nondimensional height of the first element (y þ ¼ y 1 u Ã ð Þ=ν, where y 1 is the height of the first element attached to the body, u Ã is the friction velocity, and ν is the kinematic viscosity), is y þ max ffi 6, while the mean value is y þ ffi 1:5.This satisfies the requirements of the low Reynolds wall modeling approach that is adopted in all the simulations that are reported next.
In the CFD models of all the samples studied, the corners of the geometries have been modeled as slightly rounded, with a corner radius to width ratio r=D ¼ 0:01, as in Sarwar et al. [40] or Oka and Ishihara. [41]In Figure 7, general and detailed images of the mesh created for conducting the numerical simulations of the aerodynamic behavior of the square cylinder (Basis 1) are presented.
F I G U R E 5 (a) Ternary plot representation of the basis and the samples over the design domain and (b) detail of the geometry of the tall building cross-sectional corner

| Verification and validation
A 1:1 width-to-depth ratio rectangular cylinder has been chosen as reference geometry because it is a conventional cross-sectional shape for high-rise buildings, and it is also a basic geometry extensively studied by means of experimental tests and numerical simulations.Furthermore, the square cylinder has been chosen as one of the basis geometries (Basis 1) in the application case that is being presented.One of the fundamental issues in CFD modeling is to ascertain the sensitivity of the numerical results with respect to the spatial and temporal discretizations.Consequently, the square cylinder geometry has been chosen to conduct verification studies based on the analysis of the changes in important dependent variables such as grid and time step that are systematically refined. [42] Table 3, the fundamental characteristics of the three meshes with different densities, along with the values obtained for the drag coefficient, the standard deviation of drag, lift, and moment coefficients, and the Strouhal number, are reported at a maximum Courant number, Co max ¼ 2.An equivalent simulation conducted imposing a maximum Courant number, Co max ¼ 1, for the medium mesh case has offered nearly identical values for the force coefficients and their standard deviations.In Table 3, it can be noticed that the medium and fine meshes provide very similar results, particularly for the standard deviation of the lift coefficient.Hence, the medium mesh characteristics and a maximum Courant number, Co max ¼ 2, are retained for the subsequent simulations conducted in the frame of this shape optimization research.Another fundamental issue in CFD modeling is to check the accuracy of the numerical simulations by comparing the results with equivalent wind tunnel tests.In Table 4, the numerical results obtained from the 2D URANS simulation conducted for the square cylinder are compared with the experiments at Re ffi 3:5 Â 10 4 by Luo et al. [43] and at Re ¼ 10 5 by Vickery, [44] both taken from Sohankar. [45]The agreement between the numerical results and experimental data can be judged as fairly good, since the values are similar, and a review of the published literature shows a certain scattering in the available wind tunnel data for this geometry.Furthermore, the lower value obtained for the drag coefficient may be attributed to the rounded corners (r/B = 0.01) adopted in the numerical simulations, which is consistent with the experimental results in Carassale et al. [46] conducted in smooth flow at Re ¼ 3:6 Â 10 4 .Finally, slightly higher values in the standard deviation of the lift and drag coefficients should be expected for CFD simulations adopting a two-equation turbulence model, along with the perfect correlation in the spanwise direction, which is an intrinsic feature of 2D simulations.

| Evaluation of the samples' aerodynamic response by 2D URANS simulations
The 12 geometries introduced in Sections 5.1 and 5.4 are studied by means of 2D URANS simulations aiming at obtaining their force coefficients and standard deviations.Hence, for each of these geometries, a finite volume mesh has been built, with similar characteristics to the medium mesh of the square cylinder (Basis 1) introduced in the verification study in Section 5.4.3.Consequently, the mean and maximum values of y þ are similar to the ones reported in that section for the medium mesh case.
Table 5 reports the cells count in the zones in which the flow domain has been subdivided along with the results obtained at 0 angle of attack for the mean drag coefficients and the standard deviations of drag, lift, and moment coefficients.These numerical results are used for the definition of the Kriging surrogate model (see Figure 10) that would enable to obtain the Pareto optimal set in the surrogate-based multi-objective shape optimization problem that is latter addressed.

| Discussion on the fundamental flow features over the shape design domain
Analyzing the fundamental flow features for the 12 studied cross sections, it has been found that the geometry dramatically affects the aerodynamics of each particular shape.These differences are caused by the very different behavior of BL that may be mainly attached along the T A B L E 4 Validation against wind tunnel experiments for the 1:1 rectangular cylinder windward first and second curved inclined planes (Samples 7 and 9), or separated at the corners (Samples 1 and 5, for instance).In all the cases considered, the eventual transition in the shear layers and their interaction with the leeward corners produce also very different flow features, and hence aerodynamic responses (e.g., Samples 2 and 7) .In Figures 8 and 9, images of the flow fields at the instant of maximum lift are provided For each one, a detail of the velocity magnitude field close to the upper windward corners is presented in Figure 8 to understand the effect of the modified sides in the eventual separation of the BL.Furthermore, flow fields of the turbulent kinetic energy are also given in Figure 9, aiming at representing the wake flow features and the interaction of the shear layers with each other and with the leeward geometry of the body.
According to the data in Table 5 and the flow features obtained for the set of samples considered, the aerodynamic behavior of the studied geometries may be classified into three groups: • Basis 1, Basis 2, Sample 1, Sample 2, Sample 3, and Sample 5 are characterized by a massive flow separation at the windward corner and the presence of two shear layers that strongly interact with each other in the proximity of the leeward side of the body.Hence, the standard deviations of the force coefficients show relatively high values, and the drag coefficient is also high.
• Basis 3, Sample 4, Sample 6, and Sample 9 are characterized by shear layers showing a weaker interaction, and consequently, the values of the drag coefficient are lower and the standard deviations of the force coefficients are also low.
• Samples 7 and show distinctive flow patterns, characterized by a narrow shear layer, almost attached to the body.These geometries present intermediate values for the drag coefficient and the standard deviations of the force coefficients, while the Strouhal number is remarkably higher than for the other two groups.According to Kwok, [6] a higher Strouhal number may cause resonant response at lower reduced velocities, with a shorter wind speed return period.

| Surrogate model discussion
Using the DAKOTA framework, [32] a Kriging surrogate model has been trained based on the aerodynamic responses obtained for the three bases and nine samples considered.The surrogate model has three inputs, which are the combination factors (a i , i = 1, 2, 3) of the basis that defines each particular geometry, and one output, which is the drag coefficient required to evaluate the along-wind response of candidate building designs.Hence, the aerodynamic surrogate model A can be formulated as Extended versions of this model can be built considering any of the outputs obtained in the CFD simulations when required in alternative formulations of the problem.
In Figure 10, the surface responses for the drag coefficient, the relative area, and the height, the latest ones evaluated analytically for each geometry, are depicted over the design domain.In the drag coefficient surface response, it can be noticed how the value of the coefficient changes over the entire design domain and that several minima are apparent in the vicinity of one of the boundaries of the shape design domain.
The significance of this circumstance is twofold.Firstly, the differences in the values of the drag coefficient show the sensitivity in the aerodynamic response with the combination of the three basis geometries over the relatively conservative shape design space adopted for this study.
This can also be seen in the second objective function related to the peak along-wind acceleration depicted in Figure 11b.A second issue is the added complexity in the resolution of the optimization problem due to the presence of local minima in one of the objective functions.This explains the distribution of the set of Pareto optimal points that are reported in the next section, where discontinuities in the Pareto front may be identified.A more detailed discussion explains the trends followed by both normalized objective functions over in the design space.
F I G U R E 1 1 Normalized objective functions in Equation ( 9): (a) profit and (b) along-wind peak acceleration 1.Both objective functions decrease in the interval between vertex B3 and the global minimum of the along-wind peak acceleration function, close to Sample 9. Therefore, as they perform similarly, no Pareto points exist in this region of the design space.
2. The trends in the two objective functions conflict in the interval located between the global minimum (close to Sample 9) and the local maximum (close to Samples 7 and 8) of the along-wind peak acceleration.Therefore, a portion of the Pareto front will be obtained in this region.
3. The two objective functions decrease again in the area between the local maximum and the local minimum of the along-wing peak acceleration function (close to Sample 6).Thus, as they perform with the same trend, no Pareto points will exist in this subdomain of the design space.
4. The trends of the objective functions conflict again in the interval between the local minimum of the along-wind peak acceleration and the boundary defined by vertices B1 and B2.Consequently, the remaining points of the Pareto front are obtained by the optimum combinations of both objective functions in this zone of the shape design space.
This complex behavior of the two objective functions prevents the Pareto front from being a continuous line.On the contrary, it will be nonconvex and composed of two disjoint segments.Although this circumstance is quite infrequent, it has been reported in the literature before (see, for instance, Obayashi et al. [47] ), and thus, this research describes an interesting case of multi-objective optimization.
The problem is solved in this study using gradient-based optimization algorithms, which show a good performance in finding the minima as described in the following section.Alternative approaches would involve adopting metaheuristic algorithms (see, for instance, Wahde [48] and Hernández [49] ).6] 6.2 | Non-convex and discontinuous Pareto front The calculation of the Pareto points and the Pareto front has been made using the two formulations presented previously: the weighted sum and the weighted min-max techniques.Figure 12a,b shows the results in the objective function space and in the design variables space, respectively.
The weighted sum method provided only a few points, in red color in Figure 12, which confirms the lack of efficiency of that approach when the geometry of the curve of the Pareto optima is non-convex.
On the other hand, the weighted min-max technique identified the complete Pareto front, as presented in blue color in Figure 12a.It can be observed that it has a discontinuity at a value of the along-wind peak acceleration of 0.0965 m/s 2 (Points A and B in Figure 12).The existence of discontinuous Pareto optima has been described as relatively uncommon, and it is very interesting to find a situation like this in the multiobjective optimization of the cross section of tall buildings in a formulation coherent with the requirements of the real design cases.The discontinuity is due to the existence of the two local minima found in the objective function related to the along-wing acceleration surface (Figure 11b), which is a consequence of the nonlinear nature of the aerodynamic response (Figure 10a).Point B of the Pareto front corresponds to the value of the local minima in the along-wind peak acceleration response of higher value (see Figure 11b) and Point A to another point on that surface with the same value of peak acceleration but not being a minimum.However, the profit value in point B is higher than in point A (see Figure 11a The computational demands of the proposed approach are relatively modest.Reducing the number of design variables to three has allowed the construction of the surrogate model using a limited number of samples.The 2D URANS approach adopted for obtaining the aerodynamic response of this set of geometries has proved to be accurate for the intended use, based on the verification study and the validation with experimental data reported in the literature for the square cylinder.
For the proposed optimization problem, a set of Pareto optimal points has been identified, finding that these geometries are mostly grouped in two regions of the design domain due to the presence of local minima in the along-wind peak acceleration response in the domain.The engineering significance of the Pareto optimal points has been discussed, remarking that in the criterion space, the Pareto front shows two nonconvex regions, which is a relatively uncommon characteristic, that required the use of the weighted min-max method to fully identify the Pareto set.Furthermore, the Pareto front is discontinuous due to the nonlinear nature of the aerodynamic responses, which is a rarity in engineering design.Also, the ability of the proposed approach to be applied in combination with the serviceability requirements established by international standards has been highlighted.These facts confirm the complexity inherent to wind engineering design, further stressing the convenience of introducing nonconventional design techniques such as surrogate modeling, sensitivity analysis, and optimum design that enable a more complete exploration of the design domain.
Most of the optima that form the Pareto front have been identified in regions of the design domain close to Samples 6 and 9, which are the designs with lower drag coefficient due to the low interaction of the shear layers shown in the flow structures.This fact consequently leads to lower values of the along-wind peak accelerations, which is the objective function that models human comfort.On the other hand, the parameters adopted in the profit model suggest that the highest profit is obtained for the design defined as Base 1, which has the maximum cross-sectional area and the minimum height.The Pareto fronts produced by the methodology proposed in this study are a valuable tool for designers to find the most advantageous balance between profit and aeroelastic response.
The good performance of the weighted min-max method combined with the reduced basis method showed that this strategy is a promising approach for more complex problems.Future research will build upon this framework in several directions, such as introducing additional aerodynamic responses, for example, across-wind accelerations, adopting the 3D geometry of the building in the CFD simulations, performing more advanced CFD simulations, considering turbulent properties in the incoming flow, or carrying out further refinements in the surrogate model, among many others.

F I G U R E 6
Computational flow domain and boundary conditions (square cylinder not to scale) F I G U R E 7 Finite volume mesh for the medium density case of the ratio 1:1 rectangular cylinder.(a) General view, (b) buffer zone near the body, and (c) boundary layer mesh and curved corner detail T A B L E 3 Verification study for the square cylinder Mesh Total cells Cells in buffer zone Cells in BL zone

F I G U R E 8
Flow features at the instant of maximum lift force of six geometries defining the surrogate model.Detail of the velocity magnitude field close to the upper windward corner F I G U R E 9 Flow features at the instant of maximum lift force of six geometries defining the surrogate model.Detail of the wake structure based on turbulent kinetic energy field for representative geometries.

1 0
Response surfaces: (a) drag coefficient provided by the surrogate model, (b) relative floor area computed analytically, and (c) tall building height.The design domain is represented by a ternary plot where the axes correspond to the combination factors a i for the three basis geometries.The solid blue squares represent the values obtained for the set of samples considered in the training of the surrogate model.

6 | 6 . 1 |
SOLUTION OF THE MULTI-OBJECTIVE OPTIMIZATION Response surfaces in the design spaceBoth objective functions defined in expression (9) appear in Figure11over the whole design domain herein considered.It can be seen that the profit-related normalized objective function decreases monotonically from the vertex B3 (Basis 3) to the line defined by vertex B1 (Basis 1) and B2 (Basis 2).On the other hand, the normalized objective function associated with the along-wind peak acceleration shown in Figure11bdecreases from vertex B3 reaching a global minimum in the vicinity of Sample 9 (see the colored plane view in the lower part of Figure11b).Then, the along-wind acceleration increases reaching a local maximum in the vicinity of Samples 7 and 8.This normalized objective function decreases again showing a local minimum close to Sample 6. Afterward, the response surface increases monotonically presenting the maxima values along the line linking vertices B1 and B2 in the ternary plot representing the design domain.
), which causes the discontinuity in the Pareto front.F I G U R E 1 2 Graphical representation of (a) Pareto set in the criterion space, (b) selected geometries of Pareto optima, and (c) Pareto set in the design variables space.This paper discussed the influence of nonlinear aerodynamic properties of tall buildings in the Pareto fronts obtained in the frame of multiobjective optimization problems.A bi-objective optimization problem where two competitive objective functions related to the along-wind peak acceleration response and the expected profitability of the building in a 10-year time horizon has been addressed by means of the weighted sum and weighted min-max methods.The implementation of the reduced basis method has enabled the definition of a relatively ample design domain considering only three design variables, which are the combination factors of the three basis geometries considered in the problem.
Parameters for the construction costs model in € per square meter of gross internal area (GIA) for shell and core office buildings for buildings with 15, 20, 45, and 60 stories F I U R E 4 Construction cost % relative to a 15-story low-medium rise office building Geometric definition, mesh characteristics, force coefficients, standard deviations, and Strouhal numbers for the basis geometries and the samples defining the Kriging surrogate model