Risk Assessment of Aged Concrete Gravity Dam Subjected to Material Deterioration Under Seismic Excitation

For validating the numerical model, the response surface methodology (RSM) according to Myers et al. (1995), analyzes the relationship between several variables (u) and responses (m) of the structure by the following mathematical model:where \(v\) describes the error observed in the response \(m\) and \(f\left( {u_{1} , u_{2} , \ldots u_{k} } \right)\) transmits the response of the structure due to the sets of input variables. In RSM generally, a first-order and second-order polynomial equations are used. Usually, the second order is sufficient to solve the engineering problems and in this study, which is presented as following Eq. (2):here m is the response of prediction and η is the estimated partial regression coefficient; \(u_{i}\) is the coded factor \(\left( {i,j = 1,2,3 \ldots ,k} \right)\) and v is the offset term. The polynomial equation can be used in higher order.

A design experiment tool called central composite design (CCD) (Sadhukhan et al. 2016) is used to predict the output using the equation based on central and axial points with a factorial design for optimization of the response of the structure. Using the Eq. (3), the total experimental number can be created using the CCD tool.here, k is the number of factors and \(c_{q}\) the number of center-point. For accurately amplifying the CCD method in this study, two parameters are used such as the coefficient of Young’s modulus and density.

After identifying the system of the numerical model, it is an essential factor to verify the model with previous studies. The fundamental frequency of the model will be compared with Eq. (4) according to Fenves and Chopra (1986).where \(H\) is the height of the concrete gravity dam. Besides that, the check for the crest spectral acceleration (g) is also a verification factor along with the modal shape. Also, the frequency domain decomposition (FDD) method has been adopted to verify the fundamental frequencies, which is explained more detailed in Sect. 3.4.

From the experimental results, Washa et al. (1989) developed the governing equation for considering the time effects on the compressive strength of a concrete gravity dam, which can be exhibited by Eq. (5).where \(t_{a}\) is the age of concrete in years, and \(f\left( {t_{a} } \right)\) is the compressive strength gained after time. Taking into account the gain in compressive strength of the sound concrete with age (Washa et al. 1989), the value of static elastic modulus in SI is obtained from the following expression (Mandal and Maity 2015):

However, for the external loading and surrounding environmental effects, the concrete material is damaged, and these damages will increase with the time (Kuhl et al. 2004b). This damage is manifested as the porosity of concrete and Eq. (7) shows the total porosity of concrete.here, \(\phi\) is the total porosity, \(\phi_{0}\) is the initial porosity, \(\phi_{\text{c}}\) is chemical porosity, and \(\phi_{\text{m}}\) is the apparent mechanical porosity. The mechanical porosity \(\phi_{\text{m}}\) is defined by the Eq. (8).where \(d_{e}\) is the scalar degradation parameter and the function of this parameter has been expressed (Gogoi and Maity 2007; Mandal and Maity 2015; Nayak and Maity 2013) as the following equation:here, \(k_{m}^{0}\) and \(k_{m}\) are the values of strain that represent the initial damage and is the maximum value of strain during loading history, respectively. If there is no degradation due to mechanical loading, the \(k_{m}\) may be considered as \(k_{m}^{0}\); as a result of the \(d_{e}\) and \(\phi_{m}\) is zero and \(\alpha_{m}\), \(\beta_{m}\) are parameters that have been taken here from Bangert et al. (2003). The value of \(\alpha_{s}\) will differs from 1 to 0 because of the degradation and non-degradation, respectively.

The relation between non-degraded Young’s modulus of elasticity \(E_{0}\) and degraded Young’s modulus of elasticity by the porosity effect of the concrete is \(E_{e} = \left( {1 - {\text{d}}_{e} } \right)E_{0}\) (Mandal and Maity 2015; Nayak and Maity 2013). Therefore, from Gogoi and Maity (2007), the time-varying damaged modulus of elasticity of concrete can be written by the following equation:

For assessing the seismic performance of a CGD with time-varying ageing effects, Bohyeonsan multipurpose CGD has been selected. This dam is located in the upper stream of Gohyeoncheon, which is the second tributary of the Kumho River in South Korea. Figure 2a presents the location of the sensors to get earthquake measurement data and Fig. 2b shows the sectional detailing. The dam belongs to the total crest length is 250 m and the maximum height is 57 m. This dam significantly is used for the controlling of reservoir water, the full storage capacity of the reservoir is 22.10 × 10⁶ m³ and the construction of the dam was completed in 2014. The crest width of the dam is 11.15 m and the height varies from 34.5 to 57.0 m. Table 1 shows the specification of this dam.

For seismic analysis of the Bohyeonsan dam, a two-dimensional finite element model is presented here by using ABAQUS. The FEM for the selected section (Fig. 2b) from the 3D view of the dam (Fig. 2a) is illustrated in Fig. 3. The sectional view is shown in Fig. 2b; it can be seen that this dam is built using two kinds of concrete with different elastic modulus. The dam material property was taken from Table 1 and the dimensions are shown in Fig. 3 as well as the mesh distribution. The mesh size in the model was assigned in such division that the number of finite elements for the concrete inside and outside material was 500 and 358, respectively. The FEM consists of 4 nodes, bi-linear, plane strain quadrilateral elements (CPE4R) (Fig. 3) considering reduced integration and hourglass control (Al-Shadeedi and Hamdi 2017).

The non-linear dynamic analysis was carried out by adopting the implicit integration method accounting with the gravity load due to its self-weight as a static condition and ground horizontal acceleration of selected earthquakes as the seismic condition. The upstream wall was subjected to the water pressure up to 42.82 m, where the interaction between dam and water is considered here as a dynamic effect resulting from the transverse component of ground motion. This was simply modeled as added mass at the interface of dam–water and calculated using the Westergaard (1933) formula, which is also used in several studies (Alembagheri and Ghaemian 2013; Ansari et al. 2018; Nguyen et al. 2019). According to the Westergaard (1933), in Fig. 3 to assign the water pressure, the added masses were taken at each node (25 nodes) at the interface of the dam and reservoir using the following equation:

The vertical hydrodynamic components due to the ground motion were neglected in the simulations and the rigid foundation was used for bedrock condition. To consider the free-field motion during an earthquake, it was applied at the dam base as shown in Fig. 3. Here, only the horizontal ground motion data are considered for seismic analysis of the dam (Alembagheri and Ghaemian 2013). For the optimization of the model, which will be explained later, the free-field data of the Pohang earthquake were used as the input ground motion.

The damping matrix, according to the Rayleigh method (Chopra 2011), is adopted here, applied by Mridha and Maity (2014). Considering 5% damping ratio in both inside and outside concrete, the damping coefficients can be calculated by a linear combination of the stiffness matrix [K] and the mass matrix [M] as follows:where \(\alpha\) and \(\beta\) are the mass-proportional and stiffness-proportional coefficients, respectively.

The following dynamic equation of motion can explain the above two-dimensional discretized FEM system.where \([M]\), \([C]\) and \([K]\) are the mass, damping and stiffness matrix, respectively. \(\left\{ u \right\}\) is the displacements vector of the nodal point relative to the free-field ground displacement at the dam base, \(\left\{ {\dot{u}} \right\}\) and \(\left\{ {\ddot{u}} \right\}\) are the relative velocity and acceleration vectors, respectively. \(\left\{ {\ddot{u}_{g} } \right\}\) is the free-field ground acceleration and \(\left\{ p \right\}\) is the total pre-seismic load associated with the gravity and hydrodynamic added mass.

This FEM system was taken all through the seismic analysis as well as structural system identification of this study.

For non-linear analysis of the material model, the concrete damage plasticity model (CDP) was considered. This model can be expressed as the complete inelastic potential behavior, which also can develop proper damage simulation for concrete both in tension and compression. In addition, this model can analyze the concrete structure under the loading combinations both static and dynamic and, thus, enable the transfer of results between the two (Wahalathantri et al. 2011). In the CDP model, the post-failure behavior under compression is defined by a softening stress–strain response and tension stiffening is specified either by means of post-failure stress–strain behavior in tension or by applying a fracture energy cracking criterion.

The CDP model describes that the concrete has significant volume change, when subjected to severe inelastic stress states, commonly referred to as dilation. In this study, the dilation angle has been taken as 36°, while default values were assumed for all other plasticity parameters.

The origin of the non-linearity can be introduced to various system properties such as materials, geometry, non-linear loading, and constraints. To meet the non-linear property, some material parameters are induced as the input data in Table 2.

According to the EN1992-1-1, stress–strain behavior of plain concrete in uniaxial compression is defined as the typical stress–strain relationship for nonlinear structural analysis of concrete. For introducing this the equations are followed by this for concrete compression behavior from EN1992-1-1, where the relationship between the compressive stress, \(\sigma_{c}\) and shortening strain, \(\varepsilon_{c}\) for short-term uniaxial loading is described by the following equation:where \(\sigma_{\text{c}}\) is the compression stress in concrete, \(\eta = \frac{{\varepsilon_{\text{c}} }}{{\varepsilon_{{{\text{c}}1}} }}\), \(\varepsilon_{\text{c}}\) is the compressive strain in the concrete, \(\varepsilon_{{{\text{c}}1}}\) is the compressive strain in the concrete at the peak stress \(f_{\text{c}}\) and \(k = \frac{{1.05E_{\text{cm}} \left| {\left. {\varepsilon_{{{\text{c}}1}} } \right|} \right.}}{{f_{\text{cm}} }}\). Figure 4a shows the uniaxial compression stress–strain behavior of the outside concrete material of the Bohyeonsan dam for each year (0, 10, 20, 30, 40 and 50 years) because the outer material is more vulnerable.

In the case of a tension stiffening approach for concrete exponential tension softening model was used (Cornelissen et al. 1986). This is one of the ways of concrete softening response using a fracture energy concept. Therefore, the post-failure tensile behavior is defined using the following exponential function:where \(w\) is the crack opening displacement, \(w_{c}\) is the crack opening displacement at which stress can no longer be transferred, \(c_{1}\) and \(c_{2}\) are material constants for normal concrete. Figure 4b shows the material softening behavior in tension for each year similarly as the behavior in compression and the maximum \(\sigma_{t}\) follows the splitting strength of concrete using the Eq. (11). As the cracking was started just after this tensile stress, these values were taken for the tensile damage limit states in seismic fragility analysis. From Fig. 4, it is clarified that the damage input parameters showed the effect of fracture behavior along with the effect of degradable material property. Chemo-mechanical model changes the modulus of elasticity with time and produced micro-crack propagation, which causes the tensile crack in the CGD body. As we saw the fracture behavior of the concrete in Fig. 4, it indicates how much crack displacement will dominate the concrete strength as well as the concrete durability. For the seismic capacity evaluation, the concrete tensile damage is taken for showing the failure probability with damage consideration.

As shown in Fig. 2b, the value of E for inside and outside of the dam is different. For decreasing the number of runs and to keep the same ratio of inside and outside parameters, the same coefficient of E (CoE) (i.e., a multiplying factor of E which will be used to get the original value of E) was taken for CCD. Therefore, the numerical model parameters are generated by considering two variables such as CoE and \(\rho\), respectively. The acceleration on the top of the dam under Pohang-earthquake has been counted as a structural response. Pohang earthquake is one of the strongest recorded earthquakes in the Korean Peninsula with magnitude 5.5, which occurred on November 15, 2017 (Grigoli et al. 2018).

CCD has created a total of 9 points using the Eq. (3) and after optimization, the final value of CoE and density is 0.787 and 2.32 (tone/m³), respectively. However, the seismic analysis was then carried out using the optimized parameters enlisted in Table 3.

To understand the validation by RSM method, Fig. 5 shows the response spectrum at the top of the dam before and after optimizing. It is observed that the difference between the peak acceleration and frequency is decreased when compared with the recorded data. By analyzing Fig. 5 and Table 4, it is shown that the response of the dam after RSM is not exactly matched because of many uncertain factors. However, if we consider the percentage of similarities, we can say that the results are acceptable.

After validation of the FEM, modal identification was verified here by comparing the fundamental frequencies with the previous study and existing method. The fundamental frequencies were observed from the optimized FEM simulation and the recorded data were extracted using frequency domain decomposition (FDD) methods. The FDD is a technique for the decomposition of the system response from recorded data to identify the fundamental parameters described in Brincker et al. (2000). This technique follows simple decomposition each of the estimated spectral density matrices. The singular values of the power spectral density (PSD) function matrix \(S_{yy} (\omega )\) are used to estimate the natural frequencies instead of the PSD functions themselves as follows:where ∑ is the diagonal matrix consisting of the singular values (\(\sigma_{i}^{'} s\)) and U and V are unitary matrices. Since \(S_{yy} (\omega )\) it is symmetric, U becomes equal to V (Ko et al. 2009).

From the FDD extraction, the fundamental frequencies were acquired analyzing the recorded data in Fig. 6. The analysis shows two peaks that were observed through the resonant frequencies and corresponding fundamental frequencies are listed up in Table 5.

Alongside this, the 1st fundamental frequency of the optimized FEM (Fig. 7) is compared with Ref. 1 and Ref. 2 (Table 5). From Table 5, it is observed that for the first mode, the rate of accuracy is 1.4% with Ref. 1 and 5.83% with Ref. 2, respectively. Similarly, for the second mode, the result claims an accuracy of around 4.55% with Ref. 1. Inert to be acceptance of FEM result verification, acceptance value is less than 15% (Shah 2002), where this study shows the most approvable result. Also analyzing the CoV in the last column of Table 5, it can be said that the FEM result has a good agreement with the FDD result and also with the previous study (Eq. 4). The acceptable result for CoV was taken here for verification according to Pakzad (2018). Therefore, the FEM model is validated and verified now for further analysis.

To determine the damaged modulus of elasticity (E_e), a degradation function d_e is calculated from Eq. (9). In this equation, material parameters are taken as \(\alpha_{m} = 0.9\), \(\beta_{m} = 1000\), \(\phi_{0} = 0.2\), and \(k^{0} = 0.00011\) to consider the deterioration effect of concrete material. The value of chemical porosity \(\phi_{c}\) considered in this study is 0.2 (Gogoi and Maity 2007; Kuhl et al. 2004b). The reduction of the modulus of elasticity due to porosity with the varying time has been calculated, using Eq. (10), where the sound modulus of elasticity is calculated using Eq. (6). A graphical representation is shown in Fig. 8 using the sound and damaged modulus of elasticity of the Bohyeonsan Dam. From Fig. 8, it is observed that the elastic modulus of sound concrete is increased with increasing time and after considering the chemo-mechanical damage, the elastic modulus is decreased with time.

Ground motion randomness is carried out by taking 30 unscaled earthquake datasets from k water organization in Korea. Due to the presentence of uncertainty in ground motion, the target spectrum is obtained from KDS 41 and soil class S1 (rock type soil) (KCSC 2019). The selection criteria are followed here according to the design spectra explained in Manandhar2b and Cho (2018). However, the vertical component is ignored here because of comparatively less acceleration than other components. The strongest horizontal component is taken here from the field data record based on the strong duration, i.e., 5–75% area intensities. Therefore, Table 6 outlines the set of selected normalized earthquakes along with their detailing and Fig. 9 will show the clear effect of the comparison of input data with design target spectra.

The normalization of the natural ground motion data set is the way to avoid unwarranted variability. Here, 30 ground motions data have been normalized by multiplying the factor calculated concerning PGA. Scaling of each ground motion is carried out by a scale factor according to Ansari and Agarwal (2016) and Vamvatsikos and Cornell (2002). A set of normalized earthquake data records to be collectively scaled upward or downward and the range of this scale factor depends on the failure of more than 50% damage of the structure (ATC and FEMA 2009). Approximately, 300 numerical analysis has been done for taking the output of all required IM and for each specified year.

Two different ground motion IM are used for plotting the IDA curves. These are the peak ground acceleration (PGA) and spectral acceleration at the structure’s first mode period (Sa).

Construction of fragilities for potential failure can be solved by observing the more severe limit states (Ellingwood and Tekie 2001) as mentioned above Sect. 2.3.1. Among several modes of failure criteria, two vulnerable points were taken to select the limit states for the seismic performance of this study. By analyzing the seismic result of the FEM, it was observed that the LS1 could be the neck or at the foundation zone. Even though the tensile stress is less than the compression, it may cause the crack of the dam body, which will be a significant issue that happens by investigating the relatively elastic and plastic strains on that zone. Therefore, from the parametric study in Tekie and Ellingwood (2003) and the result analysis, LS1 denotes the tensile damage at the foundation (heel of the dam) zone. The non-linear CDP consideration was captured in the cracking propagation in the dam body, where Fig. 4b presents the tensile softening behaviour of the taken CGD for the case study. The tensile damage follows the cracking length and the maximum tensile stress shows the first crack propagated identification. This is denoted as the splinting strength of concrete \(f_{\text{sp}}\).

Note that, for showing the time-dependent seismic performance, each specific year was adopted for calculating LS1, which will change according to the Eq. (11) (Mirza et al. 1979), because the splitting strength of concrete \(f_{\text{sp}}\) is correlated with the material modulus of elasticity (E). The time-dependent change in modulus of elasticity and corresponding splitting strength (LS1) of concrete is listed in Table 7.

The splitting strength is reduced with time as the modulus of elasticity is also reduced by the chemo-mechanical effect of concrete.

In the case of other limit states to get the threshold value of IM, the relative displacement on top of the dam with respect to heel is introduced here as LS2. The LS2 was taken for this dam 1.6 cm (0.028% of the monolith height of the dam), which had been remained constant throughout the seismic analysis of this study.

To assess the chemo-mechanical effect on seismic vulnerability, the fragilities are estimated at different time points for the service life of the CGD. Figure 10 show the seismic performance for 0, 10, 20, 30, 40, and 50 years of the Bohyeonsan dam in terms of peak ground acceleration (PGA) and elastic pseudo-spectral acceleration (Sa), respectively. From Fig. 10a, it can be seen that for an example, the probability of tensile damage in the dam body is about 31% for an earthquake with a PGA of 1 g. Figure 10b shows the spectral acceleration approximately 4 g for the same percentage (31%) of failure probability in case of tension damage. These values are noticed for the zero years as well. The fragilities for LS2 are delineated corresponding to relative deformations of 1.6 cm (calculated as 0.028% of the monolith height) (Sen 2018; Tekie and Ellingwood 2003).

Figure 11 presents the ageing effect by the fragility performance with the HCLPF (Reed and Kennedy 1994) point for each specified year. The result shows a significant amount of change in IM for HCLPF points in the next 50 years. The PGA even Sa looks more critical for 5% failure probability in LS1 than LS2, where the main cause remains on the non-linear material property (NLMP) for analysis. Even though this study shows the seismic fragility analysis using the 30 selected earthquakes in Korea, but it can be updated with different ground motions for other CGDs. In that case, the procedure described in the whole manuscript should be followed in the same way.

To determine the CAV capacity for the aged CGD, all earthquake data sets (taken in this paper) are scaled with the smallest HCLPF PGA (Cao et al. 2019). The estimation of the CAV is to calculate the unscaled ground motion dataset by the threshold PGA. However, this PGA value is observed from Fig. 11a, where it presents a full form of failure probability of up to 50 years. The HCLPF PGA for two limits states is observed as like LS1 < LS2, and these are 0.27 g < 0.3 g, 0.26 g < 0.29 g, 0.23 g < 0.26 g, 0.21 g < 0.23 g, 0.19 g < 0.21 g and 0.17 g < 0.19 g for 0 year, 10 years, 20 years, 30 years, 40 years and 50 years, respectively. In each year for Bohyeonsan CGD, the tensile strength failure state gives the smallest PGA with comparing the relative displacement failure state. Therefore, for calculating the CAV, this smallest PGA gives safety measurement for the structure.

In the previous study, the estimation of CAV directly from the fragility analysis of structure in terms of CAV (Mazılıgüney et al.), which will give the most conservative result to quantify the seismic risk of the dam. Here, for each year, the capacity CAV has been calculated as the mean CAV value from selected earthquakes data sets. Table 8 is presenting the variation of CAV value with the time of 30 earthquake data sets.

Using the value in Table 8, Fig. 12 shows the three-dimensional normal distribution and linear regression analysis to give a model for CAV_limit capacity of Bohyeonsan CGD for each specified year. The model comes from the linear regression (LR) analysis depends on the time (years), where this time has a relation with the degradable young’s modulus of elasticity (E_e) and the model is expressed by the following equation:

This LR model has a minimum error of \(R^{2}\) is 98%, which is acceptable fitting of the normal distribution. Using this capacity model, the result is gradually decreased with time, wherein the given entire time the capacity of CAV_limit value will vary from 0.61 to 0.36 g-s for 0 years to 50 years that is decreased by up to 40% from the present condition. As a result, the engineer or CGD operator can use this equation to predict the condition of structure at any time as well as assess the seismic risk of the degradable concrete gravity dams.