Shape design of arch dams under load uncertainties with robust optimization

Abstract

Due to an increased need in hydro-electricity, water storage, and flood protection, it is assumed that a series of new dams will be build throughout the world. The focus of this paper is on the non-probabilistic-based design of new arch-type dams by applying means of robust design optimization (RDO). This type of optimization takes into account uncertainties in the loads and in the material properties of the structure. As classical procedures of probabilistic-based optimization under uncertainties, such as RDO and reliability-based design optimization (RBDO), are in general computationally expensive and rely on estimates of the system’s response variance, we will not follow a full-probabilistic approach but work with predefined confidence levels. This leads to a bi-level optimization program where the volume of the dam is optimized under the worst combination of the uncertain parameters. As a result, robust and reliable designs are obtained and the result is independent from any assumptions on stochastic properties of the random variables in the model. The optimization of an arch-type dam is realized here by a robust optimization method under load uncertainty, where hydraulic and thermal loads are considered. The load uncertainty is modeled as an ellipsoidal expression. Comparing with any traditional deterministic optimization method, which only concerns the minimum objective value and offers a solution candidate close to limit-states, the RDO method provides a robust solution against uncertainty. To reduce the computational cost, a ranking strategy and an approximation model are further involved to do a preliminary screening. By this means, the robust design can generate an improved arch dam structure that ensures both safety and serviceability during its lifetime.

Keywords

arch dam / shape optimization / robust optimization / load uncertainty / approximation model

Introduction

Optimization techniques have been widely used in the design of arch-type dams and huge process were achieved with the application of shape optimization techniques [¹–³]. Traditional optimization techniques are realized based on a deterministic approach, which means that the material properties and loading conditions are assumed to be fixed values. As a result, the real-world structures that have been optimized by these approaches suffer from uncertainties that one needs to be aware of. Hence, in any optimization process for arch dams, it is necessary to find a methodology that is capable of considering the influences of uncertainties and generating a solution that is robust enough against the uncertainties.

To consider uncertainties, there exist several methodologies offering diverse solutions strategies. One solution is based on the probabilistic model, which is characterized in such way that the uncertainty is assumed to be known in advanced and that it can be described by appropriate probability distributions [⁴–⁶]). With this in mind, two methods were developed and applied successfully to the optimization under uncertainty: the reliability-based design optimization (RBDO) [⁷], which aims to find the minimum objective function value un-der the probabilistic constraints, generally expressed in terms of estimated failure probabilities; and the robust design optimization (RDO), which aims to minimize both the mean value of the objective function and the variance of the structural responses.

Earlier, researchers added an estimate of the model’s response variance to the objective function, so that they could select the optimal solutions, where variations in the input data led to smaller variations in the output [⁸–¹⁰]. These methods are effective and well-proven, however, since they are based on sample-based estimators of the variance, the computational efforts are high. Further, these methods also rely on a predefined probabilistic model, where again unnecessary uncertainties might be introduced due to the lack of knowledge of the statistical properties of the engaged distributions.

An alternative approach based on non-probabilistic model to this is proposed in Ben-Tal et al. [¹¹,¹²], Sundaresan et al. [¹³], Guo et al. [¹⁴] which is kind of a semi-definite programming to find the robust design candidate. This method, which considers the uncertainty as a parameter within an uncertainty data set U, focuses on searching out an optimized design candidate relatively insensible to uncertainty. We call this approach RO method (Robust Optimization).

In comparison with RO, the RDO and RBDO depend rather visibly on the accuracy of the distribution used to represent the uncertain data, and this character makes results of RDO and RBDO to be less conservative. Additionally, the RO method can find a design candidate against all possible situations generated from an uncertainty set U and makes the solution feasible enough for the general optimization problems with uncertainty. In recent years, there has been a burst in the application of RO for structural optimization problems, for example, Sun et al. [¹⁵]) describes a design in sheet metal forming using robust optimization, Kanno and Takewaki [¹⁶] and Guo et al. [¹⁴] realized a bar truss structure’s optimization under uncertainties with RO methods.

The optimization of arch-type dams is a time consuming process; considering uncertainties, shape optimization becomes an even more complex and computationally expensive challenge. To our best knowledge, the RO is applied to the shape optimization of a dam under uncertainties for the first time in this work. The objective is to minimize the volume and variations in loads (water-level and temperature). For the complex shape optimization process, it is necessary to develop methods to reduce the cost in the calculation of stress analysis. Zhu et al. [³], Schmit and Farshi [¹⁷] offered some techniques to approximate the stress and acquired huge progress in reducing the calculation consumption. In this study, the approximation model and an optimization ranking strategy are introduced to improve the efficiency of the optimization process.

The paper is organized as follows: after the introduction a brief description of the parameterization of the shape optimization process is given. This is followed by the definition of the objective and the definition of a series of constraints which guarantee the save and reliable operation of the dam. For the assessment of the stress states a coupled thermo-mechanical analysis is conducted, where the resulting partial differential equations are finally solved with the Finite Element Analysis in 3D. For further acceleration, a Kriging approximation is used. After- wards, the optimization problem is extended to the RO and linked to the previously introduced FE-model. Finally, an example is provided which compares the results of the RO with the deterministic optimization (DO).

Geometrical model

Description of the shape of an arch dam

The shape of an arch dam is determined by its central vertical section and the horizontal sections at selected elevations. There are different ways to define the shape of the vertical and horizontal sections, which are discussed in literatures [²,³]. In this study, the vertical section is determined by the curve of the upstream side and the thickness of the section; and the horizontal sections are described by central curves and their thickness at selected elevations. The parameters that determine the shapes are called design variables.

The upstream curve of vertical section(seen in Fig. 1(a) varying with the coordinate is expressed as:

(1a)y({{\displaystyle z}}_{\mathrm{1}}=\mathrm{0})=\mathrm{0},

(1b)y=({{\displaystyle z}}_{\mathrm{2}}=\alpha H)=-{{\displaystyle \beta }}_{\mathrm{1}}{{\displaystyle t}}_{\mathrm{3}};

(1c)y=({{\displaystyle z}}_{\mathrm{3}}=\alpha H)=-{{\displaystyle \beta }}_{\mathrm{2}}{{\displaystyle t}}_{\mathrm{3}};

(1d){\displaystyle \frac{\mathrm{d}y}{\mathrm{d}z}}({{\displaystyle z}}_{\mathrm{2}}=\alpha H)=\mathrm{0.}

In the same way, the thickness is:

(2a){{\displaystyle t}}_c({{\displaystyle z}}_{\mathrm{1}}=\mathrm{0})={{\displaystyle t}}_{\mathrm{1}},

(2b){{\displaystyle t}}_c({{\displaystyle z}}_2=\frac{1}{2}H)={{\displaystyle t}}_2,

(2c){{\displaystyle t}}_c({{\displaystyle z}}_{\mathrm{3}}=H)={{\displaystyle t}}_{\mathrm{3}}.

For the horizontal section (see Fig. 1(b)), the he polynomial of central curve is acquired by solving the following equations:

(3)x^{\mathrm{2}}={ay}^{\mathrm{2}}+by,

(4)\tan \varphi ={\displaystyle \frac{dy}{dx}}={\displaystyle \frac{\mathrm{2}x}{\sqrt{\mathrm{4}a{x}^{\mathrm{2}}+b^{\mathrm{2}}}}},

where the parameter \alpha and the central angle \varphi are obtained from the quadratic function of coordinate z:

(5){{\displaystyle a}}_i={\displaystyle {\sum }_{j=\mathrm{0}}^{\mathrm{2}}{{\displaystyle c}}_j}{{\displaystyle z}}_i^j,(i=\mathrm{1},\mathrm{2},\mathrm{3}),

(6){{\displaystyle \varphi }}_i={\displaystyle {\sum }_{j=\mathrm{0}}^{\mathrm{2}}{{\displaystyle \gamma }}_j}{{\displaystyle z}}_i^j,(i=\mathrm{1},\mathrm{2},\mathrm{3}).

The coefficient {{\displaystyle c}}_j,{{\displaystyle \gamma }}_j(j=\mathrm{0},\mathrm{1},\mathrm{2}) can be derived according to the value of {{\displaystyle a}}_i,{{\displaystyle \varphi }}_i,(j=\mathrm{1},\mathrm{2},\mathrm{3}) at control points \left({{\displaystyle z}}_{\mathrm{1}}=\mathrm{0},{{\displaystyle z}}_{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}H,{{\displaystyle z}}_{\mathrm{3}}=H\right) .

Based on the above, the design variables are selected and listed in the following vector:

(7)x={[\alpha ,{{\displaystyle \beta }}_{\mathrm{1}},{{\displaystyle \beta }}_{\mathrm{2}},{{\displaystyle t}}_{\mathrm{1}},{{\displaystyle t}}_{\mathrm{2}},{{\displaystyle t}}_{\mathrm{3}},{{\displaystyle \varphi }}_{\mathrm{1}},{{\displaystyle \varphi }}_{\mathrm{2}},{{\displaystyle \varphi }}_{\mathrm{3}},{{\displaystyle \alpha }}_{\mathrm{1}},{{\displaystyle \alpha }}_{\mathrm{2}},{{\displaystyle \alpha }}_{\mathrm{3}}]}^T\in R^{\mathrm{12}}.

Objective function and constraints

The objective function is considered as the volume of the arch dam, which means the purpose of this optimization is to search out a reasonable minimum dam volume subject to geometry, stress and stability constraints. The general frame of the optimization problem is expressed as:

Minimize:f(x)=volume(x)={\displaystyle \iint \left|{{\displaystyle y}}_u(x,z)-{{\displaystyle y}}_d(x,z)\right|dxdz,}

(8)Subject\hspace{0.17em}to:{{\displaystyle g}}_i(x)\le \mathrm{0}(i=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{...}),

lb\le x\le ub,

where the g_i(x) is the ith constraint, lb,ub are respectively lower and upper boundaries of design variable x, the inequality signs need to be understood in a point-wise manner.

According to Ref. [³], the selected constraints can be described as follows:

1) Stress constraints

For the assurance of safe working condition during lifetime, the maximum stress of an arch dam must be under allowable stress. The stress constraints can be expressed as follows:

(9){{\displaystyle g}}_{\mathrm{1}}={\displaystyle \frac{{{\displaystyle \sigma }}_{\mathrm{1}}}{[{{\displaystyle \sigma }}_{\mathrm{1}}]}}\le \mathrm{1},

(10){{\displaystyle g}}_{\mathrm{2}}={\displaystyle \frac{{{\displaystyle \sigma }}_{\mathrm{3}}}{[{{\displaystyle \sigma }}_{\mathrm{3}}]}}\le \mathrm{1.}

In Eqs. (9) and (10), [{{\displaystyle \sigma }}_{\mathrm{1}}],[{{\displaystyle \sigma }}_{\mathrm{3}}] are respectively the maximum allowable ten-sile stress and the maximum allowable compressive stresses.

2 Geometric constraint

For the convenient construction, the geometrical constraint is generally expressed as the overhang degrees of upstream and downstream, which is represented by s. The principle overhang degree[s] is 0.3 [³] and depending on that, the geometrical constraint could be written as:

(11){{\displaystyle g}}_{\mathrm{3}}={\displaystyle \frac{s}{[s]}}\le \mathrm{1},

where [s] is the maximum allowable overhang degree.

3 Slope stability constraint

The design for an arch dam must fulfill the requirement of the slope’s stability. Consequently, the stability against sliding must be taken into consideration by introducing the coefficients of sliding resistance, which is given by the sliding resistance K_i. The constraint condition can be written as:

(12){{\displaystyle g}}_{\mathrm{4}}={\displaystyle \frac{[{{\displaystyle K}}_i]}{{{\displaystyle K}}_i}}\le \mathrm{1},

where the [K_i ]is the minimum allowable value of sliding resistance for the ith point, and K_i is the coefficient of sliding resistance of ith point.

Finite element (FE) analysis

Thermal analysis

The thermal behavior of an arch type dam is mainly governed by the temperature distribution inside the structure and varying thermal boundary conditions. The transfer of the heat in the arch dam obeys the heat conduction equation, which is for a 3D situation is:

(13){{\displaystyle k}}_{xx}{\displaystyle \frac{{\partial }^{\mathrm{2}}T}{{\partial }^{\mathrm{2}}{x}^{\mathrm{2}}}}+{{\displaystyle k}}_{yy}{\displaystyle \frac{{\partial }^{\mathrm{2}}T}{{\partial }^{\mathrm{2}}{x}^{\mathrm{2}}}}+{{\displaystyle k}}_{zz}{\displaystyle \frac{{\partial }^{\mathrm{2}}T}{{\partial }^{\mathrm{2}}{z}^{\mathrm{2}}}}+q={{\displaystyle \rho }}_s{{\displaystyle C}}_s{\displaystyle \frac{\partial T}{\partial t}},

where k_xx,k_yy,k_zz are respectively thermal conductivities in three directions, {{\displaystyle \rho }}_s is the density of material, C_s is the specific heat. Temperature is denoted with T,x,y,z are spatial coordinates, and t denotes time. The thermal boundary of an arch dam is separated into two parts: the boundary between water and dam surface, treated as the first boundary condition,

(14)T(t,z)={{\displaystyle T}}_w(t,z),

which means the temperature of dam’s surface is equal to the temperature of water T_w.The second part is the boundary between dam surface and air, experiencing a heat transfer process. This situation can be characterized as the second boundary condition:

(15)q=-{{\displaystyle k}}_c\left({\displaystyle \frac{\partial T}{\partial x}}{{\displaystyle l}}_x+{\displaystyle \frac{\partial T}{\partial x}}{{\displaystyle l}}_y+{\displaystyle \frac{\partial T}{\partial x}}{{\displaystyle l}}_z\right),

where k_c is the surface heat transfer coefficient, l_x,l_y,l_z are the cosines of outer normal vector of boundary surface .

Figure 2 offers an example to show the distribution of the temperature in the dam’s body with the above-mentioned thermal boundary. During winter, the temperature of water is higher than the air’s, and this situation is just on the contrary in the summer time. Generally, the thermal strain existing in the dam’s body is proportional to the change in temperature, and the variance of water temperature is smaller than the air temperature [¹⁸].

Static stress analysis

The nodal forces and nodal displacements follow the relations in the FE method expressed as:

(16){{\displaystyle K}}_{g}u={{\displaystyle P}}_a,

where the P_a is the nodal force, u is the nodal displacements, and K_g is the global stiffness matrix. In FE analysis, all types of loads are transferred to the form of equivalent nodal force Pa. In this work, the equivalent nodal force is composed of body force, surface force, and thermal strain. Their element equivalent nodal forces are acquired through following equations:

(17a)\mathrm{b}\mathrm{o}\mathrm{d}\mathrm{y}\hspace{0.17em}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{c}\mathrm{e}:{{\displaystyle {{\displaystyle \mathit{P}}}_{\mathit{a}}}}_{{{\displaystyle q}}_v}^e={\displaystyle \iiint N^{\mathrm{T}}{{\displaystyle \mathit{q}}}_{\mathit{v}}}\mathrm{d}x\mathrm{d}y\mathrm{d}z,

(17b)\mathrm{s}\mathrm{u}\mathrm{r}\mathrm{f}\mathrm{a}\mathrm{c}\mathrm{e}\hspace{0.17em}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{c}\mathrm{e}:{{\displaystyle {{\displaystyle \mathit{P}}}_{\mathit{a}}}}_{{{\displaystyle q}}_s}^e={\displaystyle {\int }_S{N}^{\mathrm{T}}{{\displaystyle \mathit{q}}}_{\mathit{s}}\mathrm{d}S},

(17c)\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{m}\mathrm{a}\mathrm{l}\hspace{0.17em}\mathrm{s}\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}:{{\displaystyle {{\displaystyle \mathit{P}}}_{\mathit{a}}}}_T^e={\displaystyle \iiint {\mathit{B}}^{\mathrm{T}}\mathit{D}{{\displaystyle \mathit{ϵ}}}_{\mathrm{T}}}\mathrm{d}x\mathrm{d}y\mathrm{d}z,

where, the N^T is the shape function matrix, {q_v} is the body force, {q_s} is the surface force working on the boundary S, and {{\displaystyle \epsilon }}_T is the strain caused by the change of temperature.

The strain vector \epsilon is given as:

(18)\epsilon =\left\{\begin{array}{l}{{\displaystyle \epsilon }}_x\\ {{\displaystyle \epsilon }}_y\\ {{\displaystyle \epsilon }}_z\\ {{\displaystyle \gamma }}_{xy}\\ {{\displaystyle \gamma }}_{yz}\\ {{\displaystyle \gamma }}_{zx}\end{array}\right\}=\left\{\begin{array}{l}{\displaystyle \frac{\partial {{\displaystyle u}}_x}{\partial x}}\\ {\displaystyle \frac{\partial {{\displaystyle u}}_y}{\partial y}}\\ {\displaystyle \frac{\partial {{\displaystyle u}}_z}{\partial z}}\\ {\displaystyle \frac{\partial {{\displaystyle u}}_x}{\partial x}}+{\displaystyle \frac{\partial {{\displaystyle u}}_y}{\partial x}}\\ {\displaystyle \frac{\partial {{\displaystyle u}}_y}{\partial z}}+{\displaystyle \frac{\partial {{\displaystyle u}}_z}{\partial y}}\\ {\displaystyle \frac{\partial {{\displaystyle u}}_z}{\partial x}}+{\displaystyle \frac{\partial {{\displaystyle u}}_x}{\partial z}}\end{array}\right\}=Bu.

The stress are calculated according to Hook’s law:

(19)\sigma =D\epsilon =DBu.

Here \epsilon is the strain, \sigma is the stress, the matrix D is defined as the elastic matrix, and B is the relationship matrix between strain and displacements.

Robust shape optimization

Usually, the DO method minimizes the objective function without concerning variations and scatter of parameters in the model. As shown in Fig. 3, the design variable x is assumed to be normally distributed, f(x) is the objective function of variable x. The optimized point P acquired from DO method suffers obvious changes in f(x) caused by small tolerance of x. On the contrast, considering the uncertainty existing in the structural analysis, the robust optimization provides a feasible solution to the uncertainty problem. As is obvious, the point Q obtained by a robust optimization method shows smaller changes of f(x) with the same distribution of the input variable x.

Approximation model

Generally, the optimization procedure of an arch dam in three dimensions is computationally expensive, and the situation becomes even worse when the robust optimization is involved to consider the uncertain factors. Hence, some strategies are introduced to simplify the procedure. The proposed approach of the optimization (shown in Fig. 4) is first transferred into a two stage procedure by introducing a response surface model.

From the long list of available metamodels, also referred to as response surface methods, the Kriging has been chosen due to its known flexibility and robustness w.r.t overfitting [¹⁹], which have conducted a study on the accuracy of various metamodeling techniques under uncertainty.

Kriging is an interpolation method with the ability to provide a best estimate value of a random field [²⁰,²¹]. Its form can be written as:

(20)\widehat{y}={\displaystyle {\sum }_{j=\mathrm{1}}^k{{\displaystyle \omega }}_j}{{\displaystyle f}}_j(x)+Z(\mu (x),R(x,x^{\prime },\theta )),

where {{\displaystyle \omega }}_j is an unknown coefficient, {{\displaystyle f}}_j(x) are the basis functions of the mean value of a stochastic process, and Z(\mu (x),R(x,x^{\prime },\theta )) is assumed to be a stochastic process with zero mean. In this work, Z(\mu (x),R(x,x^{\prime },\theta )) is expressed as a stochastic Gaussian process with zero mean and unit variance. The probability space of Gaussian process is defined based on the correlation function R and correlation length \theta .

The correlation function describes the similarity of observation and new points. Generally, it is characterized in the form of R(x,x^{\prime },\theta ) , where x is the prediction point, x^{\prime } is the observation, and is a vector containing a set of parameters. For simplicity, the correlation length is assumed to be isotropic in this work. The correlation length is the distance within which points are significantly cor- related. Mathematically, is calculated from the following integration formula:

(21)\theta ={\displaystyle {\int }_{-\infty }^{\infty }\rho (\tau )}d\tau =\mathrm{2}{\displaystyle {\int }_{-\infty }^{\infty }\rho (\tau )}d\tau .

However, this formula is not convenient to be realized in real calculation process. Generally, the correlation length is approximately acquired through some optimization strategies based on maximum likelihood estimates and means of cross- validation. In this study, the Maximum Likelihood method is adopted to acquire the correlation length and the results show well enough for the preliminary analysis in the optimization. Here, we used programming UQLab [²²,²³] to acquire the approximation model. Figure 5 and Table 1 give the comparison results of tensile stress and dam volume approximation models.

Formulation of robust optimization under load uncertainty

A general frame of a robust optimization problem is to find a solution, which mostly reduces the cost as well as guaranteeing the safety of the structure even by objecting it to the worst case of uncertainty factors. This case of optimization is also denoted as worst-case design and optimization. As it guarantees reliable and safe performance of the designed structure under uncertainties, it is a deterministic approach to enhance the confidence into the structure. Thus, the notation is chosen here as the confidence robust optimization. The method is based on a pre-defined uncertainty set U, in which the uncertain parameters \xi are assumed to be defined. First, in order to accomplish the requirements of safety and reliability, the RO methodology model is formulated as:

\min imize;f(x,\xi )

(22)subiect\hspace{0.17em}to:f(x,\xi )\le {{\displaystyle f}}_{\mathrm{0}},

{{\displaystyle g}}_i(x,\xi )\le \mathrm{0},i=\mathrm{1},\mathrm{2},\mathrm{3...},m\forall \xi \in U,

where the objective function is the volume of the structure for the chosen application, is the design variables (geometric description of the dam), \xi is a collection of all variables modeling the uncertainty in the system, f₀ is the boundary of the objective function, and g_i are constraints on the safety and performance of the system.

Generally, the uncertainty factors originate from both the variations of material manufacturing and operational conditions, which can be also treated as uncertainty in the structural stiffness and load. However, for an efficient application of the present approach, the uncertainties from the material’s or resistance side are shifted toward the loading side, resulting an increase in the assumed uncertainties in the loads and keeping the material constants fixed. The calculation processes for the load uncertainty and the stiffness uncertainty are different processes. However, for now, the robust optimization is only in the capacity to solve the linear relations. Considering now only the load uncertainty as an ellipsoidal expression can be adopted to capture the range of uncertainty:

(23)U=\{{{\displaystyle P}}_a\in R^n|{({{\displaystyle P}}_a-{{\displaystyle P}}_a^{\mathrm{0}})}^T{{\displaystyle B}}_h({{\displaystyle P}}_a-{{\displaystyle P}}_a^{\mathrm{0}})\le \mathrm{1}\},

where, P_a is the applied load vector, P₀ is the nominal value of P_a , B_h form the shape matrix of the ellipsoid.

The goal is to find a robust design of an arch-type dam, considering both minimum dam’s volume and safe structural performance. The basic single-level programming including the uncertain variables does not fit for this problem anymore; instead, a bi-level programming is involved as follows:

find:x={[{{\displaystyle x}}_{\mathrm{1}},{{\displaystyle x}}_{\mathrm{2}},{{\displaystyle x}}_{\mathrm{3}},\mathrm{...},{{\displaystyle x}}_n]}^T\in R^n

(24)\min imize:f(x):=Vol(x)

subject\hspace{0.17em}to:\underset{\begin{array}{l}{{\displaystyle K}}_g(x)u={{\displaystyle P}}_a\\ {{\displaystyle P}}_a\in U\end{array}}{\max }{{\displaystyle g}}_i(x,{{\displaystyle P}}_a)-{{\displaystyle g}}_i^u\le \mathrm{0},i=\mathrm{1},\mathrm{...},m,

-\underset{\begin{array}{l}{{\displaystyle K}}_g(x)u={{\displaystyle P}}_a\\ {{\displaystyle P}}_a\in U\end{array}}{\max }{{\displaystyle g}}_i(x,{{\displaystyle P}}_a)-{{\displaystyle g}}_i^l\le \mathrm{0},i=\mathrm{1},\mathrm{...},m,

{{\displaystyle x}}_i^l\le {{\displaystyle x}}_i\le {{\displaystyle x}}_i^u,

where

{{\displaystyle g}}_i(x,{{\displaystyle P}}_a)

are the selected structural performance constraints, {{\displaystyle g}}_i^u and are upper and lower boundaries respectively.

Solving this problem equals to solving the lower level program with global optimality. In the lower level program, the optimal value of is acquired under the uncertain load with known resistance leading to a well-defined stiffness matrix K_g(x), which is naturally a function of the current design x .

Let \varsigma ={{\displaystyle P}}_a-{{\displaystyle P}}_a^{\mathrm{0}} , then the \varsigma runs around the nominal value {{\displaystyle P}}_a^{\mathrm{0}} within the boundary of the uncertainty ellipsoid. Thus, for a given design variable x, we find a combination of the vector to increase the structure’s performance. The value of reflects the resistance ability of the design x against the uncertainty. By this, the problem can be rewritten as:

find:\varsigma \in R^n,

(25)\min imizing:-{{\displaystyle g}}_i(x,\varsigma ),

subject\hspace{0.17em}to:{\varsigma }^{\mathrm{T}}{{\displaystyle B}}_h\varsigma -\mathrm{1}\le \mathrm{0}

The selected performance constraint is considered here as the maximum nodal displacement vector. For the structural analysis with FE method, it is known that the ith nodal displacement can be expressed as {{\displaystyle u}}_i={{\displaystyle d}}_i^{T}u , where d_i is a constant vector.

The RO result can now be acquired by solving the following equations:

find:x={[{{\displaystyle x}}_{\mathrm{1}},{{\displaystyle x}}_{\mathrm{2}},{{\displaystyle x}}_{\mathrm{3}},\mathrm{...},{{\displaystyle x}}_n]}^{\mathrm{T}}\in R^n,

(26)\min imizing:f(x),

subject\hspace{0.17em}to:{{\displaystyle d}}_l^T{{\displaystyle K}}_g{{\displaystyle P}}_a^{\mathrm{0}}+\sqrt{h(x)}\le \widehat{u},

In the Eq. (26), the {{\displaystyle d}}_l^T{{\displaystyle K}}_g^{-\mathrm{1}}{{\displaystyle P}}_a^{\mathrm{0}} is the displacement of the node, which owns the maximum displacement under the nominal value of load with design variable The predefined performance in the structure is denoted with \widehat{u} .

The task of solving the lower level problem is that it is usually a non-convex optimization problem and computationally challenging. To solve this problem Guo et al. [¹⁴] provided a confidence formulation, in which the dual problem is used because the dual problem of a non-convex problem with strict global optimality is always a convex problem [²⁴].

According to [¹⁴], h(x) is acquired by solving following semi-definite dual problem (SDP)

find:\delta \in R,\lambda \ge \mathrm{0}

(27)\min imize:-\delta

subject\hspace{0.17em}to\left(\begin{array}{l}-A+\lambda {{\displaystyle B}}_h\\ sym\end{array}\begin{array}{l}-A{{\displaystyle P}}_a^{\mathrm{0}}\\ {{\displaystyle P}}_a^{{\mathrm{0}}^T}A{{\displaystyle P}}_a^{\mathrm{0}}-\lambda -\delta \end{array}\right)

where \lambda is a Lagrange multiplier, the acquired \min \{-\delta \} with the condition of matrix A\ge \mathrm{0} and {{\displaystyle B}}_h\ge \mathrm{0} is smaller than the global value of the primal problem. The matrix A is obtained by calculating the algebraic value of nodal displacements of the node with the maximum displacements:

(28){{\displaystyle u}}_x^{\mathrm{2}}+{{\displaystyle u}}_y^{\mathrm{2}}+{{\displaystyle u}}_z^{\mathrm{2}}={{\displaystyle P}}_a^T({{\displaystyle K}}_g^{-\mathrm{1}}{{\displaystyle d}}_x{{\displaystyle d}}_x^T{{\displaystyle K}}_g^{-\mathrm{1}}+{{\displaystyle K}}_g^{-\mathrm{1}}{{\displaystyle d}}_y{{\displaystyle d}}_y^T{{\displaystyle K}}_g^{-\mathrm{1}}+{{\displaystyle K}}_g^{-\mathrm{1}}{{\displaystyle d}}_z{{\displaystyle d}}_z^T{{\displaystyle K}}_g^{-\mathrm{1}}){{\displaystyle P}}_a={{\displaystyle p}}_a^{T}A{{\displaystyle P}}_a.

In Eq. (28), K_g is the global stiffness matrix (see Eq. (16)), d_i(i=x,y,z) is the vector which makes that {{\displaystyle d}}_i^{T}u={{\displaystyle u}}_i(i=x,y,z), and the {{\displaystyle u}}_i is the nodal displacements on the nodes which owns the maximum displacement in x,y,z directions.

Optimization example of an arch dam

In the sequel the presented methodology shall be applied to the design of a new arch dam with the following assumptions provided only for analysis requirement:

• The height of the dam is 140 m,

• The valley has a shape of ’V’

• The basements on the both sides of valley are assumed to be rigid foundation.

• The average change of the air temperature is assumed as - 5.4°C, the temperature change of the reservoir water is relatively small and assumed to be - 4.1°C

• The water level is assumed as 135 m.

• The lower and upper boundary of design variables selected for the optimization are assumed according to empirical experiences[¹,²,¹⁸] and are shown in Table 2.

The uncertainty of applied load is the equivalent nodal force of the loads combination considered in the optimization procedure:

• self-weight

• hydrostatic pressure

• uplift pressure

• temperature load.

The uncertainty of applied load is expressed as a ellipsoid (see Eq. (23), and the range of the uncertainty is assumed to be two times of the equivalent nodal force under initial loading conditions. The nominal value of loads and material properties are chosen as initial conditions. According to [¹⁸], the nominal value of loads and material properties used in analysis are assumed as given in Table 3:

With the above-mentioned conditions, an arch dam model with the RO method and a model with deterministic optimization method are acquired to do comparison (shown in the Fig. 6). The following Tab. 4 and Tab. 5 demonstrate the results of the obtained geometrical parameters:

In Fig. 6 a visual comparison between the different qualities of the two optimized designs can be seen. As is obvious, the RO leads to more conservative designs.

Since the maximum displacement of arch dam is selected as performance constraint, the Figs. 7 and 8 show maximum summed displacements of these two arch dam models under the same load situations.

Table 6 illustrates the details of two optimal designs. From the direct comparison, it is demonstrated that the model with RO shows better performance under the same conditions, however, the volume is bigger than the model resulted in DO approach.

The question still remains that how designed structures under randomly varying loads behave. Indeed the RO structure behaves more robust, i.e. same variations on the load parameters should lead to changed mean-values and in particular lower variances in the models’ responses. Therefore 200 random load and material property combinations around the nominal values (as showed in Table 4) are generated following a range listed in Table 7 and the corresponding tensile stresses, von Mises stresses, and the sum of the displacements are recorded. The results are visualized in terms of histograms, see Fig. 9.

From Fig. 9, both of the mechanical characters of arch dam realized with RO method show better performances than with the DO methods. The gaps of the tensile stresses and Von Mises stresses are not that clear comparing with the gap of the sum of displacements, especially the von Mises stresses, it shows little difference between these two dam models. The displacements show a clear difference, and the curve’s vibration of RO model is much more smaller than the vibration of DO model. This situation may be caused by the truth that the sum of displacements are selected as the performance constraint.

Conclusions

The task of optimizing an arch-type dam by means of robust optimization has been solved by the combination of Kriging-based metamodeling and a bi-level semi-defined programming approach, which can be regarded as worst-case design optimization. By adopting a 3D finite element model of an arch-type dam under load uncertainties, an non-probabilistic-based robust design is obtained. Comparisons to the results of deterministic optimization show that the RO designs are more robust and will yield more reliable and safe structures. The price to pay is a less effective use of material, still, the RO design has visibly improved the efficient use of material than any classical design would require.

For the further application of the RO method into engineering designs, there still exist some limitations needed to be solved. The format of the constraints should fulfill specific requirements, which means not all the concerned constraints can be transformed into required form. Moreover, the analysis is only suitable for simple linear elastic structures. More efforts are still needed to improve the RO method and make it applicable to wider ranges of engineering problems.

References

[1] Li S, Ding L, Zhao L, Zhou W. Optimization design of arch dam shape with modified complex method. Advances in Engineering Software, 2009, 40(9): 804–808

[2] Seyedpoor S M, Salajegheh J. Adaptive neuro-fuzzy inference system for high speed computing in optimal shape design of arch dams subjected to earthquake loading. Mechanics Based Design of Structures and Machines, 2009, 37(1): 31–59

[3] Zhu B, Rao B, Jia J, Li Y. Shape optimization of arch dam for static and dynamic loads. Journal of Structural Engineering, 1992, 118(11): 2996–3015

[4] Youn B D, Choi K K. A new response surface methodology for reliability-based design optimization. Computers & Structures, 2004, 82(2–3): 241–256

[5] Karadeniz H, Toğan, V, Vrouwenvelder T. An integrated reliability-based design optimization of offshore towers. Reliability Engineering & System Safety, 2009, 94(10): 1510–1516

[6] Sherali H D, Ganesan V. An inverse reliability-based approach for designing under uncertainty with application to robust piston design. Journal of Global Optimization, 2006, 37(1): 47–62

[7] Tan F, Lahmer T. Shape optimization based design of arch-type dams under uncertainties. Engineering Optimization, 2017, 50(9): 1470–1482

[8] Lee K H, Park G J. Robust optimization considering tolerances of design variables. Computers & Structures, 2001, 79(1): 77–86

[9] Lee K H, Eom I S, Park G J, Lee W I. Robust Design for un- constrained optimization problems using the Taguchi method. AIAA Journal, 1996, 34(5): 1059–1063

[10] Sandgren E, Cameron T M. Robust design optimization of structures through consideration of variation. Computers & Structures, 2002, 80(20–21): 1605–1613

[11] Ben-Tal A, Ghaoui L E I, Nemirovski A. Robust Optimization, Princeton and Oxford, MA: Princeton University Press, 2009

[12] Ben-Tal A, Nemirovski A. Robust optimization methodology and applications. Mathematical Programming, 2002, 92(3): 453–480

[13] Sundaresan S, Ishii K, Houser D R. A robust optimization procedure with variations on design variables and constraints. Engineering Optimization, 1995, 24(2): 101–117

[14] Guo X, Bai W, Zhang W, Gao X. Confidence structural robust design and optimization under stiffness and load uncertainties. Computer Methods in Applied Mechanics and Engineering, 2009, 198(41–44): 3378–3399

[15] Sun G, Li G, Gong Z, Cui X, Yang X, Li Q. Multiobjective robust optimization method for drawbead design in sheet metal forming. Materials & Design, 2010, 31(4): 1917–1929

[16] Kanno Y, Takewaki I. Sequential semidefinite program for maximum robustness design of structures under load uncertainty. Journal of Optimization Theory and Applications, 2006, 130(2): 265–287

[17] Schmit L A, Farshi B. Some approximation concepts for structural synthesis. AIAA Journal, 1974, 12(5): 692–699

[18] Zhu B, Gao J, Chen Z, Li Y. Design and Research for Concrete Arch Dams. Beijing: China WaterPower Press, 2002

[19] Jin R, Du X, Chen W. The use of metamodeling techniques for optimization under uncertainty. Structural and Multidisciplinary Optimization, 2003, 25(2): 99–116

[20] Fenton G A, Griffiths D V. Risk Assessment in Geotechnical En-gineering. Hoboken: John Wiley and Sons, 2008

[21] Rasmussen C E, Williams C K I. Gaussian Processes for Machine Learning. London: MIT Press, 2006

[22] Marelli S, Lataniotis C, Sudret B. UQLab User manual- Kriging(Gaussian Process model), Report UqLab Chair of Risk, Safty and Uncertainty Quantification, ETH Zurich, 2015, 9–105

[23] Marelli S, Sudret B. UQLab: A framework for uncertainty quan-tification in Matlab. In: Proceedings of 2nd International Conf-erence on Vulnerability, Risk Analysis and Management (ICVRAM2014). Liverpool: 2014, 2554–2563

[24] Boyd S, Vandenberghe L. Convex Optimization. Cambridge: Cambridge University Press, 2004

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag GmbH Germany, part of Springer Nature