Risk-based probabilistic thermal-stress analysis of concrete arch dams

Abstract

The probabilistic risk of arch dam failure under thermal loading is studied. The incorporated uncertainties, which are defined as random variables, are associated with the most affecting structural (material) properties of concrete and thermal loading conditions. Karaj arch dam is selected as case study. The dam is numerically modeled along with its foundation in three-dimensional space; the temperature and thermal stress distribution is investigated during the operating phase. The deterministic thermal finite element analysis of the dam is combined with the structural reliability methods in order to obtain thermal response predictions, and estimate the probability of failure in the risk analysis context. The tensile overstressing failure mode is considered for the reliability analysis. The thermal loading includes ambient air and reservoir temperature variations. The effect of solar radiation is considered by an increase in the ambient temperatures. Three reliability methods are employed: the first-order second-moment method, the first-order reliability method, and the Monte-Carlo simulation with Latin Hypercube sampling. The estimated failure probabilities are discussed and the sensitivity of random variables is investigated. Although most of the studies in this line of research are used only for academic purposes, the results of this investigation can be used for both academic and engineering purposes.

Keywords

arch dams / probabilistic analysis / thermal stress / sensitivity / reliability

Introduction

Concrete arch dams are important and critical infrastructures. Their safe operation is of great concern because their damage or failure can lead to economic losses and potential loss of life. Thermal loading is one of the main loadings of thin arch dams which should be considered in their safety analysis and design [¹,²]. Temperature variation and the associated thermal stresses should be determined to define initial conditions in dam safety analysis. Temperature effects on arch dams can be studied in two distinct phases: construction and operation [³]. After the construction phase, the closure or grouting temperature is considered as free-stress temperature for the next phase. During the operation phase, thermal loadings due to variations of ambient air temperature, reservoir temperature, and solar radiation cause difference between dam’s temperature and the closure temperature. Because in this phase the dam is assumed to be a continuous system, the temperature difference would result in thermal stresses. However, for proper evaluation of the stresses, the thermal analysis of arch dams should be done in three-dimensional space [³].

The arch dams are frequently investigated using deterministic thermal analysis methods, among them [³–⁹], but because of the uncertainties related to the parameters describing the structural system and the loading conditions, probabilistic risk assessment methods are preferred. The probabilistic methods provide more comprehensive information regarding the risk levels revealed by a specific structure. Moreover, these methods give insights and perspectives that deterministic analysis cannot procure [¹⁰–¹²]. In addition to probabilistic approaches, utilizing the stochastic analysis is inevitable in complex systems to comprehensively model their behaviors. For instance, the mechanical properties of nanocomposites is influenced by miscellaneous uncertainty; hence, the stochastic impacts of the involved randomness have to be considered properly and thoroughly [¹³–¹⁵]. The different types of uncertainties are classified as aleatory (inherent randomness or natural variability) and epistemic (lack of knowledge) [¹⁶]. While aleatory randomness cannot be reduced, epistemic uncertainty can be reduced by increasing the amount of information. Comprehensive assessment of the risk of thermal-induced structural failure requires a robust numerical method that accurately models the behavior, and explicit consideration of important sources of uncertainty. A primary source of modeling uncertainties lies in the definition of numerical model parameters, specifically those related to the structural system and the loading conditions [¹⁷,¹⁸]. The uncertain parameters are defined as random variables with proper distribution functions and assigned values.

The risk in dam safety analysis is affected by the combined impact of failure scenario, probability of occurrence and the associated consequences [¹⁶,¹⁹–²¹]. Therefore, the risk analysis of dams first requires identification of potential failure scenarios, and then quantification of the conditional probability of these scenarios for different loading conditions. The failure scenario can be presented as limit-state inequality given by

(1)\Psi \left(\mathit{X}\right)\le \mathrm{0},

where Y(X) is the limit-state (performance) function with the failure condition defined by the above inequality, and X is the vector of random variables in the problem. Y(X) = 0 defines the failure or limit-state surface. The probability of occurrence of the failure scenario, P_f [Y(X)≤0], which is equal to the exceedance probability of the corresponding limit-state function, can be computed through the probabilistic framework [¹⁰]

(2)P_{\mathrm{f}}\left(\Psi \right(\mathit{X})\le \mathrm{0})={\displaystyle {\int }_{\Psi \left(\mathit{X}\right)\le \mathrm{0}}\mathrm{\Pi }\left(\mathit{X}\right)\mathrm{d}\mathit{X}},

where P(X) is the joint probability density function of the random variables. There are various statistical techniques to quantify the risk. These techniques, which are known as structural reliability methods, use different mathematical formulations. Assuming that the random variables are normally distributed, the reliability index b, is then defined as

(3)\beta =-{\mathrm{\Phi }}^{-\mathrm{1}}\left(P_{\mathrm{f}}\right),

where F is the standard normal cumulative distribution function.

The present study deals with the probabilistic assessment of the arch dam’s failure risk under thermal loading which incorporates the uncertainties associated with the structural (material) properties and the loading conditions. Karaj arch dam, located in Iran, is selected as case study. The dam is numerically modeled along with its foundation in three-dimensional space to investigate its temperature and thermal stress distribution during the operating phase. The deterministic thermal finite element analysis of the dam is combined with the structural reliability methods in order to obtain better thermal response predictions, and estimate the probability of failure in the risk analysis context. A single overstressing failure mode is defined for the reliability analysis. The thermal loading includes ambient air and reservoir temperature variations. The effect of solar radiation is considered by increasing the ambient temperatures. In the risk analysis, sources of uncertainty related to the thermal loading and material properties are treated using three reliability methods: 1) first-order second-moment method (FOSM), 2) first-order reliability method (FORM), and 3) Monte-Carlo simulation with Latin Hypercube sampling (MC-LHS).

Uncertainty treatment

The temperature distribution and consequently the thermal stresses in arch dams are affected by different parameters. These parameters can be classified either as determined variables, or as random variables whose values are not exactly known and, therefore, are reliant on uncertainty. For the variables considered as random, it is necessary to determine their probability distribution functions, for which, results of research, experiments, and field tests should be taken into account [¹⁶]. The uncertainties must reflect change in environmental loads and the actual condition of the structure. The two main sources of uncertainty captured in this study are related to the thermal (environmental) actions and the material (structural) parameters.

A variety of approaches have been used to study the effects of modeling uncertainties on the structural response. In this research, for propagating modeling uncertainties and quantifying their effect on the structural response, the following methods are implemented: FOSM, FORM, and MC-LHS. These methods have different accuracy and computational cost.

The FOSM has the least computational cost, so it is first applied in order to obtain an initial estimation of the occurrence probability of failure scenario, and to determine the contribution of each random variable to the variance of the limit-state function, Y(X) [¹⁶]. If the mean values of the random variables are defined through vector M, in the FOSM, Y(X) is linearized using a Taylor series expansion about the mean point. The mean value of Y is determined by m_Y = Y(M), and the standard-deviation of the limit-state function σ_Y, is computed from the gradients of Y(X). The linear approximation may be problematic when the limit-state function is highly nonlinear. The reliability index, b, is computed by dividing m_Y to σ_Y. A sensitivity importance vector λ, can be also calculated by the FOSM whose components measure the relative importance of random variables [²²]. The greater absolute value of λ_i shows the greater importance of the corresponding ith random variable. The sign of λ_i indicates whether the random variable is a resistance or load variable, i.e., whether an increase in the realization of the random variable yields a larger or smaller reliability index, respectively.

The FORM is another class of reliability-based methods that uses linear approximation of the failure or limit-state surface, Y(X) = 0. The approximation is centered around a design point which is the point on the failure surface associated with the highest probability of failure. The distance between the design point and the origin is defined as the reliability index. This method is very effective at handling large numbers of random variables, and the approximation is very good at low failure probabilities. Its detailed formulation can be found elsewhere [¹⁷,²²]. The gradient computations required in the FORM would result in the sensitivity measures. The importance vector, called as α-vector, is related to the sensitivity measures and reveals the relative importance of the random variables. The contribution of each random variable to the total variance of the limit-state function is represented by the α-vector. The sign of the components of the α-vector indicates that the random variable is load or resistance if the sign is positive or negative, respectively. There are two other importance vectors called as δ-vector and η-vector. The δ-vector and the η-vector represent the influence of the mean and the standard deviation of the random variables on the reliability index, respectively [²²].

The benchmark and most reliable method is sampling, such as Monte-Carlo method, which is used for its ease of application in problems formulated in terms of a performance function [²³]. In the Monte-Carlo method realizations of each random variable are generated, which form a simulation model, and then the model is analyzed to determine the thermal response. By repeating the process for thousands of sets of realizations, distribution of thermal response results associated with the input random variables is determined [¹⁷]. The exceedance probability of failure scenario, P_f, is estimated by dividing the number of simulations where failure scenario occurred to the total number of simulations. While conceptually straightforward, the Monte-Carlo procedure can become computationally very intensive because the number of simulations performed should be large enough to capture the searched probability. The remedy can be using more efficient sampling methods, among them Latin Hypercube Sampling (LHS) is selected for this research because of its efficiency and simplicity. The LHS is a variance reduction sampling method that stratifies variable marginal distributions in order to fully cover the range of each variable in a more efficient way than pure Monte-Carlo sampling [²⁴]. Unfortunately, the appropriate sample size N, cannot be formerly determined to achieve a certain confidence level. However, using a relatively high N that is substantially larger than the number of random variables will result in reasonably accurate estimates for practical purposes.

Thermal-stress analysis of arch dams

System and boundary conditions

The system considered here is a concrete arch dam which is affected by environmental thermal actions. The governing equation of conductive heat transfer within body of a homogenous and isotropic arch dam in a three-dimensional Cartesian system is defined by [²⁵]

(4)k{\nabla }^{\mathrm{2}}T=\rho C_p\dot{T}-Q,

in which k is the isotropic thermal conductivity coefficient (W/(m²·°C)), T is the temperature (°C), {\nabla }^{\mathrm{2}} is Laplacian operator, r is the concrete density (kg/m³), C_p is the specific heat (J/(kg·°C)), Q is the internal heat generation per unit volume such as hydration process (W/m³), and ‘dot’ represents derivative with respect to time. The boundary conditions of the above equation are

(5a)T={{\displaystyle T}}_{\mathrm{w}}\mathrm{o}\mathrm{n}{{\displaystyle \mathrm{\Lambda }}}_{\mathrm{w}},

(5b)k{\displaystyle \frac{\partial T}{\partial n}}=q_{\mathrm{a}}-q_{\mathrm{c}}-q_{\mathrm{r}}\mathrm{o}\mathrm{n}{{\displaystyle \mathrm{\Lambda }}}_{\mathrm{a}}\mathrm{.}

The first boundary is a prescribed temperature on boundary L_w. In problem of arch dams, this boundary is considered at the dam-water interface, and T_w is the water temperature. On the exposure surface or dam-air interface L_a, the heat can be exchanged through the convective heat flux q_c, the radiative heat flux q_r, or the solar radiation flux q_a; all fluxes are in W/m². n represents the outflow normal direction of L_a. The convective heat flux is the result of temperature difference between L_a and ambient air temperature. It is formulated by Newton’s cooling law [⁶]

(6){{\displaystyle q}}_c=h\left(T-{{\displaystyle T}}_{\mathrm{A}}\right),

where h is the convection coefficient (W/(m²·°C)), T is the temperature of dam surface (°C), and T_A is the ambient air temperature (°C). The heat transfer by convection is a complex phenomenon and it can be influenced by many parameters such as shape and rugosity of the surface, viscosity, and speed of wind. However, some formulas have been proposed which linearly relate h to the wind speed [³].

The heat can be also exchanged by electromagnetic radiation from the exposure surfaces of the dam. The radiative heat flux, q_r, can be measured by the Stefan-Boltzman law

(7)q_{\mathrm{r}}=e{C}_{\mathrm{s}}\left(T^{\mathrm{4}}-T_{\mathrm{A}}^{\mathrm{4}}\right),

where e is the emissivity of surface, and C_s is the Stefan-Boltzman constant (C_s = 5.669×10⁻⁸ W/m²).

The solar radiation heat flux absorbed by the body of the dam is dependent on the total amount of solar energy reaching the surface and solar absorptivity of the surface [³]. The amount of solar radiation reaching the concrete surface follows cyclical seasonal variation. Unlike gravity dams, in concrete arch dams, the total solar radiation reaches have spatial variation through the whole exposed surfaces of arch dam [³]. To simplify the model, in this research, this effect is considered by increasing the air, and consequently, the water temperatures.

The daily air temperature variation, T_air, can be approximated by a sinusoidal function [⁶]

(8)T_{\mathrm{a}\mathrm{i}\mathrm{r}}={\overline{T}}_{\mathrm{a}\mathrm{i}\mathrm{r}}+A_{\mathrm{a}\mathrm{i}\mathrm{r}}\cos \left[{\displaystyle \frac{\mathrm{2}\mathrm{\pi }t}{\mathrm{365}}}\right],

where A_air is the amplitude (°C), {\overline{T}}_{\mathrm{a}\mathrm{i}\mathrm{r}} is the mean annual temperature (°C), and t is the time (days). The amplitude A_air, is dependent on the maximum and minimum average monthly temperatures and the yearly mean average temperature. The reservoir temperature is estimated at various depths of the reservoir. In this research, the Bofang’s experimental-analytical formula is used to predict water temperature stratification in the deep reservoir.

Finite element formulation

By descretizing the temperature domain as T = N^TT, where N^T is the matrix of shape functions and T is the vector of unknown temperatures, the finite element formulation of Eq. (4) considering the defined boundary conditions of Eq. (5), would be in the form of

(9)\mathit{H}\dot{\mathit{T}}\left(t\right)+\mathit{K}\mathit{T}\left(t\right)=\mathit{Q}\left(t\right),

where H is the system heat capacity matrix, K is the system thermal conductivity matrix, and Q(t) is the system heat flux vector which is the sum of fluxes due to internal heat generation, solar, convection, and radiation using the following equations:

(10)\mathit{H}=\underset{\mathrm{\Omega }}{\overset{}{{\displaystyle \int }}}\rho C_p\mathit{N}{\mathit{N}}^{\mathrm{T}}\mathrm{d}\mathrm{\Omega },

(11)\mathit{K}=\underset{\Omega }{\overset{}{{\displaystyle \int }}}k\left({\displaystyle \frac{\partial \mathit{N}\partial {\mathit{N}}^{\mathrm{T}}}{\partial x^{\mathrm{2}}}}+{\displaystyle \frac{\partial \mathit{N}\partial {\mathit{N}}^{\mathrm{T}}}{\partial y^{\mathrm{2}}}}+{\displaystyle \frac{\partial \mathit{N}\partial {\mathit{N}}^{\mathrm{T}}}{\partial z^{\mathrm{2}}}}\right)\mathrm{d}\Omega +\underset{{\mathrm{\Omega }}_s}{\overset{}{{\displaystyle \int }}}\mathit{N}h\mathrm{d}s+\underset{{\mathrm{\Omega }}_s}{\overset{}{{\displaystyle \int }}}\mathit{N}{h}_r\mathrm{d}s,

(12)\mathit{Q}\left(t\right)=\underset{\mathrm{\Omega }}{\overset{}{{\displaystyle \int }}}\mathit{N}Q\mathrm{d}\Omega +\underset{{\mathrm{\Omega }}_s}{\overset{}{{\displaystyle \int }}}\mathit{N}{q}_{\mathrm{a}}\mathrm{d}s+\underset{{\mathrm{\Omega }}_s}{\overset{}{{\displaystyle \int }}}\mathit{N}h{T}_{\mathrm{A}}\mathrm{d}s+\underset{{\mathrm{\Omega }}_s}{\overset{}{{\displaystyle \int }}}\mathit{N}{h}_{\mathrm{r}}{T}_{\mathrm{A}}\mathrm{d}s,

where Ω and Ω_s are the total domain of the system and its surface, respectively.

In the presence of radiation, K is dependent on the unknown temperatures. So the above equations are a system of nonlinear ordinary differential equations. The solution at each time step may be carried out by iteration [²⁵]. Other sources of nonlinearity in the finite element implementation of thermal conduction within the dam body are: (1) the time or temperature dependence on mechanical-thermal properties of concrete, (2) the presence of latent heat effects during phase change, and (3) the presence of anisotropic conductivities. In the operation phase, the heat transfer process will not be affected by the latent heat effects. Since the temperature varies in a relative narrow region of −18°C–36°C during this phase, the temperature dependence of thermal and mechanical properties of concrete are negligible [³]. Therefore, the only source of nonlinearity in Eq. (9) is the presence of radiative boundary conditions.

Thermal stresses

The difference between calculated nodal temperatures T and the nodal closure (grouting) temperatures T₀ for each element of the dam domain would result in a thermal strain vector as

(13){{\displaystyle \mathit{ϵ}}}_{th}=\alpha \left(\mathit{T}-{\mathit{T}}_{\mathrm{0}}\right),

where α is the coefficient of thermal expansion (1/°C). The thermal strain generates the thermal force vector in the following form

(14){\mathit{F}}_t={\displaystyle {\int }_{\mathrm{\Omega }}{\mathit{B}}^{\mathrm{T}}\mathit{D}}{\mathit{ϵ}}_{th}\mathrm{d}\mathrm{\Omega },

in which B is the strain-displacement matrix, and D is the stress-stain matrix. This force vector contributes in the force balance of the system and leads to nodal displacement vector δ and element stresses through the following equation:

(15){\mathit{F}}_t={\mathit{K}}_s\mathit{\delta },

(16){\mathit{K}}_s=\underset{\mathrm{\Omega }}{\overset{}{{\displaystyle \int }}}{\mathit{B}}^{\mathrm{T}}\mathit{D}\mathit{B}\mathrm{d}\mathrm{\Omega }\mathrm{.}

Under long-term action of stress, creep phenomenon would result in an increase in the concrete strain, so relieves some of the induced static thermal stresses. The stress-creep-temperature interaction analysis is beyond the scope of this study and is not considered.

Karaj arch dam model

Karaj dam, which is located on Karaj River, was the first multi-purpose dam in Iran completed in 1961 (Fig. 1(a)). It is 180m high double curvature arch dam with the crest length of 384m. The dam width is about 32m at the base and 8m at the crest level. The complete description of its geometry can be found in Ref. [³]. The dam is modeled along with a part of its rock foundation using the finite element method, the finite element mesh is shown in Fig. 1(b). The only source of loading is the thermal effects during its operation phase; the other static loads due to self-weight or hydrostatic pressure is not considered. This is because of the fact that tensile stress distribution due to temperature, self-weight, and hydrostatic load is very similar with that of only temperature load conditions. So thermal loading has a significant effect on tensile stress response of this concrete arch dam [³].

Because of small thickness of the arch dam, it can be assumed that the majority of the hydration heat has been dissipated in a few years after the dam completion [³]. Therefore, the generation of internal heat by hydration of concrete is not considered. The mean monthly reservoir depth of Karaj dam shows 40m difference between winter and summer [³]. But in this study, the reservoir depth is considered constant same as the dam height. At the dam-water interface, neither convection nor radiation occurs, and the concrete temperature is considered equal to the water temperature [⁶]. At the concrete-air interface, the considered boundary conditions are convection between concrete and air, and radiation from concrete surface. The effects of solar radiation are taken into account as an increase in the ambient air and water temperatures. The foundation thermal effects are ignored. The closure temperature recorded at the time of grouting of the dam monoliths was different along the dam height, but the same closure temperature is considered for the entire dam body.

The concrete and foundation rock are assumed to be homogeneous and linear elastic. The effective mechanical and thermal properties of concrete are assumed as random variables, but isotropic and temperature-independent. The effects of uncertain random variables on the predicted limit-state function are incorporated employing the FOSM, FORM, and MC-LHS methods in the probabilistic risk-based thermal analysis of the dam. The random variables and the limit-state function are defined in the next sections. In the MC-LHS method, there is no predefined sample size N to achieve a certain confidence level [²⁶]. Some formulas have been presented for various applications [²⁷]; one of the simplest formulas suggests that [²⁸]

(17)N>-{\displaystyle \frac{\ln \left[\mathrm{1}-C_l\right]}{P_{\mathrm{f}}}},

where C_l is the confidence level, and P_f is the exceedance probability of failure. Based on this convenient formula, about 4000 samples are required for a 98% confidence level and P_f = 10⁻³. This P_f value is reasonable in dam engineering [²⁹]. In this study, total of 1000, 2000, and 4000 realizations of the dam system are separately generated for the MC-LHS method. These numbers are relatively high such that allow acceptable level of accuracy and lead to statistically significant results. In all above-mentioned procedures, the impact of modeling uncertainties on the thermal response may be significantly affected by the correlations between the input random variables [³⁰]. There is insufficient data to quantify these correlations, so values are typically based on previous experiments and expert judgment.

Random variables

The assessment of modeling uncertainties focuses on those determine the structural material properties and the loading (boundary) conditions or environmental parameters which are defined as random variables. The uncertainties of concrete material properties show the extent to which the idealized model accurately represents real behavior. The thermal and mechanical properties of concrete which are considered as random variable are listed in Table 1 along with their probability density function, mean and standard deviation. It should be noted that some of these variables have large uncertainties due to the lack of data (epistemic). The large uncertainty justifies using the probability distributions such as uniform or triangular with large coefficients of variation. This is the case for thermal conductivity, emissivity, coefficient of thermal expansion, and specific heat. For the other random variables, some information is available from dam sites or empirical data of similar dams, so they are assumed to be lognormally distributed, where the mean and standard deviation are obtained from previous research [³,⁶,²⁵,²⁹–³³]. The material’s random variables are all assumed to be uncorrelated.

The closure temperature or the concrete temperature at time zero for transient computations of Eq. (9) is usually determined by a steady-state heat transfer analysis with applying mean annual air and water temperatures. But in practice, it is common to relate the closure temperature T₀ of the dam with the mean annual air temperature {\overline{T}}_{\mathrm{a}\mathrm{i}\mathrm{r}}. In this study, the following correlation is considered [³⁰]

(18)T_{\mathrm{0}}={\overline{T}}_{\mathrm{a}\mathrm{i}\mathrm{r}}-{\mathrm{3}}^{°}\mathrm{C}\mathrm{.}

Environmental uncertainties in the thermal loading are related to the variation of ambient air and water temperatures. The ambient air temperature can be predicted using Eq. (8). In this equation, the amplitude and the mean annual air temperature are considered as independent random variables listed in Table 2.

The effect of solar radiation is considered as an increase in the ambient air temperature [⁶]

(19)\mathrm{\Delta }{T}_{\mathrm{s}\mathrm{o}\mathrm{l}\mathrm{a}\mathrm{r}}\left(t\right)={\mathrm{B}}_{\mathrm{s}}+{\mathrm{A}}_{\mathrm{s}}\cos \left({\displaystyle \frac{\mathrm{2}\mathrm{\pi }t}{\mathrm{365}}}\right)=\mathrm{3.5}-\mathrm{1.5}\cos \left({\displaystyle \frac{\mathrm{2}\mathrm{\pi }t}{\mathrm{365}}}\right),

where A_s and B_s are two constants that are assumed as 3.5 and −1.5, respectively.

So the applied air temperature would be

(20)T_{\mathrm{a}\mathrm{i}\mathrm{r}+\mathrm{s}\mathrm{o}\mathrm{l}\mathrm{a}\mathrm{r}}\left(t\right)=\left({\overline{T}}_{\mathrm{a}\mathrm{i}\mathrm{r}}+\mathrm{3.5}\right)+\left(A_{\mathrm{a}\mathrm{i}\mathrm{r}}-\mathrm{1.5}\right)\cos \left({\displaystyle \frac{\mathrm{2}\mathrm{\pi }t}{\mathrm{365}}}\right)\mathrm{.}

The reservoir water temperature T_water, which is function of the depth of water y and time t, is correlated to the air temperature using the Bofang’s method with some empirical modifications [³,²⁵].

(21)T_{\mathrm{w}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{r}}(y,t)={\overline{T}}_{\mathrm{w}}\left(y\right)+A_{\mathrm{w}}\left(y\right)\cos \left[{\displaystyle \frac{\mathrm{2}\mathrm{\pi }(t-\zeta )}{\mathrm{365}}}\right],

(22){\overline{T}}_{\mathrm{w}}\left(y\right)=c+\left({\overline{T}}_{\mathrm{s}\mathrm{u}\mathrm{r}\mathrm{f}}-c\right){\mathrm{e}}^{\left(-\mathrm{0.04}y\right)},

(23)A_{\mathrm{w}}\left(y\right)=A_{\mathrm{s}\mathrm{u}\mathrm{r}\mathrm{f}}{\mathrm{e}}^{\left(-\mathrm{0.018}y\right)},

(24){\overline{T}}_{\mathrm{s}\mathrm{u}\mathrm{r}\mathrm{f}}=\mathrm{3.84}+\mathrm{0.76}{\overline{T}}_{\mathrm{a}\mathrm{i}\mathrm{r}},

(25)A_{\mathrm{s}\mathrm{u}\mathrm{r}\mathrm{f}}=\mathrm{3.84}+\mathrm{0.76}{A}_{\mathrm{a}\mathrm{i}\mathrm{r}},

(26)c={\displaystyle \frac{T_{\mathrm{b}\mathrm{o}\mathrm{t}}-{\overline{T}}_{\mathrm{s}\mathrm{u}\mathrm{r}\mathrm{f}}{\mathrm{e}}^{\left(-\mathrm{0.04}H\right)}}{\mathrm{1}-{\mathrm{e}}^{\left(-\mathrm{0.04}H\right)}}},

(27)\zeta =\mathrm{30}\times \left(\mathrm{2.15}-\mathrm{1.3}{\mathrm{e}}^{-\mathrm{0.085}}\right)=\mathrm{28.68},

where T_w (y) and A_w (y) are the mean annual and amplitude of the water temperature at depth y, respectively, z is the phase lag between air and water temperatures, T_surf and A_surf are the mean annual and amplitude of the water surface temperature, respectively, T_bot is the temperature of the deep reservoir which is assumed constant as 6°C, and H is the depth of reservoir. To be applied in numerical model, the mean annual and amplitude of the water temperature are computed and applied separately at nine different levels below the water surface. Therefore, the number of thermal loading random variables is 23, and the total number of random variables in the model is 30.

Limit-state function

The proper definition of the limit-state function and the associated failure scenario is a key issue to ensure the soundness of the solutions obtained with the probabilistic risk analysis [¹⁶]. This study simply focuses on overstressing limit-state as it is one the most important failure scenario when evaluating the static performance of arch dams [¹⁰]. This overstressing is assessed in tension, because the previous analyses of Karaj arch dam have shown that the compressive stresses are far below the compressive strength of concrete [³]. The tensile overstressing would result in tensile cracking of the dam body which is one of the most important failure modes of concrete dams [³⁰]; however it is dependent on the tensile strength of the concrete. So, the limit-state function Y, is defined as

(28)

where f_t is the tensile strength of concrete, and τ_max is the maximum principal thermal stress of the dam body during the transient thermal analysis. The state of Y<0 indicates that the limit-state function is exceeded or the failure scenario (tensile cracking) occurs, however, this failure scenario is not coincident with the total failure of the arch dam.

Probabilistic analysis results

In this section, the results of the probabilistic analysis methods are separately presented. Also the sensitivity of the random variables is investigated.

FOSM

First, the exceedance probability of limit-state function, P_f [f_t– τ_max(X)<0], is computed using the FOSM for varying values of concrete tensile strength f_t. The obtained P_f - f_t plot considering all random variables (30 R.V.) is shown in Fig. 2. From this figure, the probability of development of thermal tensile stress within the dam body is about 70%. But, the probability of exceeding the thermal tensile stress of 1.5 MPa within the dam body is below 1%. Moreover, it is observed from this figure that the exceedance probabilities obtained for the model with 30, 7, and 5 random variables are the same. Consequently, it can be concluded that it is reasonable to consider the model with only 5 random variables for further analysis to reduce the computational cost of the analyses and also the complexity of the model. In addition, this figure reveals the exceedance probabilities of the concrete of the dam body with allowable tensile varying from 1 to 3 MPa which can be used by not only the researchers but also practical engineers.

The contribution of all random variables to the variance of the limit-state function is computed using the FOSM importance vector; some of its measures are presented in Fig. 3 with blue color with their related values. In this figure, T_w,h is the mean annual water temperature at depth h which is dependent random variable. So T_w,0 is equal to T_surf in Eq. (24). The variables with highest contributions or the most important variables of the structural system properties, from Fig. 3(a), are the concrete Young’s modulus E, and coefficient of thermal expansion α, with the positive sign which shows that an increase in these random variables yields a larger probability of failure. The importance measures for the other structural properties are near or almost zero. About the thermal loading parameters in Fig. 3(b), the most important variables are the closure temperature T₀, and the mean annual air temperature {\overline{T}}_{\mathrm{a}\mathrm{i}\mathrm{r}}, with opposite signs which means that T₀ is a resistance variable while {\overline{T}}_{\mathrm{a}\mathrm{i}\mathrm{r}} is a load variable. From Fig. 3(b), it is observed that the mean annual water temperature in top 20m depth of the reservoir is the most important, below that the annual mean temperatures of the reservoir can be assumed deterministic at their mean values. The amplitude of air temperature is of low importance, and the amplitudes of water temperatures have almost zero importance measure. Also, the importance of loading variables is much higher than the structural properties.

Based on the importance measures, it is decided to reduce the number of random variables to the most important ones, making the values of the rest constant and equal to their mean values. Hence, two new sets are defined: (1) 7 random variables, 7 R.V., including α, E, {\overline{T}}_{\mathrm{a}\mathrm{i}\mathrm{r}}, A_air, T₀, T_w,0, T_w,20; and (2) 5 random variables, 5 R.V., including {\overline{T}}_{\mathrm{a}\mathrm{i}\mathrm{r}}, A_air, T₀, T_w,0, T_w,20. In the latter set, the parameters of structural properties are all considered deterministic. By repeating the FOSM analysis for these two sets, their P_f -f_t plots and importance measures are shown in Figs. 2 and 3, respectively. As it is observed, the results of the two new sets are very similar to those calculated with 30 random variables, so it is concluded that the number of random variables can be reduced to just 5 loading random variables; the other random variables including all of structural variables can be considered constant.

The contour of exceedance probability of the dam body for tensile strength of f_t = 2 MPa, normalized with respect to the maximum value, considering all 30 random variables is illustrated in Fig. 4. The probabilities are calculated for maximum envelope of tensile stress during the transient thermal analysis of the dam. Higher exceedance probabilities are observed along the dam-foundation interface specifically on the downstream face. The highest values belong to the crest level near abutments. The lowest probabilities occur at the middle regions of the dam body on both faces. From this figure, the regions in which the cracking is most probable can be predicted.

FORM

The estimated values of reliability index and exceedance probability obtained from the FORM for varying concrete tensile strength are illustrated in Fig. 5. Also shown are the results of the FORM considering 7 and 5 important random variables, as previously defined. The probability of exceedance of the maximum thermal tensile stress of the dam body from 1.5 MPa, considering all of random variables, is about 0.2%. As it is obvious, the results of 30 and 7 random variables are almost identical. It means that other unimportant variables can be eliminated from the analysis in order to decrease computational cost of analysis. However, they show more probability of failure with respect to 5 random variables set. Therefore, neglecting randomness due to important structural properties α and E, underestimates the probability of failure especially for higher values of f_t, it is because of dispersion for lower exceedance probabilities. Since the FORM sensitivity results are more reliable than what was provided by the FOSM, the important system properties α and E should be considered as random variables in the thermal reliability analysis of concrete arch dams.

The components of the FORM’s sensitivity α-vector considering 30, 7, and 5 random variables are depicted in Fig. 6. Also reported are the values of the 30 R.V. set. The results are computed for f_t = 2.0 MPa; they are approximately the same for other f_t values. The important random variables are the same as those predicted by the FOSM with the same sign, but with a significant difference that the importance of structural properties α and E are in the same order with the loading important variables. So the results of the 5 random variables set differ from 30 and 7 random variables sets which are approximately the same. Therefore, it is found out that the structural properties α and E play important rule in reliability analysis and it is mandatory to be assumed as random variables. The other structural properties are of very low importance. The amplitude of water temperature is again not important at all. The mean annual water temperature is important just for top 20m of the reservoir. From the sign of sensitivity measures in Fig. 6, it can be inferred that as all important random variables except T₀ increase, reliability index decreases. On the other hand, by increasing T₀, the structure response which is based on maximum tensile stress decreases, and the system becomes more reliable (higher value of reliability index).

MC-LHS

The random variables are separately sampled 1000, 2000, and 4000 times, each time generating a set of realizations that are consistent with their assumed probability distribution functions. For each set of realizations, the maximum (tensile) principal stress during the transient thermal analysis is computed. The histograms of the obtained thermal stresses, considering all 30 random variables, are shown in Fig. 7 along with the best distribution fit and its second-moment characteristics determined by the Anderson-Darling test [³³]. Regardless of the number of samples, the best-fitted probability distribution is lognormal with approximate mean and standard deviation of 0.61 and 0.47 MPa, respectively. Although increasing the number of samples would lead to more accurate results, but large coefficient of variation (standard deviation divided by mean) of 0.77 shows high dispersion in the results due to the assumed random variables. It means that large uncertainties are controlling the results, no matter what number of simulations is being used. So it is expected that the assumed distribution functions and the related parameters, specifically those for important random variables, have high impact on the probabilistic and reliability analysis as it has been shown in other applications [¹⁰]. Furthermore, by increase number of realizations, the probability characteristics of the outputs do not change significantly indicating that lower number of samples will provide the same result accuracy as greater number of samples. Hence, for complex system which has high computational cost, similar to what is dealt in this investigations, it is recommended to use efficient and also practical number of realizations. It is worth mentioning that performing more realizations may be impractical as the computational cost of deterministic model is high. In addition, there are several methods in the literature which have dealt with complex engineering models using meta-models, chaos polynomial expansion method, and surrogate models [¹¹–¹⁵]. They can be the subject of future studies.

The failure exceedance probabilities P_f can be consequently computed from the histograms and plotted for different f_t values as shown in Fig. 8. As it is clear, the P_f - f_t curves are very close for different number of model evaluations, N. So, approximately the same estimation for P_f would be obtained by using lower number of samples that demands lower computational cost. The exceedance probability drops significantly for f_t<1 MPa and its decline rate decreases for higher f_t value. The exceedance probability for f_t = 1.5 MPa is about 6%.

Because approximately the same exceedance probability is obtained using lower number of samples, the P_f - f_t curves are re-generated including 7 and 5 random variables, as described before, using 1000 samples. The curves are illustrated in Fig. 9 along with that one considering 30 random variables. From this figure, the curves of 30 and 7 random variables are very close to each other showing more P_f values in every f_t level with respect to the curve of 5 random variables, similar to what was concluded from the FORM. Therefore, the structural parameters of elastic modulus and coefficient of thermal expansion should be considered as random variables while the other structural parameters can be assumed deterministic. The same conclusion can be drawn for the thermal loading parameters.

Comparison of the reliability methods

The P_f - f_t curves generated through various reliability methods considering all 30 random variables are compared in Fig. 10. In each f_t level, the highest exceedance probability of failure belongs to the MC-LHS method while the FORM provides the lowest value. Although the FORM is time-consuming and demands high computational cost, but it is observed the more expensive the reliability method such as the MC-LHS method used, the higher the values of the exceedance probability obtained. Because the MC-LHS method is a more reliable method, which provides more precise results, it is concluded that simpler methods such as the FOSM systematically underestimates the failure probability in entire f_t range. Considering f_t = 2 MPa as a common lower bound for static tensile strength of dams’mass concrete [³²], the probability of exceeding this f_t value is 0.036%, 0.004%, and 1.82% for the FOSM, FORM, and MC-LHS methods, respectively. So the probability of tensile cracking of the dam body under thermal loading in the worst case is below 2%. This is a high failure probability in dam engineering [²³], but it is noteworthy that this thermal cracking does not lead to total collapse of the dam.

In the risk context, the exceedance probabilities conditioned on specific and deterministic thermal loading situation do not provide much information about the dam safety [¹⁶]. To be applicable, these probabilities should be multiplied by the annualized probability of the thermal loadings. In this way, the effects of large scatter in the obtained results can be evaluated on the global results and on the decisions related to the dam safety. The total annualized probability of thermal overstressing failure can be calculated as the sum of the products of the probability of loading by the conditional probability of failure. However, in this study, the probability of the most affecting thermal loadings was directly considered as random variables with proper distributions and parameters. The total probability can be meaningfully compared with the published guidelines [³⁴–³⁶], according to which, the individual risk, considering consequences of the failure computed as the product of the total annualized probability of failure by 1 fatality should be less than 10⁻⁴ yr⁻¹ [¹⁶].

The FOSM importance measures expressed that structural properties are not important which is in contrast with the sensitivity α-vector calculated by the FORM. The results of the MC-LHS method confirmed the FORM results which revealed that only 7 of 30 assumed random variables are important and the other parameters should be omitted as random variable in order to lower computational cost of analyses. Between all random variables, T₀ and {\overline{T}}_{\mathrm{a}\mathrm{i}\mathrm{r}}are the loading parameters that affect results the most in opposite direction. The mean annual temperature of top 20m of the reservoir is important as well, however among these temperatures just {\overline{T}}_{\mathrm{a}\mathrm{i}\mathrm{r}} varies independently. Among structural properties, concrete elastic Young’s modulus and coefficient of thermal expansion are the most important parameters that should be considered as random variable in thermal reliability analysis of arch dams.

Summary and conclusions

The thermal response of concrete arch dams was probabilistically investigated considering uncertainties associated with the structural (material) properties and the thermal loadings. Considering Karaj arch dam as case study, it was numerically modeled along with its foundation in three-dimensional space to investigate its temperature and thermal stress distribution during dam operating phase. The thermal loading included ambient air temperature variation using a sinusoidal function, and reservoir temperature variation employing Bofang’s method. The effect of solar radiation was considered by increasing the ambient temperatures. The tensile overstressing that leads to tensile cracking of the dam body was taken into account as the potential failure mode. The thermal finite element analysis of the dam was combined with three reliability methods: FOSM, FORM, and MC-LHS in order to obtain better thermal response predictions, and estimate the probability of failure in the risk analysis context. The obtained reliability curves showed that in each level of concrete tensile strength, the highest exceedance probability of failure belongs to the MC-LHS method while the FORM provides the lowest value. Since the MC-LHS method is a more precise method, the more expensive the reliability method used, the higher the values of the exceedance probability obtained. It was concluded that simpler methods such as the FOSM systematically underestimates the failure probability.

Investigating the contribution of random variables to the variance of the limit-state function using the FOSM importance vector showed that the most important variables of the structural system properties are the concrete Young’s modulus E, and coefficient of thermal expansion α, such that an increase in these random variables yields a larger probability of failure. About the thermal loading parameters, the most important variables are the closure temperature T₀, and the mean annual air temperature {\overline{T}}_{\mathrm{a}\mathrm{i}\mathrm{r}}, such that T₀ is a resistance variable while {\overline{T}}_{\mathrm{a}\mathrm{i}\mathrm{r}} is a load variable. It was observed that the mean annual water temperature in top 20m depth of the reservoir is important, below that the annual mean temperatures of the reservoir can be assumed deterministic. The FOSM showed that the importance of loading variables is much higher than the structural properties, so the number of random variables can be essentially decreased from 30 to just 5 loading random variables; the other random variables including all of structural variables can be considered deterministic. It means that the other unimportant variables can be eliminated from the analysis in order to decrease computational cost of analysis. But the FORM and MC-LHS methods showed that neglecting randomness due to important structural properties, α and E, underestimates the probability of failure. Since the FORM and MC-LHS are more reliable than the FOSM, the important system properties α and E should be considered as random variables in the thermal reliability analysis of concrete arch dams.

The results of the MC-LHS showed that regardless of the number of samples, the best-fitted probability distribution of the output results is lognormal with large coefficient of variation of 0.77 meaning that large uncertainties are controlling the results, no matter what number of simulations is being used. Also, approximately the same estimation of exceedance probability would be obtained by using lower number of samples, in order of 1000 samples in this type of analysis, that demands lower computational cost.

References

[1] U.S. Army Corps of Engineering. Engineering and Design: Arch dam design, Engineer manual 1110-2-2201, 1994

[2] Wieland M, Kirchen G F. Long-term dam safety monitoring of Punt dal Gall arch dam in Switzerland. Frontiers of Structural and Civil Engineering, 2012, 6(1): 76–83

[3] Sheibany F, Ghaemian M. Effects of environmental action on thermal stress analysis of Karaj concrete arch dam. Journal of Engineering Mechanics, 2006, 132(5): 532–544

[4] Agullo L, Aguado A. Thermal behavior of concrete dams due to environmental actions. Dam Engineering, 1995, VI(1): 3–21

[5] Daoud M, Galanis N, Ballivy G. Calculation of the periodic temperature field in a concrete dam. Canadian Journal of Civil Engineering, 1997, 24(5): 772–784

[6] Léger P, Venturelli J, Bhattacharjee S S. Seasonal temperature and stress distributions in concrete gravity dams (Parts I and II). Canadian Journal of Civil Engineering, 1993, 20(6): 999–1017

[7] Meyer T, Mouvet L. Behavior analysis of the Vieux-Emosson arch gravity dam under thermal loads. Dam Engineering, 1995, VI(4): 275–292

[8] Zhang Z, Garga V K. State of temperature and thermal stress in mass concrete structures subjected to thermal shock. Dam Engineering, 1996, VIII(4): 336–350

[9] Jin F, Chen Z, Wang J, Yang J. Practical procedure for predicting non-uniform temperature on the exposed face of arch dams. Applied Thermal Engineering, 2010, 30(14–15): 2146–2156

[10] Bernier C, Padgett J E, Proulx J, Paultre P. Seismic fragility of concrete gravity dams with spatial variation of angle of friction: case study. Journal of Structural Engineering, 2016, 142(5): 05015002

[11] Vu-Bac N, Lahmer T, Zhuang X, Nguyen-Thoi T, Rabczuk T. A software framework for probabilistic sensitivity analysis for computationally expensive models. Advances in Engineering Software, 2016, 100: 19–31

[12] Vu-Bac N, Silani M, Lahmer T, Zhuang X, Rabczuk T. A unified framework for stochastic predictions of mechanical properties of polymeric nanocomposites. Computational Materials Science, 2015, 96: 520–535

[13] Hamdia K M, Silani M, Zhuang X, He P, Rabczuk T. Stochastic analysis of the fracture toughness of polymeric nanoparticle composites using polynomial chaos expansions. International Journal of Fracture, 2017, 206(2): 215–227

[14] Vu-Bac N, Lahmer T, Keitel H, Zhao J, Zhuang X, Rabczuk T. Stochastic predictions of bulk properties of amorphous polyethylene based on molecular dynamics simulations. Mechanics of Materials, 2014, 68: 70–84

[15] Vu-Bac N, Lahmer T, Zhang Y, Zhuang X, Rabczuk T. Stochastic predictions of interfacial characteristic of polymeric nanocomposites (PNCs). Composites Part B, Engineering, 2014, 59: 80–95

[16] Altarejos-García L, Escuder-Bueno I, Serrano-Lombillo A, de Membrillera-Ortuño M G. Methodology for estimating the probability of failure by sliding in concrete gravity dams in the context of risk analysis. Structural Safety, 2012, 36–37: 1–13

[17] Vu-Bac N, Rafiee R, Zhuang X, Lahmer T, Rabczuk T. Uncertainty quantification for multiscale modeling of polymer nanocomposites with correlated parameters. Composites Part B, Engineering, 2015, 68: 446–464

[18] Liel A B, Haselton C B, Deierlein G G, Baker J W. Incorporating modeling uncertainties in the assessment of seismic collapse risk of buildings. Structural Safety, 2009, 31(2): 197–211

[19] Luísa M, Farinha B, de Lemos J V, Neves E M. Analysis of foundation sliding of an arch dam considering the hydromechanical behavior. Frontiers of Structural and Civil Engineering, 2012, 6(1): 35–43

[20] Raychowdhury P, Jindal S. Shallow foundation response variability due to soil and model parameter uncertainty. Frontiers of Structural and Civil Engineering, 2014, 8(3): 237–251

[21] Nariman N A, Lahmer T, Karampour P. Uncertainty quantification of stability and damage detection parameters of coupled hydrodynamic-ground motion in concrete gravity dams. Frontiers of Structural and Civil Engineering, doi: 10.1007/s11709-018-0462-x

[22] Haukaas T, Der Kiureghian A. Parameter sensitivity and importance measures in nonlinear finite element reliability analysis. Journal of Engineering Mechanics, 2005, 131(10): 1013–1026

[23] Val D, Bljuger F, Yankelevsky D. Reliability evaluation in nonlinear analysis of reinforced concrete structures. Structural Safety, 1997, 19(2): 203–217

[24] Mckay M D, Beckman R J, Conover W J. Comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics, 2000, 42(1): 55–61

[25] Reddy J N, Gartling D K. The Finite-element Method in Heat Transfer and Fluid Dynamics. Boca Raton: CRC, 2001, 5–108

[26] Iman R. Latin hypercube sampling. In: Encyclopedia of Statistical Sciences. Wiley: New York. 1980

[27] Melchers R E. Structural Reliability Analysis and Prediction. 2nd ed. JohnWiley & Sons, 1999

[28] Broding W C, Diederich F W, Parker P S. Structural optimization and design based on a reliability design criterion. Journal of Spacecraft, 1964, 1(1): 56–61

[29] Chauhan S S, Bowles D S. Dam safety risk assessment with uncertainty analysis. ANCOLD Bulletin, 2004, 73–88

[30] Haselton C B. Assessing seismic collapse safety of modern reinforced concrete frame buildings. Dissertation for the Doctoral Degree. San Francisco: Stanford University, 2006

[31] EPRI (Electrical Power Research Institute). Uplift pressures, shear strengths, and tensile strengths for stability analysis of concrete gravity dams. Report No. EPRI TR-100345, 1992

[32] Lo K Y, Lukajic B, Wang S, Ogawa T, Tsui K K. Evaluation of strength parameters of concrete-rock interface for dam safety assessment. In: Canadian Dam Safety Conf, Toronto: Canadian Dam Association, 1990, 71–94

[33] Champagne K, Rivard P, Quirion P M. Shear strength parameters of concrete gravity dams in Quebec. In: CDA 2013 Annual Conf. Toronto: Canadian Dam Association, 2012

[34] Bureau of Reclamation. Guidelines for achieving public protection in dam safety decision making. Technical report. 2003

[35] Australian Committee on Large Dams. Guidelines on risk assessment, 2003

[36] Munger D F, Bowles D S, Boyer D D, Davis D W, Margo D A, Moser D A, Regan P J, Snorteland N. Interim tolerable risk guidelines for US Army Corps of Engineering dams, 2009

RIGHTS & PERMISSIONS

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