J. Cent. South Univ. (2017) 24: 178-189

DOI: 10.1007/s11771-017-3419-7

Nonlinear finite-element-based structural system failure probability analysis methodology for gravity dams considering correlated failure modes

HU Jiang(胡江), MA Fu-heng(马福恒), WU Su-hua(吴素华)

State Key Laboratory of Hydrology-Water Resources and Hydraulic Engineering, Nanjing Hydraulic Research Institute, Nanjing 210029, China

 Central South University Press and Springer-Verlag Berlin Heidelberg 2017

Abstract: The structural system failure probability (SFP) is a valuable tool for evaluating the global safety level of concrete gravity dams. Traditional methods for estimating the failure probabilities are based on defined mathematical descriptions, namely, limit state functions of failure modes. Several problems are to be solved in the use of traditional methods for gravity dams. One is how to define the limit state function really reflecting the mechanical mechanism of the failure mode; another is how to understand the relationship among failure modes and enable the probability of the whole structure to be determined. Performing SFP analysis for a gravity dam system is a challenging task. This work proposes a novel nonlinear finite-element-based SFP analysis method for gravity dams. Firstly, reasonable nonlinear constitutive modes for dam concrete, concrete/rock interface and rock foundation are respectively introduced according to corresponding mechanical mechanisms. Meanwhile the response surface (RS) method is used to model limit state functions of main failure modes through the Monte Carlo (MC) simulation results of the dam-interface-foundation interaction finite element (FE) analysis. Secondly, a numerical SFP method is studied to compute the probabilities of several failure modes efficiently by simple matrix integration operations. Then, the nonlinear FE-based SFP analysis methodology for gravity dams considering correlated failure modes with the additional sensitivity analysis is proposed. Finally, a comprehensive computational platform for interfacing the proposed method with the open source FE code Code Aster is developed via a freely available MATLAB software tool (FERUM). This methodology is demonstrated by a case study of an existing gravity dam analysis, in which the dominant failure modes are identified, and the corresponding performance functions are established. Then, the dam failure probability of the structural system is obtained by the proposed method considering the correlation relationship of main failure modes on the basis of the mechanical mechanism analysis with the MC-FE simulations.

Key words: gravity dam; structural system failure probability; nonlinear finite element; response surface method; computational platform

1 Introduction

A large amount of field, laboratory and analytical investigations has been conducted on the structural failure mechanism of gravity dams on rock foundations. It is shown that the structural system failure of these dam-interface-foundation systems are controlled by several factors including the properties of dam concrete, interface and foundation, the relative strengths and states etc [1-3]. Sliding failure is only possible when the contact surfaces with low shear resistance, such as intensely weathered surfaces, large tectonic joints, faults or clay interlayers. If a suitable site with adequate shear parameters cannot be found, then, prior to construction of the dam, a rock mass with such discontinuities are either removed or reliably strengthened by grouting, concreting, anchorage and other methods. Investigations demonstrate that neither sliding nor overturning causes the failure of gravity dams in the absence of weakened foundation surfaces. Instead, failures usually occur due to upstream rupturing and downstream crushing of the foundation, followed by overturning of the structure [1-3].

Structural failure probability analysis based on multi-source information is a valuable tool for safety evaluation and management decision [4, 5]. Traditionally, the deterministic analytical limit state function approach is used for such analyses. Generally, these analytical limit state functions were based on empirical formulae or lots of assumptions [6, 7]. But this is only sufficient to consider the overly simplified structural models, for example, frame structures due to the easily understandable mechanical mechanisms. Comparatively, however, for 3-D continuum structures, for example gravity dams, the system failure probability problem becomes much more difficult and complicated, and is still in its infancy up to now to the authors’ knowledge. This is owing to two reasons. Firstly, it is more difficult to obtain approximate performance functions (or safety margins) for dominant failure modes and, secondly, correlation may exist among these failure modes. Meanwhile this methodology is becoming increasingly used worldwide due to more and more deteriorated and obsolete engineering structures. After the previous two difficulties have been overcome successfully, it is possible to evaluate overall system failure probabilities (SFP) of 3-D continuum structures, for example, gravity dams. Over the past years, SUDRET et al [8, 9] coupled the finite element (FE) analysis with the probabilistic analysis tools to obtain structural failure probability for 3-D continuum structures. The former is to obtain their performance functions through FE simulations. However common simulation methods, even those incorporating variation reduction techniques, usually involve a very large number of FE analyses to achieve acceptable accuracy. SACHDEVA et al [10, 11] developed a directional approach significantly improves the efficiency of directional simulation by utilizing deterministic point sets to preserve the underlying joint probability distribution of the random vector describing the uncertain factors. Comparatively less effort has been made on gravity dams regarding the concrete dam, the dam/rock interface and the foundation as a structural system.

The objective of this work is to develop an efficient procedure that can be used to estimate the SFP of gravity dams to obtain approximate performance (or limit state) functions considering the dam-interface-foundation interaction with the help of nonlinear FE analysis using reasonable nonlinear constitutive modes and the probabilistic analysis tool.

2 Nonlinear FE-based SFP analysis methodology for gravity dams

2.1 Failure causes and modes for gravity dams

Structure failure is the state or condition of not meeting a desirable or intended objective. As summarized by International Commission on Large Dams (ICOLD), concrete dams fail for one or combinations of the following factors: deliberate acts of sabotage, structural failure of materials used in dam construction, movement and/or failure of the foundation supporting the dam, settlement and cracking of concrete dams, overtopping caused by floods that exceed the discharge capacity, inadequate maintenance and upkeep [12]. Among these factors, the internal erosion of foundation, the deterioration of dam concrete and these induces lower capacity are the most causes dam failures. A concrete dam may withstand significant overloading or overtopping and the limiting factor in such conditions is the erosion of the foundation or dam concrete. Structural failure is usually due to structural deficiencies, weak foundation or sabotage. Correspondingly, the common failure modes of gravity dams are: a) Sliding. Sliding of the whole dam section or a monolith, or part thereof along the concrete-rock interface, lift joints (construction joint) in the dam body or along weak planes in the foundation. b) Overstressing. Ultimate stresses exceed ultimate strength in foundation or dam body. c) Overturning. Overturning of the whole monolith or part thereof. This failure mode is not explicitly accounted for in the current guidelines in China. Overstressing in itself is, strictly speaking, not a global failure mode. Local crushing or rupturing does not necessarily lead to global failure, but it is the cause and initiating factor leading to a global failure. Therefore, sliding as well as overstressing of several elements forming a global failure path is analyzed as common failure modes herein.

2.2 Nonlinear constitutive models for numerical analysis for structural systems of gravity dams

The rigid body analysis, where the dam body, concrete-rock interface, and rock foundation are considered and analyzed as separate parts, does not account for the interaction effect. Thus, it can only roughly reflect the behavior and possible failure modes for gravity dams. The nonlinear effects including material ageing deterioration and foundation material dissolution may impact the performance of structural systems to withstand various challenges from operation, environment and natural events such as hydraulic fracture and earthquake. Therefore, reasonable nonlinear constitutive modes need to be developed to deal with this problem by extending the analysis to an overall structural problem, which is known as system problems.

2.2.1 Developed Mazars damage model for dam concrete

In quasi-brittle materials such as dam concrete, failure is preceded by a gradual development of a nonlinear fracture process zone and by strain localization [13]. Numerically obtained results are not objective with respect to the discretization, and the total energy consumed by the fracture process tends to zero as the computational grid is refined. In past, simple remedies are dependent on the discretization parameters, e.g., on the size and type of finite elements, on the orientation of the crack with respect to the elements, and on the order of integration scheme. This leads to globally realistic responses, but no information about the strain and damage distribution within the fracture process zone is obtained. Mathematically, the localization returns the problem to be solved badly posed because softening causes a loss of ellipticity of the differential equations which describe the process of strains. The numerical solutions do not converge towards physically   acceptable solutions in spite of refinements of mesh.

Advanced regularization methods introduce an additional material parameter—the characteristic length, which is related to the size and spacing of major inhomogeneities and controls the width of the numerically resolved fracture process zone. Such methods can be based on nonlocal averaging, higher-order gradient theories, continua with microstructure, or viscous regularization. Thus, a regularization method is necessary. In this work, the developed Mazars damage model using the nonlocal continuum damage theory is represented herein. The state of the stresses in a material point cannot be any more only described by the characteristics at the point but must also take into account its environment.

The simple damage model with an equal degradation of the bulk and shear moduli postulates the stress–strain law in the form

                               (1)

where σ is the stress; ε is the strain; A is the elastic material stiffness matrix, assumed to be isotropic; D is a scalar damage variable. For undamaged material, D is set to zero, and the response is linear elastic; for D=1, the stiffness vanishes, which corresponds to a completely damaged material.

In the Mazars damage model, a damage function D, controlling the evolution of damage, is introduced similar to the plasticity theory. The damage function has the form

                            (2)

where β is a coefficient which is used to describe shear behavior, and is fixed at 1.06 here to improve behavior in shears; the coefficients α_t and α_c define a restrain between the damage caused by tension and compression.

To realize the nonlocal continuum damage model, the deformation gradient is chosen here, and a regularized strain tensor  is thus used to check the characteristic equation

                               (3)

where the scalar L_c is the characteristic length, and is determined to the size of the domain and the shape of the nonlocal weight function, and related to the material microstructure. In the nonlocal version of the model, the ‘‘local’’ equivalent strain, ε_eq, is replaced by its weighted average.

               (4)

where  is a given nonlocal weight function; l_c is an internal length and can be estimated at three times size of the largest aggregate of dam concrete), from the practical point of view, lc4Lc; is the representative volume at the point x; Ω is the volume of structure.

2.2.2 Cohesive zone model for concrete-rock interface fracture

The mechanical behavior of dam-foundation interfaces plays a key role in concrete gravity dam engineering since it is probably the weakest part of the structure and therefore the evolutionary crack process occurring along this interface determines the global load- bearing capacity. The energy dissipated during crack propagation should be considered during crack propagation at the interface between concrete dam and rock foundation. Meanwhile, water pressure plays a significant role in the crack propagation since the fluid pressure distribution can promote further propagation, and this hydraulic loading depends on the degradation of the fracturing surfaces that increases the crack opening. To determine the stress state and the crack propagation at bimaterial interfaces between dam base and foundation rock more accurately, the problem is discussed below in the framework of the cohesive crack models considering the uplift pressure effects for quasi-brittle materials [14, 15].

As shown in Fig. 1, when the cohesive surfaces separate, traction first increases until a maximum is reached, and then subsequently reduces to zero, which results in complete separation. The variation in traction in relation to displacement is plotted on a curve, which is called the traction-displacement curve. The area under this curve equals to the energy needed for separation. The traction-displacement curve gives the constitutive behavior of the fracture. The amount of fracture energy Ф dissipated in the region depends on the shape of the model (proportional to its length in 2D and its surface in 3D) considered. Thus, the surface fracture energy can beexpressed by this function:

                              (5)

where Г is the geometrical shape (length in 2D and surface area in 3D) of the crack model.

Fig. 1 Cohesive zone model for concrete-rock interface fracture

The displacement field u in a mechanical equilibrium state can be solved by minimizing the total energy of the elastic strain energy Ψ, the surface fracture energy Ф, and the work of the external forces W_ext. The solution can be obtained by using a variational approach:

                   (6)

The interface surface is discretized in 2D or 3D by thin finite elements in this paper. The opening displacement of the element is a linear function of nodal displacements, i.e.,(δ_n and δ_t are respectively the normal and the tangential components), and the cohesive force exerted on the crack lips, noted by σ, is defined by derivative of the surface energy density compared to the opening displacement.

A simple model, based on linear softening for the two pure-mode decohesion laws, as schematically depicted in Fig. 2, is adopted. Denoting k₀ and σ_max as the first cracking relative displacement component and the peak value of the traction component associated to the fracture respectively. The interface behaves initially in a linear elastic way, then as soon as the normal stress σ_n reaches the breaking value σ_max, it has a nonlinear behavior characterized by a slope of softening of the fracture –K_n/P_rupt, where K_n is the normal stiffness. The first cracking relative displacement component is determined by elastic stiffness namely k₀=σ_max/K_n; the critical relative displacement k_rupt is evaluated as by P_rupt namely k_rupt=σ_max(1+ P_rupt)/K_n. The dissipation energy G_f  increases with P_rupt.

Fig. 2 Normal cohesive stress vs normal displacement for interface

In the proposed model, a distribution of an uplift pressure is assumed on the crack. When a crack propagates at the concrete-rock interfaces, the water pressure effect is taken into account by considering a further static pressure acting on the real crack and in the process zone. The value of this pressure is set as a function of the crack opening displacement and it varies from a prescribed initial pressure to the external water pressure.

Assume that the hydrostatic energy is explicitly dependent on the water pressure gradient  namely where C is a parameter. According to Fick’s first law, the hydraulic flow is proportional to the gradient of pressure, and can be given through the pressure gradient field:

                            (7)

The pressure field of the equilibrium state can be obtained by the energy minimization criteria

                         (8)

Based on the Eq. (8), the quasi-static mechanical balance equation is modified, namely  Finally the pressure gradient is coupled with the strain field, namely

At the points where the water penetrates, the normal strain σ_n is changed to σ_n-p. The water pressure distribution can promote further propagation. Meanwhile, the hydraulic flow is a function of the crack opening displacement and it varies from the external water pressure gradient

                               (9)

According to Poiseuille’s law, the flow in the crack opening can be written in the following way:

                (10)

where  and ρ are the dynamic viscosity and the mass density, respectively.

In the case of flow through the concrete-rock interface, significant flows exist even for the closed interfaces. Then it is necessary to define the interface thickness ε_min for significant flows, thus Eq. (10) becomes:

         (11)

Once a crack propagates at the interfaces and at the dam base, a further static pressure acts on the real crack. Meanwhile, due to the opening crack, the flow is easy. Finally, the following coupled relationships are obtained.

 (12)

where the generalized stress is introduced, thus the resolution of the mechanical- hydraulic coupled balance equations is equivalent to the resolution of the quasi-static mechanical problem.

It can be noted that the static water pressure is considered in the model presented by Eq. (12); thus, the proposed approach can be adopted to analyze the behavior of dam subjected to static hydraulic pressure.

2.2.3 Elastoplastic Mohr-Coulomb model for concrete- rock interface sliding

Many studies concluded that some tensile strength exist on the concrete-rock interface before shear sliding failure [16, 17]. The Mohr-Coulomb criterion is used as the yield function and the plastic potential function for concrete-rock interface sliding. The maximal tensile strength R_t for interfaces are assumed to c/μ, where c and μ are the cohesion and the kinetic friction coefficient respectively. With this assumption, the tensile strength is simply related to the cohesion. In addition, an isotropic hardening parameter (noted as K), which makes it possible to regularize the tangential slope in the phase of sliding, is introduced. The introduced elastoplastic model is only related to the tangent part of the constitutive law. There is no plastic part for the normal displacement, the aforementioned is always elastic. The tangential displacement is broken up into an elastic part and a plastic part  the cumulated tangential displacement is noted as λ, and then the elastoplastic model is given as

              (13)

In the elastic region the relation between the opening and the stress is linear, and the plastic tangential parameter does not evolve, namely  However, once it reaches the edge of the sliding cone defined by  the evolution of the tangential plastic opening is governed by the non-aligned flow model represented by Eq. (13). The interval of the tensile strength is (0, R_t), and it is the function of the shear stress, it is null for a shear stress higher than the adhesion parameter c and is worth the maximum if the shear stress is null. It is also supposed that once reached the value of the maximal tensile strength, the normal stress does not evolve any more.

2.2.4 Drucker–Prager model for foundation rock

The Drucker-Prager yield criterion is a pressure- dependent model for determining whether a material has failed or undergone plastic yielding, and has been applied to rock and other pressure-dependent materials [17]. The Drucker–Prager yield criterion has the form:

                     (14)

where  is the equivalent stress (or von Mises stress); A is the given coefficient of the material and can be obtained by the friction angle φ, namely A=2sinφ/(3-sinφ); k(γ) defines the Drucker–Prager yield surface, and is the hardening function of the cumulated plastic strain γ of linear or parabolic type.

For isotropic linear hardening materials, the function k(γ) has the following form:

                 (15)

where σ_Y is the given coefficient representing the material yield strength related to the cohesion c and the friction angle φ, and can be given by σ_Y=(6cosφ)/(3-sinφ), and H is the hardening modulus coefficient.

2.3 Response surface method for FE-based failure probability analysis for complex structures

Response surface (RS) method is a collection of mathematical and statistical techniques that are useful for the modeling and analysis of problems in which a response of interest is influenced by several variables and the aim is to optimize this response [18, 19]. The RS method is used to aid in the solution of certain types of problems which are pertinent to scientific or engineering processes. For complex structures like gravity dams, it is difficult to use analytical limit function to reflect the structural nonlinear response of the system under various loads. The random FE method can be introduced to obtain the typical responses under different loads. In this way, the actual limit state function g(x) is replaced by a polynomial type of function as an example, a quadratic polynomial with cross terms:

   (16)

where a, b_i, c_i and d_ij are the coefficients of the polynomial; n is the number of the variables x; ε is the random error that contains the error due to neglecting the higher order terms.

In practice, it is suggested that polynomials are up to second order for the response function, often without the mixed terms because they are uneconomic for problems involving large numbers of random variables. In some cases, the simplest form of the linear type polynomial is used instead of the actual function. To obtain unknown coefficients, the experimental points should be chosen. There are many sampling methods such as star, full factorial, central composite and Box– Behnken designs to obtain the experimental points. In this study, the full factorial method is used. The experimental points are selected around the mean values according to their standard deviation, namely x_i=μ_iσ_i, where μ_i and σ_i are the mean and the standard deviation of the variables, respectively. The coefficients of the RS function including a, b_i, c_i and d_ij are achieved by using the least square method corresponding to the experimental points.

Then, the negative normalized gradient vector at the so-called Most Probable Point (MPP) of failure which is the limit-state surface closest point to the origin is computed as

                       (17)

where  and u^* are respectively the gradient vector of the limit state function g(·) and the MPP of independent standard normal variables u.

In SFP analysis, these normalized gradient vectors are used to quantify the statistical dependence among the correlated failure modes.

3 Numerical integration method based SFP calculation for correlated failure modes

SU et al [5] developed an analytical method for calculating the SFP of gravity dams with correlated failure modes, which is applicable to series, parallel and general systems. However this method is not flexible in incorporating amount of available information on multi failure modes. For these reasons, the numerical integration method, which is capable of solving general system events with various merits, is developed to solve the SFP easily.

First, consider a system event is the vector of the basic component event e_j (j=1, …, m). Note that the component event e_j should be the mutually exclusive and collectively exhaustive event (MECE), thus it has two distinct states, e.g. the failure and survival. Then, any system event can be represented by an event vector c whose the j-th element is 1 if e_j belongs to the system event and 0 otherwise. Due to the mutual exclusiveness of e_j, the system event probability P_sys is the sum of the probabilities p_j of e_j that belongs to the system event, namely

                 (18)

where c and p are respectively the event vector whose element has 1 or 0 depending on whether e_j belongs to E_sys or not and the probability vector. Both c and p are given as column vectors.

Equation (18) can be generalized to calculate the probabilities of multiple system events under multiple conditions of component failures by a single matrix multiplication, i.e.,

                                 (19)

where P is the matrix whose element at the i-th row and the j-th column is the probability of the i-th system event under the j-th condition, C=[c₁, c₂, …, c_Nsys] is the matrix containing the event vectors of the N_sys system events considered, and P=[p₁, p₂, …, p_Ncond] is the matrix that has the probability vectors of the N_cond different component conditions.

Let X denote the vector containing common source random variables (CSRV), the SFP can be computed by the total probability theorem as

                   (20)

where is the conditional probability of the system event, X=x, and f_x(x) is the joint probability density function (PDF) of X.

Using the proposed formulation in Eq. (19), the SFP in Eq. (20) can be computed as

                  (21)

where p(x) is the vector of the conditional probabilities of the basic MECE events given X=x.

For CSRV, they can be represented by the Dunnett-Sobel (DS) correlation coefficient matrix. The system event definition, represented by the event vector c, is not affected by the outcome of X. Therefore, the SFP is alternatively computed as

                 (22)

Suppose that Z_i (i=1, …, n) are DS class standard normal random variables. This means that the correlation coefficient between Z_i and Z_j is specified as ρ_ij=r_i·r_j (i≠j) and ρ_ii=1. Then, Z_i can be represented by (n+1) independent random variables:

                (23)

where X and U_i (i=1, …, n) are independent standard normal random variables. For a given outcome X=x, therefore, Z_i is conditionally independent of each other. If the safety margins or the limit state functions of component events are represented as or transformed to DS class random variables, the common source effect can be identified by the single random variable X. One special example is that the components are equally correlated, i.e., ρ_ij=ρ, for which,

For a general correlation coefficient matrix, one can try to fit it with a DS class by finding a set of r_i that minimizes the errors between the actual correlations ρ_ij and r_i·r_j. When the differences between the correlation coefficients cause significant errors in SFP estimations, one can generalize the DS class by adding more CSRVs. When m CSRVs, X_k, k=1, …, m is used for accuracy, Z_i is represented as

                   (24)

In this case, the correlation coefficient between Z_i and Z_j is derived as

                   (25)

Applying the condition X=x into Eq. (24), the probabilities of conditionally independent components  or P_i(x) in Eq. (22) can be computed.

For dam safety management, it is helpful to estimate the sensitivities of failure probability P_f with respect to major parameters. For example, the sensitivities with respect to the standard deviations are critical in efforts for uncertainty safety management or danger-removal reinforcement. The following sensitivity-based importance measures are often used to quantify the importance of the uncertain design variables, respectively

                                 (26)

                                 (27)

where μ_i and σ_i respectively denote the mean and standard deviation of the i-th design variable.

4 SFP calculation and analysis platform for gravity dams

The open-source Matlab toolbox FERUM (finite element reliability using Matlab) [20] and the open-source FE analysis software package CA (Code_Aster) [21] are respectively selected as the general purpose structural failure probability calculation and numerical simulation solver tool in structural mechanics and failure probability analysis. The FERUM toolbox was developed by HAUKAAS and KIUREGHIAN from the University of California, Berkeley primarily for structural reliability and failure probability calculation with stochastic FE methods, and can be easily developed for research purposes. The CA was originally developed as an in-house application by the French company EDF, is capable of multi-physics coupling solver in various fields beyond the standard functionalities of a FEM software for solid mechanics, and is especially attractive since it has strong interaction capacity during input and output processes. By using these two software packages specialized in their own areas, it becomes possible to take full advantages of them and solve challenging problems. The efficient transfer of data and information between the two computational tools is facilitated via an intervening interface. Here, the intervening interface is the scripting language MATLAB. Then integrating these two open source codes and the newly developed interface, the FERUM-CA platform is finally developed. The development of an FERUM-CA platform enables to perform failure probability analysis for complex structural system failure in conjunction with sophisticated nonlinear FE simulations. FERUM-CA is employed so that FERUM can obtain the output responses that appear in the limit-state functions (e.g., strain, displacement, stress, force, and energy results) and their gradients from an FE-based computational simulation. Figure 3 outlines an overview of the steps involved in the coupled software platform. It shows various levels of interactions between the two computational tools.

A system command in MATLAB is issued to execute the CA so that the user can easily conduct the structural response analysis using the nonlinear constitutive models for the dam concrete, the foundation rock and the concrete-rock interface described previously. Firstly, the CA input file is updated through a series of dynamic data exchanges using the data obtained from the FERUM. It should be noted that the CA input file is generated using template file containing the dam system model data, parameters of nonlinear constitutive models shown in Section 2.2 and solution information for the deterministic analysis. Thus if any of given parameters is changed, the new numerical model will be regenerated accordingly. Secondly, the input file is fed into CA to compute the structural behavior for the dam-interface- foundation system which is eventually used by the RS method given by Eq. (16) to determine the limit state functions for all potential failure modes summarized in Section 2.1. Then, the numerical integration code for correlated failure modes, also developed using MATLAB based the derivation in Section 3, is used for calculating the dam SFP.

5 Case study of SFP analyses of gravity dam

5.1 Project brief and considered random variables

The proposed platform for FERUM-CA and the interface code are demonstrated by a numerical example of an existing hydropower project in southeastern China. The water retaining structure of this project is a concrete gravity dam and belongs to Grade 1 projects. The elevation at the dam crest is 179.00 m while 80.00 m at the base. The normal storage water level, the design flood level, and the check flood level are 173.00, 174.76 (0.2%), and 177.80 m (0.02%), respectively. This is an annual regulation reservoir. Its tailwater depth is not considered herein. One typical cross section of water retaining dam is selected for investigation, shown in  Fig. 4. Its upstream face is vertical. Its downstream face is vertical above the elevation of 168.00 m and below with the ratio of 1:0.75. Its crest width and base width are 7.00 and 73.00 m, respectively.

Fig. 3 FERUM-CA platform

Considering that this project is in operation, its geometry size and configuration are both regarded as fixed parameters. The mean value and standard deviation of the headwater level H₁ can be obtained from the observed data. The result of the K-S test shows that H₁ is a normal distribution variable, 91.5 m and 1.21 m. Many documents and published literature show that there is a significant correlation among the elastic modulus, the tensile strength and the compression strength of the dam concrete, and the correlation coefficient is assumed to be 0.425. Based on monitoring data, the correlation coefficient between headwater level and the reduction factor of uplift pressure is 0.900. The material properties of the dam, the foundation and the interface used in non-linear analysis are given in Table 1. The hardening modulus coefficient H is fixed at -200 MPa here.

Fig. 4 Cross section of studied dam

Table 1 Variables and COV used in reliability analysis

5.2 SFP calculation based on FERUM-CA platform

Numerical solutions are performed for case analysis which contain geometrically and materially non-linear analyses accounting for kinds of uncertainties by use of various material properties of the dam, the foundation and the interface such as the elastic modulus, the Poisson ratio and the volume-weight of dam concrete, the elastic modulus of foundation rock, the headwater level, the reduction factor of the uplift pressure under the interface, the cohesion and the friction angle (reflected by the yield strength in Eq. (14)), the tensile strength, the friction angle and the cohesion of the concrete-rock interface. The full factorial method is used with the direct coupled method and the RS method in the failure probability analysis. The FE model of the typical dam section is shown in Fig. 5.

CA is performed for determining individual failure mode and single limit-state functions with the RS method. FERUM-CA is employed so that FERUM can obtain the output responses that appear in the limit-state functions and their gradients from FE-based computational simulations. This calculation focuses on main failure modes summarized in Section 2.1. The results of the preliminary deterministic analysis help identify important system failure modes to be considered during system reliability analysis. For individual failure mode, the limit state function is estimated with more than 100 times structural analyses. After structural analyses, dominant system failure modes, consisting of the foundation rock strength failure and the tensile crack along the concrete-rock interface, are selected, and noted as NM1, NM2, NM3 and NM4. The foundation rock strength failure include two types namely the rock under dam toe and under dam upstream heel failures, the failure paths are respectively shown in Fig. 6. For example, the failure path of NM3 could be the failure of element 1→2→3→4→5. For completeness, the sliding failure along the concrete-rock interface is also considered herein, and the failure initiates from the upstream side, and gradually propagates to downstream. Using the RS method, the limit state functions of the tensile crack of the interface, the rock strength failure A (under the dam upstream heel), the rock strength failure B (under the dam downstream toe) and the sliding failure can then be respectively represented by the following functions:

where x₁, x₂, …, x₁₀ are respectively represent the headwater level, the elastic modulus of the foundation rock, the tensile strength of the interface, the friction angle of the foundation rock, the cohesion of the foundation rock, the reduction factor of the uplift pressure, the friction angle of the interface, the cohesion of the interface and the elastic modulus of the dam concrete. The unit, the mean and the standard deviation values of the variables used in the failure probability analysis of the dam are given in Table 1.

Fig. 5 FE model for of dam system including its foundation and interface

Based on the FE analysis results, the sensitivities of the failure probability of individual with respect to statistical parameters failure mode FERUM analysis can be provided. Meanwhile, the normalized negative gradient vectors of the corresponding limit-state function can also be obtained, the results are listed in Table 2. Using these sensitivities, it is possible to compute the sensitivity-based importance measures in Eq. (26) and Eq. (27). In order to compute the SFP in Eq. (22), the correlation coefficient matrix ρ is constructed by the inner product of the negative normalized gradient vectors:

Finally, the overall SFP of this dam under considered loads and the current state is computed as 8.94×10^-3 using the numerical integration method for correlated failure modes.

5.3 Result analyses and discussion

The findings of numerical gravity dam models, where the water retaining structure is placed on a relatively soft foundation with an interactional interface, is that failure occurred in two ways: either by shear fracture along the interface or by the tensile (rupturing) and compressive crack (crushing) failure under principal compressive stresses at upstream and downstream sides of rock foundation. The failure mechanism is the same as the experiment investigations by FISHMAN [1] and KALUSTYAN [22].

Referring to the Unified Design Standard for Reliability of hydraulic engineering structures in China, the values of the tolerable failure probability of the bearing capacity for Grades 1-3 (safety grade) during design reference periods are 1.03×10^-4, 6.8×10^-4 and 3.4×10^-3, respectively. Based on the results, it is observed that the current SFP is higher than the previously mentioned tolerable failure probability. It must be pointed out that the owner of the dam has already taken a series of reinforcement measurements such as foundation anchorage and anti-seepage treatment of dam upstream surface to improve its safety level.

Correlations affect the structural SFP of the dam-interface-foundation system, involving multiple series failure modes. Without consideration of correlations among failure modes, the SFP may be relatively rough to some extent.

Fig. 6 Dominant failure path of NM3

Table 2 Results of failure probability analysis for dominant failure modes

Table 3 Sensitivities of failure probabilities of various modes with respect to major parameters

From Tables 2-3, it is observed that the dominant failure modes are determined by statistical characteristic value of material properties of the dam concrete, the concrete-rock interface, the rock foundation and the relative values. Among these, the mechanical properties of the foundation rock and the integrality of the concrete-rock interface are the key influence factors. In the condition of lower strength of the foundation rock, the system failure is controlled by the strength failure of the rock (rupturing or crushing) to a great extent. Specially, when the rock strength is poor, the stability of the dam-interface-foundation system is determined by the location under the dam upstream heel. From the results of the investigated engineering, both the rock elements under the upstream dam heel and the dam downstream toe are vulnerable when the rock strength is lower.

6 Conclusions

1) A novel method of FE-based SFP analysis for gravity dams subject to multiple correlated failure modes is developed in this paper. As the first step, a nonlinear FE-based structural SFP analysis method is proposed by integrating nonlinear FE simulations with the RS method. In addition, a numerical method based SFP calculation through matrix integration operations for correlated failure modes is developed to perform SFP analysis while accounting for the same factors. Lastly, a new SFP calculation and analysis platform for gravity dams namely FERUM-CA is developed. The proposed method is applied to an existing gravity dam and successfully demonstrated.

2) As a computational platform of the proposed FE based framework, an interface code, FERUM-CA, is developed. In the platform, the failure probability analysis package FERUM repeatedly calls the code CA to obtain structural responses during a nonlinear numerical simulation analysis, and a SFP analysis using the numerical matrix integration method is performed by use of the results of the failure probability analyses in the individual failure mode level. FERUM-CA is a versatile computing platform, because both FERUM and CA are respectively specialized in failure probability and FE analyses. The proposed platform allows computing the probabilities of general systems with respect to parameters based on the monitoring or measured results.

3) The results of the investigated gravity dam show that the structure, properties, and state of the rock mass (natural state, grouted rocks) have a considerable effect on the system response of a gravity dam. Continuous low strength rock material is readily to fail with the formation of tensile cracks and compressive cracks (crushing and rupturing zones), which develop, respectively, along the principal planes of tensile and compressive stresses. The bearing capacity of gravity dams on poor rock foundations is determined by the material traction/ crushing resistance.

4) Although only a gravity dam is investigated, the proposed FE-based SFP analysis method is also readily applicable to other types of geometrically and mechanically constitutive nonlinear concrete hydraulic structures. The FERUM-CA can also be extended to other structural problems if the response under corresponding loads can be computed by the FE schemes. Such an extension would make FERUM-CA a powerful tool for structural safety evaluation in the dam engineering field.

5) Most of the input parameters of variables of the failure probability model are based on the measured data or engineering archives. Because of insufficient available data, some parameters are still being determined, as referred to in the related research achievements and experiences in this paper. In future work, it will be essential to determine these parameters from field tests and nondestructive inspection.

References

[1] FISHMAN Y A. Stability of concrete retaining structures and their interface with rock foundations [J]. International Journal of Rock Mechanics & Mining Sciences, 2009, 46(6): 957-966.

[2] ZHU Hong-hu, YIN Jian-hua, ZHANG Lin. Physical modelling of sliding failure of concrete gravity dam under overloading condition [J]. Geomechanics and Engineering, 2010, 2(2): 89-106.

[3] GU Chong-shi, WANG Shao-wei, BAO Teng-fei. Influence of construction interfaces on dynamic characteristics of roller compacted concrete dams [J]. Journal of Central South University, 2015, 22(4): 1521-1535.

[4] WANG Chao, ZHANG She-rong, SUN Bo, WANG Gao-hui. Methodology for estimating probability of dynamical system’s failure for concrete gravity dam [J]. Journal of Central South University, 2014, 21(2): 775-789.

[5] SU Huai-zhi, HU Jiang, WEN Zhi-ping. Service life predicting of dam systems with correlated failure modes [J]. Journal of Performance of Constructed Facilities, 2013, 27(3): 252-269.

[6] PEYRAS L, CARVAJAL C, FELIX H, BACCONNET C, ROYET P, BECUE J P, BOISSIER D. Probability-based assessment of dam safety using combined risk analysis and reliability methods– Application to hazards studies [J]. European Journal of Environmental and Civil Engineering, 2012, 16(7): 795-817.

[7] WESTBERG W M, JOHANSSON F. System reliability of concrete dams with respect to foundation stability: Application to a spillway [J]. Journal of Geotechnical and Geoenvironmental Engineering, 2013, 139(2): 308-319.

[8] SUDRET B, KIUREGHIAN A D. Comparison of finite element reliability methods [J]. Probabilistic Engineering Mechanics, 2002, 17(4): 337-348.

[9] HAUKAAS T, KIUREGHIAN A D. Methods and object-oriented software for FE reliability and sensitivity analysis with application to a bridge structure [J]. Journal of Computing in Civil Engineering, 2007, 21(3): 151-163.

[10] SACHDEVA S K, NAIR P B, KEANE A J. On using deterministic FEA software to solve problems in stochastic structural mechanics [J]. Computers and Structures, 2007, 85(5, 6): 277-290.

[11] NIE J S, ELLINGWOOD B R. Finite element-based structural reliability assessment using efficient directional simulation [J]. Journal of Engineering Mechanics, 2005, 131(3): 259-267.

[12] ICOLD. Bulletin 99: Dam failures–statistical analysis [S]. International Commission on Large Dams, Paris, 1995.

[13] LEMAITRE J. A course on damage mechanics [M]. New York: Springer Science & Business Media, 2012.

[14] CARRIER B, GRANET S. Numerical modeling of hydraulic fracture problem in permeable medium using cohesive zone model [J]. Engineering Fracture Mechanics, 2012, 79: 312-328.

[15] ALFANO G, MARFIA S, SACCO E. A cohesive damage–friction interface model accounting for water pressure on crack propagation [J]. Computer Methods in Applied Mechanics and Engineering, 2006, 196: 192-209.

[16] STEFANISHIN D V. Certain theoretical aspects of evaluating aging of hydraulic structures [J]. Hydrotechnical Construction, 1996, 30(9): 534-538.

[17] PAN Peng-zhi, FENG Xia-ting, HUDSON J A. The influence of the intermediate principal stress on rock failure behaviour: A numerical study [J]. Engineering Geology, 2012, 124: 109-118.

[18] KARTAL M E, BASAGA H B, BAYRAKTAR A. Probabilistic nonlinear analysis of CFR dams by MCS using response surface method [J]. Applied Mathematical Modelling, 2011, 35(6): 2752-2770.

[19] LI Dian-qing, CHEN Yi-feng, LU Wen-bo, ZHOU Chuang-bing. Stochastic response surface method for reliability analysis of rock slopes involving correlated non-normal variables [J]. Computers and Geotechnics, 2011, 38(1): 58-68.

[20] FERUM. FERUM User's Guide, Version 3.0 [EB/OL]. [2014-05-01]. http://www.ce.berkeley.edu/projects/ferum/User_s_ Guide/user_s_guide.html.

[21] EDF Group. The aster-full packages (Code_Aster stable and Code_ Aster testing). [EB/OL]. [2014-05-18]. http://www.code-aster.org/ V2/spip.php?rubrique21.

[22] KALUSTYAN E S. Lessons from failures of concrete dams on rock foundations [J]. Hydrotechnical Construction, 1995, 29(2): 13-17.

(Edited by DENG Lü-xiang)

Cite this article as: HU Jiang, MA Fu-heng, WU Su-hua. Nonlinear finite-element-based structural system failure probability analysis methodology for gravity dams considering correlated failure modes [J]. Journal of Central South University, 2017, 24(1): 178-189. DOI: 10.1007/s11771-017-3419-7.

Foundation item: Projects(51409167, 51139001, 51179066) supported by the National Natural Science Foundation of China, Projects(201401022, 201501036) supported by the Ministry of Water Resources Public Welfare Industry Research Special Fund, China; Projects(GG201532, GG201546) supported by the Scientific and Technological Research for Water Conservancy, Henan Province, China

Received date: 2015-06-23; Accepted date: 2015-12-10

Corresponding author: HU Jiang, PhD; Tel: +86-25-858288847; E-mail: huj@nhri.cn