Uncertainty quantification of stability and damage detection parameters of coupled hydrodynamic-ground motion in concrete gravity dams

Abstract

In this paper, models of the global system of the Koyna dam have been created using ABAQUS software considering the dam-reservoir-foundation interaction. Non coupled models and the coupled models were compared regarding the horizontal displacement of the dam crest and the differential settlement of the dam base in clay foundation. Meta models were constructed and uncertainty quantification process was adopted by the support of Sobol’s sensitivity indices considering five uncertain parameters by exploiting Box-Behnken experimental method. The non coupled models results determined overestimated predicted stability and damage detection in the dam. The rational effects of the reservoir height were very sensitive in the variation of the horizontal displacement of the dam crest with a small interaction effect with the beta viscous damping coefficient of the clay foundation. The modulus of elasticity of the clay foundation was the decisive parameter regarding the variation of the differential settlement of the dam base. The XFEM approach has been used for damage detection in relation with both minimum and maximum values of each uncertain parameter. Finally the effects of clay and rock foundations were determined regarding the resistance against the propagation of cracks in the dam, where the rock foundation was the best.

Keywords

massed foundation / hydrodynamic pressure / Box-Behnken method / meta model / Sobol’s sensitivity indices

Introduction

Concrete gravity dams are known to having a complex seismic response because of the interaction of the structure with water from the reservoir and the foundation area. Where most of this type of dams is constructed on rocky foundation, but numerous of concrete gravity dam need to be constructed on soft soil, so in this case, a particular attention to be given to the problems of soil-structure interaction. In 1928, following the failure of St. Francis Dam in California, problems that fail to attract the attention of many large dams led to extensive research in the field. For several events since then for a variety of dams occurred in many parts of the world shows the importance of this issue. Sinfengking Dam in China, Konya dam in India noted that in 1960 it was entered in the quake damage. Another example Sefidrud dam was injured in the earthquake of 1990 Rudbar [¹–³]. When subjected to strong ground motions, mass concrete dams are likely to experience cracking due to the low tensile resistance of concrete. Meanwhile, the potential crack initiation and propagation would adversely affect the static and dynamic performance of dams. As cracks penetrate deep inside a dam, its structural resistance may be considerably weakened, thereby endangering the safety of the dam. In order to meet the ever increasing demand for power, irrigation, drinking water, etc., numerous high concrete dams are being built or to be built and the majority of them are located in active seismic regions. Considering that the possible dam failure due to seismic activities could result in heavy loss of human life and substantial property damages, seismic safety evaluation of high dams remains a crucial problem in dam construction [⁴–⁶] which resulted in development of efficient techniques for damage and crack detection [⁷–¹⁷].

As cracks penetrate deep inside a dam, its structural resistance may be considerably weakened, thereby endangering the safety of the dam. In order to assess the reliability and structural safety of dams, advanced computational methods are commonly applied which are capable of capturing the cracks. These methods can be classified into continuous and discrete fracture approaches. Continuous approaches to fracture smear the crack over a certain width. They include non-local damage models [¹⁸,¹⁹] and gradient models [²⁰,²¹]. Also the introduction of a viscosity [²²–²⁶] smears the crack over a certain width. With the seminal work of Miehe et al. 2010 [²⁷], phase field approaches have become another alternative to non-local and gradient models. Instead of relating the non-local damage to internal state variables, the damage is obtained by solving a partial differential equations. Phase field approaches have meanwhile been applied to numerous interesting problems [²⁸–³⁴] including fracture in thin shells [³³]. Discrete crack approaches include meshfree methods [³⁵–⁴⁰], extended finite element methods [⁴¹–⁴³], smoothed extended finite element methods [⁴⁴–⁴⁶], phantom node methods [⁴⁷–⁵⁰], extended meshfree methods [⁵¹–⁶¹], extended IGA [⁶²–⁶⁵], cracking particle methods [⁶⁶–⁷³], the numerical manifold method [⁷⁴–⁸⁰], efficient remeshing techniques [⁸¹–⁸⁷], multiscale models for fracture [⁸⁸–¹⁰⁰], peridynamics [¹⁰¹–¹⁰³] and dual-horizon peridynamics [¹⁰⁴,¹⁰⁵]. One challenge in these methods is tracking the crack path which can become a daunting task especially in 3D and for dynamic fracture.

Safety evaluation of dynamic response of dams is important for most of researchers. When such system is subjected to an earthquake, hydrodynamic pressures that are in excess of the hydrostatic pressures occur on the upstream face of the dam due to the vibration of the dam and impounded water. Consequently, the prediction of the dynamic response of dam to earthquake loadings is a complicated problem and depends on several factors, such as interaction of the dam with clay or rock foundation and reservoir, and the computer modeling and material properties used in the analysis. In fluid-structure interaction one of the main problems is the identification of the hydrodynamic pressure applied on the dam body during earthquake excitation. Seismic analysis of dam-reservoir-foundation system has a great importance from the vulnerability assessment point of view. The first rigorous analysis of hydrodynamic forces on dam faces during earthquakes was reported by westergaard in 1933 .The modeling and analysis of concrete gravity dam is required for safety evaluation. Reservoir can be modeled by the techniques namely, Wastergaard, Lagrangian and Eulerian approaches [¹⁰⁶–¹¹⁶].

The sloshing oscillations are characterised by the presence of gravity surface waves, which behave in a different manner than the acoustic waves. Gravity waves are nonconservative and their velocity depends on their wavelength. It has been shown that for most concrete gravity dams free-surface waves are negligible. However, in cases where the duration of the excitation is long, the surface wave effect has to be taken into account by imposing an appropriate boundary condition at the free surface of the reservoir [¹¹⁷].

For the structure on the rigid foundation, the input seismic acceleration gives rise to an overturning moment and transverse base shear. As the rock is very stiff, these two stress resultants will not lead to any (additional) deformation or rocking motion at the base. For the structure founded on flexible soil, the motion of the base of the structure will be different from the free-field motion because of the coupling of the structure-soil system. This process, in which the response of the soil influences the motion of the structure and response of the structure influences the motion of the soil, is referred to as soil-structure interaction [¹¹⁸]. The transmitting boundaries or non reflecting boundary conditions (NRBC) are required if one wishes to use exclusively FEM to simulate the semi- infinite soil region. There are two main types of NRBCs: approximate local NRBCs and exact nonlocal NRBCs. The first local boundary indeed the first transmitting boundary was proposed by Lysmer and Kuhlemeyer (1969). It is also known as viscous or absorbing boundary since it places viscous dashpots at the boundary to absorb the energy of the traveling waves [¹¹⁹].

Literature review

Using the basic principles of earthquake engineering, Chopra[¹²⁰] suggested an analytical solution to calculate the hydrodynamic pressure during an earthquake. In his calculations, he assumed a dam as a rigid structure and, for the first time, presented his calculations for actual earthquakes. Besides, unlike previous researchers, he included both the dam’s response to the earthquake and the hydrodynamic pressure together and as an interaction of the dam and reservoirsystem in the frequency domain, and he also included water compressibility effect in his study [¹²¹].

Berrabah et al. [¹²²] concluded that from the dynamic analyses for Brezina concrete arch dam for the three studied cases: dam without soil, dam with mass soil foundation and dam with massless soil foundation show that the presence of the soil in the model develops more stresses in the dam body, especially when the soil mass is considered in the model.

Pekan and Cui [¹²³] performed a comprehensive study on the dynamic behavior of the fractured dam during earthquakes using the distinct element method. Their results showed that the safety of the dam is ensured if the crack shape is horizontal or upstream-sloped; it is very dangerous if the crack slopes downstream.

Some methods have been introduced for the dam–reservoir-foundation interaction. Fenves and Chopra studied the dam-water-foundation rock interaction in a frequency domain linear analysis and developed a computer finite element program called EAGD-84. Then, Leger and Bhattacharjee presented a method which is based on frequency-independent models to approximate the representation of dam-reservoir-foundation interaction. In the work presented by Gaun et al. [¹⁰⁶] an efficient numerical procedure has been described to study the dynamic response of a reservoir-dam-foundation system directly in the time domain.

Khosravi and Heydari [¹¹²] investigated developing 2D Finite Element model different geometrical shape of concrete gravity dam of considering dam-reservoir-foundation interaction effect by using ANSYS modal analysis they proved that Dam soil water interaction is very important for safety design of dams [¹²⁴].

Chopra and Chakrabarti [¹²⁵] studied seismic behavior of dam’s crack path by using linear elastic analysis. The analysis shows, places that are in damage or and risk of the concerning stability of structure. Pal [¹²⁶] was the first researcher who examined Koyna dam by using non-linear analysis. In this research, assuming no effect of reservoir, being rigid foundation, smeared crack model use for crack expansion and strength criteria to crack growth, Koyna dam was analyzed and was shown that the results of material properties and element size are very sensitive.

Most of the previous researches related to the response of the concrete gravity dams were supporting on the theoretical application of the hydrodynamic pressure of the reservoir in the upstream region of the dam system without adopting the real coupling effect of the ground motion with the hydrodynamic pressure due to earthquakes in the analysis and simulations of the finite element applications, where the sloshing behavior of the dynamic water mass is the critical case for the stability and damage detection of the dam especially when the concrete gravity dam is constructed on clay foundation. Thus in this paper the effect of the couplimg is considered and the sloshing of water mass is utilized in the analysis and the simulations by considering the dam-reservoir-foundation interaction. Uncertainty quantification for five decisive uncertain parameters related to the reservoir and the clay foundation would be adopted to determine the rational effects for each parameter on the predicted response of the dam system by utilizing Box-Behnken experimental method to construct the meta models of the dam system response in addition to determining the effect of clay and rock foundations in the control of propagation of cracks in the dam using XFEM approach.

Coupled dam-reservoir-foundation-interaction

Applying the standard Galerkin’s method to simulate of dam-reservoir-foundation interaction is performed. The discretized structural dynamic equation including the dam and foundation rock subjected to ground motion can be formulated using the finite-element approach:

(1){{\displaystyle \mathit{M}}}_{\mathit{s}}{{\displaystyle \ddot{\mathit{u}}}}_{\mathit{e}}+{{\displaystyle \mathit{C}}}_{\mathit{s}}{{\displaystyle \dot{\mathit{u}}}}_{\mathit{e}}+{{\displaystyle \mathit{K}}}_{\mathit{s}}{{\displaystyle \mathit{u}}}_{\mathit{e}}=-{{\displaystyle \mathit{M}}}_{\mathit{s}}{{\displaystyle \ddot{\mathit{u}}}}_{\mathit{g}}+\mathit{Q}{{\displaystyle \mathit{P}}}_{\mathit{e}},

where M_s, C_s and K_s are the structural mass, damping and stiffness matrices, respectively, u_e is the nodal displacement vector with respect to ground and the term QP_e represents the nodal force vector associated with the hydrodynamic pressure produced by the reservoir. In addition, ü_e and ü_g are the relative nodal acceleration and nodal ground acceleration vectors, respectively. The term Q is referred to as the coupling matrix. The discretization of hydrodynamic pressure equation to get the matrix form of the wave equation:

(2){{\displaystyle \mathit{M}}}_{\mathit{f}}{{\displaystyle \ddot{\mathit{P}}}}_{\mathit{e}}+{{\displaystyle \mathit{C}}}_{\mathit{f}}{{\displaystyle \dot{\mathit{P}}}}_{\mathit{e}}+{{\displaystyle \mathit{K}}}_{\mathit{f}}{{\displaystyle \mathit{P}}}_{\mathit{e}}+{{\displaystyle \rho }}_{\mathit{W}}{\mathit{Q}}^{\mathbf{T}}\left({{\displaystyle \ddot{u}}}_{\mathit{e}}+{{\displaystyle \ddot{\mathit{u}}}}_{\mathit{g}}\right)=\mathrm{0},

where M_f, C_f, and K_f are the fluid mass, damping and stiffness matrices, respectively and P_eü_e and ü_g are the nodal pressure, relative nodal acceleration and nodal ground acceleration vectors, respectively. The term is also referred to as the transpose of the coupling matrix. The dot represents the time derivative. Eqs. (1) and (2) describe the complete finite-element discretized equations for the dam-water-foundation rock interaction problem and can be written in an assembled form:

(3)\left[\begin{array}{l}{{\displaystyle \mathit{M}}}_{\mathit{s}}\\ {{\displaystyle \mathit{M}}}_{\mathit{f}\mathit{s}}\end{array}\begin{array}{l}\mathrm{0}\\ {{\displaystyle \mathit{M}}}_{\mathit{f}}\end{array}\right]\left\{\begin{array}{l}{{\displaystyle \ddot{\mathit{u}}}}_{\mathit{e}}\\ {{\displaystyle \ddot{\mathit{P}}}}_{\mathit{e}}\end{array}\right\}+\left[\begin{array}{l}{{\displaystyle \mathit{C}}}_{\mathit{s}}\\ \mathrm{0}\end{array}\begin{array}{l}\mathrm{0}\\ {{\displaystyle C}}_{\mathit{f}}\end{array}\right]\left\{\begin{array}{l}{{\displaystyle \dot{u}}}_{\mathit{e}}\\ {{\displaystyle \dot{\mathit{P}}}}_{\mathit{e}}\end{array}\right\}+\left\{\begin{array}{l}{{\displaystyle \mathit{u}}}_{\mathit{e}}\\ {{\displaystyle \mathit{P}}}_{\mathit{e}}\end{array}\right\}=\left\{\begin{array}{l}-{{\displaystyle \mathit{M}}}_{\mathit{s}}{{\displaystyle \ddot{\mathit{u}}}}_{\mathit{g}}\\ -{{\displaystyle \mathit{M}}}_{\mathit{f}\mathit{s}}{{\displaystyle \ddot{\mathit{u}}}}_{\mathit{g}}\end{array}\right\},

where K_fs = −Q and M_fs = {{\displaystyle \mathit{\rho }}}_{\mathit{W}}{\mathit{Q}}^{\mathit{T}}. Eq. (3) expresses a second order linear differential equation having unsymmetrical matrices and can be solved by means of direct integration methods [¹,¹⁰⁷,¹¹⁹].

Finite element model

A model of the Koyna concrete gravity dam, foundation and the reservoir is created in ABAQUS using 2D model. The dimensions of the global model is as follows:

1- The reservoir is (206*91.75) m²;

2- The dam is (103*70.2) m²;

3- The foundation is (379.2*103) m².

The global model is assigned with massed boundary condition type which consists of roller supports at the base in addition to dashpots in the sides of the soil to count for the effect of the transmitting boundaries or approximate local non reflecting boundary conditions (NRBC) to simulate the semi- infinite soil region to absorb the energy of the traveling waves which was proposed by Lysmer and Kuhlemeyer (1969) see Fig. 1.

The concrete dam part has been meshed using the standard linear plane stress (CPS4R) A4-node quadrilateral, reduced integration hourglass control elemnts. The water has been meshed using the explicit linear plain strain (CPE4R) A4-node quadrilateral, reduced integration, hourglass control elements with element deletion option. The foundation has been meshed using the standard linear plain strain (CPE4R) A4-node quadrilateral, reduced integration, hourglass control elemnts.

A surface to surface tangential behavior contact interactions have been used between the concrete dam and the water reservoir, between the water reservoir and the foundation and finally between the concrete dam and the foundation simultaneously but in the first the tangential type is frictionless and the second and the latter are rough type contact interactions. A self contact interaction was assigne to the left side of the water reservoir with frictionless tangential type. The material properties of the parts which assembles the global model are mainly the mass density, elastic and damping for both the concrete dam and the foundation, but fro the reservoir water only mass density was used. The viscosity and US-UP equation of state were used for the water and it is worthy to mention that the concrete damage plasticity and Mohr coulomb plasticity were used for the concrete dam and the foundation simultaneously.

For the non-coupled model the hydrostatic pressure of the water was assigned as shown in Eq. (4).

(4){{\displaystyle P}}_{hs}=\mathrm{0.5}{{\displaystyle \gamma }}_w{Z}^{\mathrm{2}},

where {{\displaystyle \gamma }}_w is the specific weight of the water and Z is the height of the reservoir.

The uplift pressure of the of the soil was assigned to both the coupled and non-coupled models as shown in Eq. (5).

(5){{\displaystyle P}}_u=\mathrm{0.5}{{\displaystyle \gamma }}_{w}ZB,

where B is the width of the dam base.

A ground motion acceleration with vertical and horizontal amplitudes from Koyna dam earthquake time history were applied in a dynamic explicit step of 10 seconds to simulate the coupling of hydrodynamic water pressure and the earthquake. The XFEM approach has been used to detect the cracks and damages in the concrete dam due to the earthquake ground motion supporting on the maximum principal stress in the concrete.

Mesh convergence

The global model of the Koyna concrete gravity dam has been prepared after performing mesh convergence analysis to obtain the most efficient model that is not suffering from the change of the outputs due to further mesh refinement. The criteria for the output was both the horizontal displacement of the base and the vertical differential settlement of the dam base due to coupling of the ground motion originated from the Koyna earthquake and the hydrodynamic pressure of the reservoir water mass for the duration of 10 seconds. Non-uniform mesh convergence was adopted for the assembly by considering 13 element seeding size cases starting from 626 to 1364 total elements. The selected total elements case was 1170 elements where the results are stable at all cases almost for both the horizontal displacement and the vertical differential settlement of the dam base see (Fig. 2).

Earthquake amplitudes

The Koyna earthquake occurred on December 11, 1967 near Koyna Dam in India. A real time history of the Koyna earthquake horizontal and vertical amplitudes are used for the duration of 10 seconds inthe simulation of coupled hydrodynamic and ground motion for the global model of the Koyna dam see both (Fig. 3) and (Fig. 4).

Coupling effects

When the coupling effect is considered in the simulations of the global model of the dam, it is important to compare and analyze the results for both non coupled and coupled models so that to determine the effects of this feature on the stability and the damage detection in the dam.

Results of stability

Considering the stability of the dam for horizontal displacement and differential settlement of the dam base, it is obvious that after the coupling the displacement of the crest is decreasing after 4 seconds from the coupled ground motion simulation and being stable till the end of the simulation where the decrease is 0.35 m in the result which is considered too high compared to the non coupled model. It is an indication that the non coupled model is not manifesting the actual behavior of the dam during the earthquake and the response would be overestimated see (Fig. 5).

Regarding the results of the differential settlement of the dam base. The response of the dam in both cases is the same almost, but after 3 seconds from the ground motion in the coupled case, a decrease in the differential settlement of the dam base can be seen by a range starting from 0.05 m to 0.15 m compared to the case of non coupled case. In the same way, the response of the dam in the non coupled case would be overestimated too, see (Fig. 6).

Results of damage detection

The effect of coupling in the earthquake simulation of the global model of the dam is shown in 10 screen shots in (Fig. 7) and the non coupled model is shown also in 10 screen shots for comparison during the earthquake. In the coupled model at 0.1 s time from the earthquake, the damaged regions are clear in the downstream part of the dam at the base end and the neck, but in the non coupled model the damage is at the upstream part at the end of the dam base only. At 0.4 s the damage region in the coupled model at the neck disappears and appears in the upstream near the neck and the damage in the downstream end of the dam base increases, while in the non coupled the damaged region will disappear from the end base and appear in the downstream neck region. The situation for damage in coupled model continues in the same region resting at the downstream part of of the dam end base and the damage rate increases and decreases till the end of the earthquake, while for the non coupled model the damage region changes from the neck to the base ends with increase and decrease in rate resting at the downstream part of the dam base end. This is an indication that there is a vital difference between the two models , and the damaged regions are commonly at the neck and the base regions of the dam, but the rate of damage and cracks propagation are widely different in both cases.

Uncertainty quantification

The uncertainty in the stability and damage detection in concrete gravity dams exposed to the coupling effect of both the earthquake and the hydrodynamic pressure of the reservoir is a critical status.The main effects and the interaction effects of many uncertain parameters are very important to quantify where many damage and failure situations in dams are still need to be studied. Using experimental methods with the help of computer simulations of such problems are essential factors for constructing the meta models. The meta models are efficient tools for uncertainty quantifications of the system response by predicting the global and local behaviors of the system. The contribution of [¹⁵] in the area of uncertainty quantification is a firm background for future analyses. In this study, the Box-Behnken experimental method is utilized in addition to Sobol’s sensitivity indices to quantify the uncertainties of five uncertain parametrs on the prediction of the stability and damage detection in concrete gravity dams.

Box-Behnken experimental method

Box–Behnken experimental design was based on to identify the relationship between the response functions (lateral displacement and differential settlement) of the concrete gravity dam model and the five design variables of the reservoir (reservoir height) and soil properties (density, modulus of elasticity, uplift pressure and beta coefficient of viscous damping). Box–Behnken design is rotatable second-order designs based on three-level incomplete factorial designs. The special arrangement of the Box-Behnken design levels allows the number of design points to increase at the same rate as the number of polynomial coefficients. For three factors, for example, the design can be constructed as three blocks of four experiments consisting of a full two-factor factorial design with the level of the third factor set at zero.

Box–Behnken design requires an experiment number according to N=k²+ k+ c_p, where (k) is the factor number and (c_p ) is the replicate number of the central point. Box-Behnken is a spherical, revolving design. Viewed as a cube (Fig. 8(a)), it consists of a central point and the middle points of the edges. However, it can also be viewed as consisting of three interlocking 2² factorial design and a central point (Fig. 8(b)). It has been applied for optimization of several chemical and physical processes.

For the three-level three-factorial Box–Behnken experimental design, a total of 50 experimental runs are needed. Table 1 shows the levels of variables for the experimental design. The meta model is of the following form:

(6)\begin{array}{l}y={{\displaystyle \beta }}_{\mathrm{0}}+{{\displaystyle \beta }}_{\mathrm{1}}{{\displaystyle x}}_{\mathrm{1}}+{{\displaystyle \beta }}_{\mathrm{2}}{{\displaystyle x}}_{\mathrm{2}}+{{\displaystyle \beta }}_{\mathrm{3}}{{\displaystyle x}}_{\mathrm{3}}+{{\displaystyle \beta }}_{\mathrm{4}}{{\displaystyle x}}_{\mathrm{4}}+{{\displaystyle \beta }}_{\mathrm{5}}{{\displaystyle x}}_{\mathrm{5}}+{{\displaystyle \beta }}_{\mathrm{11}}{{\displaystyle x^{\mathrm{2}}}}_{\mathrm{1}}+{{\displaystyle \beta }}_{\mathrm{22}}{{\displaystyle x^{\mathrm{2}}}}_{\mathrm{2}}+{{\displaystyle \beta }}_{\mathrm{33}}{{\displaystyle x^{\mathrm{2}}}}_{\mathrm{3}}+{{\displaystyle \beta }}_{\mathrm{44}}{{\displaystyle x^{\mathrm{2}}}}_{\mathrm{4}}+{{\displaystyle \beta }}_{\mathrm{55}}{{\displaystyle x^{\mathrm{2}}}}_{\mathrm{5}}\\ +{{\displaystyle \beta }}_{\mathrm{12}}{{\displaystyle x}}_{\mathrm{1}}{{\displaystyle x}}_{\mathrm{2}}+{{\displaystyle \beta }}_{\mathrm{13}}{{\displaystyle x}}_{\mathrm{1}}{{\displaystyle x}}_{\mathrm{3}}+{{\displaystyle \beta }}_{\mathrm{14}}{{\displaystyle x}}_{\mathrm{1}}{{\displaystyle x}}_{\mathrm{4}}+{{\displaystyle \beta }}_{\mathrm{15}}{{\displaystyle x}}_{\mathrm{1}}{{\displaystyle x}}_{\mathrm{5}}+{{\displaystyle \beta }}_{\mathrm{23}}{{\displaystyle x}}_{\mathrm{2}}{{\displaystyle x}}_{\mathrm{3}}+{{\displaystyle \beta }}_{\mathrm{24}}{{\displaystyle x}}_{\mathrm{2}}{{\displaystyle x}}_{\mathrm{4}}\\ +{{\displaystyle \beta }}_{\mathrm{25}}{{\displaystyle x}}_{\mathrm{2}}{{\displaystyle x}}_{\mathrm{5}}+{{\displaystyle \beta }}_{\mathrm{34}}{{\displaystyle x}}_{\mathrm{3}}{{\displaystyle x}}_{\mathrm{4}}+{{\displaystyle \beta }}_{\mathrm{35}}{{\displaystyle x}}_{\mathrm{3}}{{\displaystyle x}}_{\mathrm{5}}+{{\displaystyle \beta }}_{\mathrm{45}}{{\displaystyle x}}_{\mathrm{4}}{{\displaystyle x}}_{\mathrm{5}},\end{array}

where y is the predicted response, β₀ model constant; x₁, x₂, x₃, x₄ and x₅ are independent variables; β₁, β₂β₃, β₄ and β₅ are linear coefficients; β₁₂, β₁₃, β₁₄β₁₅ β₂₃ β₂₄ β₂₅ β₃₄ β₃₅and β₄₅ are cross product coefficients and β₁₁, β₂₂, β₃₃β₄₄and β₅₅ are the quadratic coefficients. The main effect (β_i) and two factors interactions (βij) have been estimated from the experimental results by computer simulation programming applying least square method using MATLAB codes [¹²⁷–¹³²].

Meta models results

The regression coefficients have been calculated from the second-degree meta model both for the displacement of the dam crest and the differential settlement of the dam base, and they were obtained using Box-Behnken experimental method. Total of 50 runs used to formulate the meta models based on quadratic orders. The meta model for the case of displacement of the dam crest has a coefficient of determination R² of 99.21% (see Fig. 9) which is an excellent representation of the predicted coupled model response, and just 0.79% of the system response still unexplained.

The meta model for the case of the differential settlement of the dam base has a coefficient of determination R² of 99.67% (see Fig. 10) which is also an excellent representation for the prediction of the coupled model response, and just 0.33% of the system response still unexplained. According to the good coefficient of determination we did not test further meta models.

Convergence of the results

A global sensitivity analysis supporting on Sobol's sensitivity indices requires certain or adequate samples of experiments to find out the predicted effect of the parameters on the response of the coupled model system. The most efficient number of samples is being identified through the convergence of the sum of first orders and total sensitivity indices of the parameters. All the sensitivity indices (first orders, interaction orders and total orders) for each parameter have been calculated using MATLAB codes. Two outputs have been utilized in the process of convergence which are the displacement of the dam crest and the differential settlement of the dam base. Fig. 11 and Fig. 12 show the relation between the number of samples and the sum of first orders sensitivity indices, in the same time between the number of samples and the sum of total sensitivity indices of all parameters. For the case of displacement of the dam crest (see Fig. 11), the two curves of the sum of first orders and total orders of sensitivity indices at the beginning are not coinciding to reach convergence till 150 samples. After this stage the two curves are starting to converge at the 150 number of samples, where the two curves continue to remain in a stable position after many times of changing the number of samples.

While for the differential settlement of the dam base (see Fig. 12) the two curves are trying to reach convergence where it starts at 200 number of samples and the process continues in the stability pattern. The convergence results of the two cases necessitate utilizing 150 samples of design experiments to efficiently get the predicted rational effects of each parameter on the variation of the displacement of the dam crest, and to consider 200 samples to obtain the predicted rational effects of the design parameters on the variation of the differential settlement of the dam base.

Sensitivity indices results and discussion

The main orders of sensitivity indices for each parameter in addition to their interaction orders were calculated considering the convergence process recommending the use of 150 samples to calculate the displacement of the dam crest (see Fig. 11). Furthermore, the main and interaction orders of sensitivity indices for each parameter were calculated supporting on 200 samples to calculate the differential settlement of the dam base (see Fig. 12). Supporting on the calculated results, the total sensitivity indices for each parameter have been calculated (see Table 2).

In relation with the displacement of the dam crest, the total order sensitivity index of the parameter reservoir height X₁ is 0.7845, this value is much bigger than the total order sensitivity index of the parameter density of soil X₂ which is 0.0089, also it is much bigger than the total order sensitivity index of the parameter soil modulus of elesticity X₃ which is 0.0770, also much bigger than the parameter uplift pressure of soil X₄ which is 0.0108 and finally much bigger than the parameter soil beta viscous damping coefficient X₅ which is 0.1407, this means that the displacement of the dam crest is 78.45% due to variation in the reservoir height, and it is 0.89% due to the variation in soil density, also it is 7.7% due to the variation in the soil modulus of elasticity, and it is 1.08% due to the variation in soil uplift pressure, also it is 14.07% due to the variation in the soil beta viscous damping coefficient. While the interaction index between the design parameter X₁ and X₂ is 0.0002 and between X₁ and X₃ is 0.0016, and between X₁ and X₄ is 0.0002 and between X₁ and X₅ is 0.0151 , while between X₂ and X₃ is 0.0002, and between X₂ and X₄ is 0.0001, and between X₂ and X₅ is 0.0004, and between X₃ and X₄ is 0.0002, and between X₃ and X₅ is 0.0012, and between X₄ and X₅ is 0.0002, which means that there is a small interaction between the input parameters taking part in the variation of the displacement of the dam crest especially the interaction between the reservoir height and the soil beta viscous damping coefficient. While considering the differential settlement of the dam base, the total order sensitivity index of the parameter reservoir height X₁ is 0.3433, this value is much bigger than the total order sensitivity index of the parameter density of soil X₂ which is 0.0531, also much bigger than the parameter uplift pressure of soil X₄ which is 0.1142 and finally much bigger than the parameter soil beta viscous damping coefficient X₅ which is 0.1179, but it smaller a little bit than the total order sensitivity index of the parameter soil modulus of elesticity X₃ which is 0.3881, this means that the differential settlement of the dam base is 34.33% due to variation in the reservoir height, and it is 5.31% due to the variation in soil density, also it is 38.81% due to the variation in the soil modulus of elasticity, and it is 11.42% due to the variation in soil uplift pressure, also it is 11.79% due to the variation in the soil beta viscous damping coefficient. While the interaction index between the design parameter X₁ and X₂ is 0.0001 and between X₁ and X₃ is 0.0028, and between X₁ and X₄ is 0.0005 and between X₁ and X₅ is 0.0093 , while between X₂ and X₃ is 0.0002, and between X₂ and X₄ is 0.0001, and between X₂ and X₅ is 0.0001, and between X₃ and X₄ is 0.0005, and between X₃ and X₅ is 0.0008, and between X₄ and X₅ is 0.0007, which means that there is a small interaction between the input parameters taking part in the variation of the differential settlement of the dam base especially the interaction between the reservoir height and the soil beta viscous damping coefficient.

Effects of parameters

In the following sections the effect of each parameter on the stability and damage detection of the Koyna dam models would be detailed with the support of XFEM approach and MATLAB code using the uncoupled models for simulating the effect of the earthquake on the response of the Koyna dam models, where the effect of the hydrostatic pressure is considered only and the foundation is clay.

Parameter X₁- results of stability

The effect of the reservoir height parameter on the displacement of the dam crest for both minimum and maximum values is shown in (Fig. 13). When the reservoir height is at the maximum value 91.75 m the displacement of the dam crest decreases upto 0.5 m compared to the case of minimum reservoir height 41.75 m. This indicates that when the reservoir is full the stability of the dam regarding horizontal displacement would be better achieved because the weight of the water counteracts the horizontal displacement of the dam.

Regarding the differential settlement of the dam base, when the maximum value of the reservoir height is assigned 91.75 m the differential settlement of the dam base is decreased 0.2 m compared to the minimum case of the reservoir height 41.75 m see (Fig. 14), which is an indication that the dam is better in safety situation against differential settlement when it is at maximum case because the weight of the water helps to prevent the side saturated clay to deform as a result better stability is provided.

Parameter X₁- results of damage detection

When the reservoir height is at minimum value 41.75 m, the cracks appear in the base of the dam near to the upstream part and propagate vertically toward the crest, see (Fig. 15(a)). While the reservoir height is at maximum value 91.75 m, the cracks appear at the same region as for previous case but they propagate to the upstream face of the dam see (Fig. 15(b)). This indicates that maximum hydrostatic pressure of the water creates larger tensile stresses in that region which helps the cracks to propagate toward the tension zone.

Parameter X₂ - results of stability

The effect of the soil density on the displacement of the dam crest for both minimum and maximum values is shown in (Fig. 16). When the soil density is at the maximum value 2200 kg/m³_, the displacement of the dam crest increases with a small rate of 0.025 m compared to the case of minimum soil density 1800 kg/m³ with anonlinear behavior . This indicates that variation of this parameter affects the stability of the dam regarding horizontal displacement of the dam crest a little bit because lower density soil type helps the dam base to settle totally (not differentially) more which supports decreasing the horizontal displacement of the dam crest.

Relating to the differential settlement of the dam base, when the maximum value of the soil density is assigned 2200 kg/m³ the differential settlement of the dam base is increased with a maximum value of 0.07 m compared to the minimum case of the soil density 1800 kg/m³ see (Fig. 17), which is an indication that the dam is better in safety situation against differential settlement when it is at minimum case because the low density of the caly makes it easier for the dam base to settle tptally before settling differentially.

Parameter X₂- results of damage detection

When the soil density is at minimum value 1800 kg/m3, the cracks appear in the base of the dam near to the upstream part and propagate toward the upstream face of the dam, see (Fig. 18(a)). In the case of the soil density at maximum value 2200 kg/m3, the cracks appear at the same region as for previous case and propagate to the upstream face of the dam too, see (Fig. 18(b)). This indicates that the variation of soil density does not affect the cracks appreciably and the way of its propagation , this is due to total settlement of the dam base which prevents the generation of additional appreciable stresses in all parts of the dam, see both (Fig. 18(a)) and (Fig. 18(b)).

Parameter X₃ - results of stability

The effect of the soil modulus of elasticity on the displacement of the dam crest for both minimum and maximum values is shown in (Fig. 19). When the soil modulus of elasticity is at the maximum value 2.5E+08 Pa, the displacement of the dam crest decreases upto 0.25 m compared to the case of minimum soil modulus of elasticity 2.0E+08 Pa. This indicates that variation of this parameter affects the stability of the dam regarding horizontal displacement of the dam crest because the foundation is stronger in the maximum case which prevents larger horizontal displacement of the dam.

Regarding the differential settlement of the dam base, when the maximum value of the soil modulus of elasticity is assigned 2.5E+08 Pa, the differential settlement of the dam base is decreased maximum 0.2 m compared to the minimum case of the soil modulus of elasticity 2.0E+08 Pa see (Fig. 20), which is an indication that the dam is better in safety situation against differential settlement when it is at maximum case because the foundation would be stronger and it prevents the dam base to settle larger as a result better stability is provided.

Parameter X₃- results of damage detection

When the soil modulus of elasticity is at minimum value 2.0E+08 Pa, the cracks appear in the base of the dam near to the upstream part and propagate toward the upstream face of the dam, see (Fig. 21(a)). In the case of the soil modulus of elasticity at maximum value 2.5E+08 Pa, the cracks appear at the same region as for previous case and propagate to the upstream face of the dam too but the length of the crack in this case is shorter than the previous case, see (Fig. 21(b)). This indicates that the variation of soil modulus of elasticity affects the cracks appreciably, this is because the stronger soil foundation prevents the generation of additional appreciable stresses in the cracked region of the dam toenhance the continue propagate cracks compared to the propagated cracks in the minimum case.

Parameter X₄ - results of stability

The effect of the soil uplift pressure on the displacement of the dam crest for both minimum and maximum values is shown in (Fig. 22). When the soil uplift pressure is at the maximum value 3.16E+07 Pa, the displacement of the dam crest decreases maximum of 0.15 m compared to the case of minimum soil uplift pressure 1.44E+07 Pa with anonlinear behavior . This indicates that variation of this parameter affects the stability of the dam regarding horizontal displacement of the dam crest with a fair rate because the uplift pressure of the soilcounteracts vertically the horizontal displacement of the dam which results in decreasing the horizontal displacement.

Regarding the differential settlement of the dam base, when the maximum value of the soil uplift pressure is assigned 3.16E+07 Pa, the differential settlement of the dam base is decreased maximum 0.1 m compared to the minimum case of the soil uplift pressure 1.44E+07 Pa, see (Fig. 23), which is an indication that the dam is better in safety situation against differential settlement when it is at maximum case because the soil uplift pressure counteracts the differential settlement of the dam base , thus it decreases the overall differential settlement.

Parameter X₄ - results of damage detection

When the soil uplift pressure is at minimum value 1.44E+07 Pa, the cracks appear in the base of the dam near to the downstream part and propagate with a slope manner toward the upstream dam face with showing a long path, see (Fig. 24(a)). While the soil uplift pressure is at maximum value 3.16E+07 Pa, the cracks appear at the middle part of the dam base but propagate toward the upstream face of the dam with a shorter length, see (Fig. 24(b)). This indicates that the maximum soil uplift pressure counteracts the differential settlement of the dam base and stabilizes it better , thus it decreases the tensile stresses in that region which helps the cracks to stop propagate anymore.

Parameter X₅ - results of stability

The effect of the soil beta viscous damping coefficient on the displacement of the dam crest for both minimum and maximum values is shown in (Fig. 25). When the soil beta viscous damping coefficient is at the maximum value 0.13, the displacement of the dam crest decreases upto maximum 0.25 m compared to the case of minimum soil beta viscous damping coefficient 0.03. This indicates that variation of this parameter affects the stability of the dam regarding horizontal displacement of the dam crest because the foundation damping accomodates the movements of the dam resulting from the ground motion.

Regarding the differential settlement of the dam base, when the maximum value of the soil beta viscous damping coefficient is assigned 0.13, the differential settlement of the dam base is decreased maximum 0.08 m compared to the minimum case of the soil beta viscous damping coefficient 0.03, see (Fig. 26), which is an indication that the dam is better in safety situation against differential settlement when it is at maximum case because the damping of the foundation accommodates the movement effect of the dam resulting from the ground motion, thus in resultant it decreases the overall differential settlement.

Parameter X₅ - results of damage detection

When the soil beta viscous damping coefficient is at minimum value 0.03, the cracks appear in the base of the dam near to the downstream part and propagate toward the upstream face of the dam (the uplift pressure parameter kept maximum) this why the pattern of the damage is similar to the previous maximum case for soil uplift pressure parameter, see (Fig. 27(a)). In the case of the soil beta viscous damping coefficient at maximum value 0.13, the cracks appear at the same region as for minimum case and propagate to the upstream face of the dam too with same pattern and length, see (Fig. 27(b)). This indicates that the variation of soil beta viscous damping coefficient does not affect the cracks in this case because the effect of maximum uplift pressure is overcoming where the interaction between the the reservoir height and the soil beta viscous damping coefficient has been determined in the sensitivity analysis section to have a small effect also.

Effect of foundation- results of damage detection

In this section the effect of the foundation type is studied regarding the damage detection in the dam by considering two effective parameters which are reservoir height and soil uplift pressure. It is worthy to mention that only these two parameters have been considered because they are directly and non directly not related to the soil parameters where the reservoir height parameter is directly related to the water in the upstream region and non directly the soil uplift pressure parameter which depends on the reservoir height in its action.

Parameter X₁

Related to using the minimum value of the reservoir height 41.75 m for both the clay and rock foundations see (Fig. 28(a)) and (Fig. 28(b)) and when considering the clay foundation the cracks appear in the base of the dam near to the upstream region and propagates toward the crest vertically with a recognized length relatively. But when the foundation is rock for the same conditions, the cracks appear at the middle of the dam base with a short length propagating in a zigzag style but the stresses increase because the rock foundation does not absorb the response of the dam like the clay foundation when exposed to the ground motion this by settlement action. When the maximum value of the reservoir height 91.75 m is used , in the clay foundation case the cracks appear in the dam base near the upstream region and propagate toward the dam crest with a certain length see(Fig. 28(c)). But when the foundation is rock for the same conditions, the cracks appear at the same region but propagate towards the near region of the upstream face of the dam near the base with a shorter length relatively, see (Fig. 28(d)). This behavior is due to absorbtion of movement of the dam by the clay foundation through settlement action with a resulted tensile stresses in the cracked region.

Parameter X₄

Related to using the minimum value of the soil uplift pressure 1.44E+07 Pa for both the clay and rock foundations see (Fig. 29(a)) and (Fig. 29(b)) and when considering the clay foundation the cracks appear in the base of the dam near to the downstream region and propagates toward the crest vertically with a recognized length relatively. But when the foundation is rock for the same conditions, the cracks appear at the middle of the dam base with a short length propagating towards the upstream face of the dam. This is because the rock foundation supports the dam from settling relatively compared to the clay foundation case and also the resulted tensile stresses are in a condition which does not enhance propagation of the cracks such as it was in the clay foundation. When the maximum value of the soil uplift pressure 3.16E+07 Pa is used , in the clay foundation case the cracks appear in the middle of the dam base and propagate horizontally toward the upstream region of the dam face with a short length see (Fig. 29(c)). But when the foundation is rock for the same conditions, the cracks appear at the dam base near the upstream face and propagate towards the left side of the dam base with a shorter length relatively, see (Fig. 29(d)). This behavior is due to absorbtion of movement of the dam by the clay foundation through settlement action also the uplift pressure counteracts the vertical movement and in the same time the resulted tensile stresses at this region of the dam is supporting this pattern of the cracks.

Conclusions

The following conclusions have been constructed.

1) The coupling feature between the ground motion and the hydrodynamic pressure of the reservoir water mass at the upstream region of the Koyna dam has been simulated and analyzed showing the sloshing of the water due to the earthquake which represents the actual behavior of the system regarding astability and the damage detection in the dam compared to the classic method used by the previous researchers by using the theoretical value of the hydrodynamic load without real coupling in the numerical analysis and simulations.

2) In the non coupled analysis , the responses of the dam regarding the horizontal displacement of the dam crest and the differential settlement of the dam base due to the earthquake are overestimated compared to the coupled analysis case, thus the design would be more complex and the cost of the execution would be much more. Also for the constructed dam structure the solutions would be not reliable and the the safety of the dam would be having a non logic big factor of safety.

3) The damage detection in the dam for the non coupled system would be not real because the global response of the dam would be overestimated, so the solutions also would be overestimated which makes the this process complex and costly and migh lead to wrong decisions regarding the rehabilitation of damaged dams due to earthquakes.

4) The meta models which represent the coupling effect of the dam system due to the earthquake showed a very good prediction of the dam system response, hence the controlling of the stability and failure of the dam in the design stage would be efficiently achieved.

5) The uncertainty of five parameters related to the reservoir and the foundation were quantified , and the rational effects of each parameter were determined regarding the stability and damage detection in the dam. The reservoir height parameter was the most important parameter that plays the biggest role in the variation of the horizontal displacement of the dam crest in addition to a small effect of interaction between the reservoir height and the beta viscous damping coefficient of the soil. While regarding the differential settlement of the dam base, the soil modulus of elasticity plays the biggest role in the variation of this response, also a little interaction between the reservoir height and the beta viscous damping coefficient of the soil.

6) The effects of reservoir height and uplift pressure of the foundation in both clay and rock foundations on the damage detection in the dam were simulated and analyzed, where the simulations showed that the rock foundation is much better for constructing concrete gravity dams because it helps to better resisting the propagation of cracks due to earthquakes. The minimum value of the reservoir height is behaving better regarding the resistance against crack propagation in the dam. Furthermore, the maximum uplift pressure of the foundation was the optimum situation to prevent the propagation of cracks in the dam.

References

[1] Ilinca C, Varvorea R, Popovici A. Influence of dynamic analysis methods on seismic response of a buttress dam. Mathematical Modelling in Civil Engineering, 2014, 10(3): 1–15

[2] Mohsin A Z, Omran HA, Al-Shukur A H K. Dynamic response of concrete gravity dam on random soil. International Journal of Civil Engineering and Technology, 2015, 6(11): 21–31

[3] Mehdipour B. Effect of foundation on seismic behavior of concrete dam considering the interaction of dam- Reservoir. Journal of Basic and Applied Scientific Research, 2013, 3(5): 13–20

[4] Wang G, Wang Y, Lu W, Zhou C, Chen M, Yan P. XFEM based seismic potential failure mode analysis of concrete gravity dam–water–foundation systems through incremental dynamic analysis. Engineering Structures, 2015, 98: 81–94

[5] Lu L, Li X, Zhou J, Chen G, Yun D. Numerical simulation of shock response and dynamic fracture of a concrete dam subjected to impact load. Earth Sciences Research Journal, 2016, 20(1): 1–6

[6] Huang J. Seismic Response Evaluation of Concrete Gravity Dams Subjected to Spatially Varying Earthquake Ground Motions. Dissertation for PhD degree. Drexel University, Philadelphia, USA, 2011

[7] Ghasemi H, Park H S, Rabczuk T. A level-set based IGA formulation for topology optimization of flexoelectric materials. Computer Methods in Applied Mechanics and Engineering, 2017, 313: 239–258

[8] Ghasemi H, Brighenti R, Zhuang X, Muthu J, Rabczuk T. Optimum fiber content and distribution in fiber-reinforced solids using a reliability and NURBS based sequential optimization approach. Structural and Multidisciplinary Optimization, 2015, 51(1): 99–112

[9] Zhang C, Nanthakumar S S, Lahmer T, Rabczuk T. Multiple cracks identification for piezoelectric structures. International Journal of Fracture, 2017, 206(2): 151–169

[10] Nanthakumar S, Zhuang X, Park H, Rabczuk T. Topology optimization of flexoelectric structures. Journal of the Mechanics and Physics of Solids, 2017, 105: 217–234

[11] Nanthakumar S, Lahmer T, Zhuang X, Park H S, Rabczuk T. Topology optimization of piezoelectric nanostructures. Journal of the Mechanics and Physics of Solids, 2016, 94: 316–335

[12] Nanthakumar S, Lahmer T, Zhuang X, Zi G, Rabczuk T. Detection of material interfaces using a regularized level set method in piezoelectric structures. Inverse Problems in Science and Engineering, 2016, 24(1): 153–176

[13] Nanthakumar S, Valizadeh N, Park H, Rabczuk T. Surface effects on shape and topology optimization of nanostructures. Computational Mechanics, 2015, 56(1): 97–112

[14] Nanthakumar S S, Lahmer T, Rabczuk T. Detection of flaws in piezoelectric structures using extended FEM. International Journal for Numerical Methods in Engineering, 2013, 96(6): 373–389

[15] 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

[16] 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

[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] Rabczuk T, Akkermann J, Eibl J. A numerical model for reinforced concrete structures. International Journal of Solids and Structures, 2005, 42(5-6): 1327–1354

[19] Bažant Z P. Why continuum damage is nonlocal: micromechanics arguments. Journal of Engineering Mechanics, 1991, 117(5): 1070–1087

[20] Thai T Q, Rabczuk T, Bazilevs Y, Meschke G. A higher-order stress-based gradient-enhanced damage model based on isogeometric analysis. Computer Methods in Applied Mechanics and Engineering, 2016, 304: 584–604

[21] Fleck N A, Hutchinson J W. A phenomenological theory for strain gradient effects in plasticity. Journal of the Mechanics and Physics of Solids, 1993, 41(12): 1825–1857

[22] Rabczuk T, Eibl J. Simulation of high velocity concrete fragmentation using SPH/MLSPH. International Journal for Numerical Methods in Engineering, 2003, 56(10): 1421–1444

[23] Rabczuk T, Eibl J, Stempniewski L. Numerical analysis of high speed concrete fragmentation using a meshfree Lagrangian method. Engineering Fracture Mechanics, 2004, 71(4–6): 547–556

[24] Rabczuk T, Xiao S P, Sauer M. Coupling of meshfree methods with nite elements: basic concepts and test results. Communications in Numerical Methods in Engineering, 2006, 22(10): 1031–1065

[25] Rabczuk T, Eibl J. Modelling dynamic failure of concrete with meshfree particle methods. International Journal of Impact Engineering, 2006, 32(11): 1878–1897

[26] Etse G, Willam K. Failure analysis of elastoviscoplastic material models. Journal of Engineering Mechanics, 1999, 125(1): 60–69

[27] Miehe C, Hofacker M, Welschinger F. A phase field model for rate-independent crack propagation: robust algorithmic implementation based on operator splits. Computer Methods in Applied Mechanics and Engineering, 2010, 199(45–48): 2765–2778

[28] Amiri F, Millan D, Arroyo M, Silani M, Rabczuk T. Fourth order phase-field model for local max-ent approximants applied to crack propagation. Computer Methods in Applied Mechanics and Engineering, 2016, 312(C): 254–275

[29] Areias P, Rabczuk T, Msekh M. Phase-field analysis of finite-strain plates and shells including element subdivision. Computer Methods in Applied Mechanics and Engineering, 2016, 312(C): 322–350

[30] Msekh M A, Silani M, Jamshidian M, Areias P, Zhuang X, Zi G, He P, Rabczuk T. Predictions of J integral and tensile strength of clay/epoxy nanocomposites material using phase-eld model. Composites. Part B, Engineering, 2016, 93: 97–114

[31] Hamdia K, Msekh M A, Silani M, Vu-Bac N, Zhuang X, Nguyen-Thoi T, Rabczuk T. Uncertainty quantication of the fracture properties of polymeric nanocomposites based on phase field modeling. Composite Structures, 2015, 133: 1177–1190

[32] Msekh M A, Sargado M, Jamshidian M, Areias P, Rabczuk T. ABAQUS implementation of phase-field model for brittle fracture. Computational Materials Science, 2015, 96: 472–484

[33] Amiri F, Millán D, Shen Y, Rabczuk T, Arroyo M. Phase-field modeling of fracture in linear thin shells. Theoretical and Applied Fracture Mechanics, 2014, 69: 102–109

[34] Hamdia K M, Zhuang X, He P, Rabczuk T. Fracture toughness of polymeric particle nanocomposites: evaluation of Models performance using Bayesian method. Composites Science and Technology, 2016, 126: 122–129

[35] Rabczuk T, Belytschko T, Xiao S P. Stable particle methods based on Lagrangian kernels. Computer Methods in Applied Mechanics and Engineering, 2004, 193(12–14): 1035–1063

[36] Rabczuk T, Belytschko T. Adaptivity for structured meshfree particle methods in 2D and 3D. International Journal for Numerical Methods in Engineering, 2005, 63(11): 1559–1582

[37] Nguyen V P, Rabczuk T, Bordas S, Duflot M. Meshless methods: a review and computer implementation aspects. Mathematics and Computers in Simulation, 2008, 79(3): 763–813

[38] Zhuang X, Cai Y, Augarde C. A meshless sub-region radial point interpolation method for accurate calculation of crack tip elds. Theoretical and Applied Fracture Mechanics, 2014, 69: 118–125

[39] Zhang X, Zhu H, Augarde C. An improved meshless Shepard and least square method possessing the delta property and requiring no singular weight function. Computational Mechanics, 2014, 53(2): 343–357

[40] Zhuang X, Augarde C, Mathisen K. Fracture modelling using meshless methods and level sets in 3D: framework and modelling. International Journal for Numerical Methods in Engineering, 2012, 92(11): 969–998

[41] Chen L, Rabczuk T, Bordas S, Liu G R, Zeng KY, Kerfriden P.Extended finite element method with edge-based strain smoothing (ESm-XFEM) for linear elastic crack growth. Computer Methods in Applied Mechanics and Engineering, 2012, 209‒212: 250–265

[42] Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. International Journal for Numerical Methods in Engineering, 1999, 45(5): 601–620

[43] Moës N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. International Journal for Numerical Methods in Engineering, 1999, 46(1): 131–150

[44] Vu-Bac N, Nguyen-Xuan H, Chen L, Lee C K, Zi G, Zhuang X, Liu G R, Rabczuk T. A phantom-node method with edge-based strain smoothing for linear elastic fracture mechanics. Journal of Applied Mathematics, 2013( 2013): 978026

[45] Bordas S P A, Natarajan S, Kerfriden P, Augarde C E, Mahapatra D R, Rabczuk T, Dal Pont S. On the performance of strain smoothing for quadratic and enriched finite element approximations (XFEM/GFEM/PUFEM). International Journal for Numerical Methods in Engineering, 2011, 86(4–5): 637–666

[46] Bordas S P A, Rabczuk T, Hung N X, Nguyen V P, Natarajan S, Bog T, Quan D M, Hiep N V. Strain Smoothing in FEM and XFEM. Computers & Structures, 2010, 88(23–24): 1419–1443

[47] Rabczuk T, Zi G, Gerstenberger A, Wall W A. A new crack tip element for the phantom node method with arbitrary cohesive cracks. International Journal for Numerical Methods in Engineering, 2008, 75(5): 577–599

[48] Chau-Dinh T, Zi G, Lee P S, RabczukT, Song J H. Phantom-node method for shell models with arbitrary cracks. Computers & Structures, 2012, 92– 93: 242–256

[49] Song J H, Areias P M A, Belytschko T. A method for dynamic crack and shear band propagation with phantom nodes. International Journal for Numerical Methods in Engineering, 2006, 67(6): 868–893

[50] Areias P M A, Song J H, Belytschko T. Analysis of fracture in thin shells by overlapping paired elements. Computer Methods in Applied Mechanics and Engineering, 2006, 195(41–43): 5343–5360

[51] Rabczuk T, Areias P M A. A meshfree thin shell for arbitrary evolving cracks based on an external enrichment. CMES-Computer Modeling in Engineering and Sciences, 2006, 16(2): 115–130

[52] Zi G, Rabczuk T, Wall W A. Extended meshfree methods without branch enrichment for cohesive cracks. Computational Mechanics, 2007, 40(2): 367–382

[53] Rabczuk T, Bordas S, Zi G. A three-dimensional meshfree method for continuous multiple-crack initiation, propagation and junction in statics and dynamics. Computational Mechanics, 2007, 40(3): 473–495

[54] Rabczuk T, Zi G. A meshfree method based on the local partition of unity for cohesive cracks. Computational Mechanics, 2007, 39(6): 743–760

[55] Rabczuk T, Areias P M A, Belytschko T. A meshfree thin shell method for nonlinear dynamic fracture. International Journal for Numerical Methods in Engineering, 2007, 72(5): 524–548

[56] Bordas S, Rabczuk T,Zi G. Three-dimensional crack initiation, propagation, branching and junction in non-linear materials by extrinsic discontinuous enrichment of meshfree methods without asymptotic enrichment. Engineering Fracture Mechanics, 2008, 75(5): 943–960

[57] Rabczuk T, Zi G, Bordas S, Nguyen-Xuan H. A geometrically non-linear three dimensional cohesive crack method for reinforced concrete structures. Engineering Fracture Mechanics, 2008, 75(16): 4740–4758

[58] Rabczuk T, Gracie R, Song J H, Belytschko T. Immersed particle method for fluid-structure interaction. International Journal for Numerical Methods in Engineering, 2010, 81(1): 48–71

[59] Rabczuk T, Bordas S, Zi G. On three-dimensional modelling of crack growth using partition of unity methods. Computers & Structures, 2010, 88(23–24): 1391–1411

[60] Amiri F, Anitescu C, Arroyo M, Bordas S, Rabczuk T. XLME interpolants, a seamless bridge between XFEM and enriched meshless methods. Computational Mechanics, 2014, 53(1): 45–57

[61] Talebi H, Samaniego C, Samaniego E, Rabczuk T. On the numerical stability and mass- lumping schemes for explicit enriched meshfree methods. International Journal for Numerical Methods in Engineering, 2012, 89(8): 1009–1027

[62] Nguyen-Thanh N, Zhou K, Zhuang X, Areias P, Nguyen-Xuan H, Bazilevs Y, Rabczuk T. Isogeometric analysis of large-deformation thin shells using RHTsplines for multiple-patch coupling. Computer Methods in Applied Mechanics and Engineering, 2017, 316: 1157–1178

[63] Nguyen-Thanh N, Valizadeh N, Nguyen M N, Nguyen-Xuan H, Zhuang X, Areias P, Zi G, Bazilevs Y, De Lorenzis L, Rabczuk T. An extended isogeometric thin shell analysis based on Kirchho-Love theory. Computer Methods in Applied Mechanics and Engineering, 2015, 284: 265–291

[64] Jia Y, Anitescu C, Ghorashi S, Rabczuk T. Extended isogeometric analysis for material interface problems. IMA Journal of Applied Mathematics, 2015, 80(3): 608–633

[65] Ghorashi S, Valizadeh N, Mohammadi S, Rabczuk T. T-spline based XIGA for fracture analysis of orthotropic media. Computers & Structures, 2015, 147: 138–146

[66] Rabczuk T, Belytschko T. Cracking particles: a simplified meshfree method for arbitrary evolving cracks. International Journal for Numerical Methods in Engineering, 2004, 61(13): 2316–2343

[67] Rabczuk T, Areias P M A. A new approach for modelling slip lines in geological materials with cohesive models. International Journal for Numerical and Analytical Methods in Geomechanics, 2006, 30(11): 1159–1172

[68] Rabczuk T, Belytschko T. Application of particle methods to static fracture of reinforced concrete structures. International Journal of Fracture, 2006, 137(1–4): 19–49

[69] Rabczuk T, Belytschko T. A three dimensional large deformation meshfree method for arbitrary evolving cracks. Computer Methods in Applied Mechanics and Engineering, 2007, 196(29–30): 2777–2799

[70] Rabczuk T, Areias P M A, Belytschko T. A simplied meshfree method for shear bands with cohesive surfaces. International Journal for Numerical Methods in Engineering, 2007, 69(5): 993–1021

[71] Rabczuk T, Samaniego E. Discontinuous modelling of shear bands using adaptive meshfree methods. Computer Methods in Applied Mechanics and Engineering, 2008, 197(6–8): 641–658

[72] Rabczuk T, Song J H, Belytschko T. Simulations of instability in dynamic fracture by the cracking particles method. Engineering Fracture Mechanics, 2009, 76(6): 730–741

[73] Rabczuk T, Zi G, Bordas S, Nguyen-Xuan H. A simple and robust three-dimensional cracking-particle method without enrichment. Computer Methods in Applied Mechanics and Engineering, 2010, 199(37–40): 2437–2455

[74] Cai Y,Zhuang X Y, Zhu H. A generalized and ecient method for nite cover generation in the numerical manifold method. International Journal of Computational Methods, 2013, 10(5): 1350028

[75] Liu G, Zhuang X, Cui Z. Three-dimensional slope stability analysis using independent cover based numerical manifold and vector method. Engineering Geology, 2017, 225: 83–95

[76] Nguyen H B, Zhuang X, Wriggers P, Rabczuk T, Mears M E, Tran H D. 3D Isogeometric symmetric Galerkin boundary element methods. Computer Methods in Applied Mechanics and Engineering, 2017, 323: 132–150

[77] Nguyen B H, Tran H D, Anitescu C, Zhuang X, Rabczuk T. An isogeometric symmetric Galerkin boundary element method for two-dimensional crack problems. Computer Methods in Applied Mechanics and Engineering, 2016, 306: 252–275

[78] Zhu H, Wu W, Chen J, Ma G, Liu X, Zhuang X. Integration of three dimensional discontinuous deformation analysis (DDA) with binocular photogrammetry for stability analysis of tunnels in blocky rock mass. Tunnelling and Underground Space Technology, 2016, 51: 30–40

[79] Wu W, Zhu H, Zhuang X, Ma G, Cai Y. A multi-shell cover algorithm for contact detection in the three dimensional discontinuous deformation analysis. Theoretical and Applied Fracture Mechanics, 2014, 72: 136–149

[80] Cai Y, Zhu H, Zhuang X. A continuous/discontinuous deformation analysis (CDDA) method based on deformable blocks for fracture modelling. Frontiers of Structural and Civil Engineering, 2013, 7(4): 369–378

[81] Nguyen-Xuan H, Liu G R, Bordas S, Natarajan S, Rabczuk T. An adaptive singular ES-FEM for mechanics problems with singular eld of arbitrary order. Computer Methods in Applied Mechanics and Engineering, 2013, 253: 252–273

[82] Areias P, Rabczuk T. Steiner-point free edge cutting of tetrahedral meshes with applications in fracture. Finite Elements in Analysis and Design, 2017, 132: 27–41

[83] Areias P, Reinoso J, Camanho P, Rabczuk T. A constitutive-based element-by-element crack propagation algorithm with local mesh refinement. Computational Mechanics, 2015, 56(2): 291–315

[84] Areias P M A, Rabczuk T, Camanho P P. Finite strain fracture of 2D problems with injected anisotropic softening elements. Theoretical and Applied Fracture Mechanics, 2014, 72: 50–63

[85] Areias P, Rabczuk T, Dias-da-Costa D. Element-wise fracture algorithm based on rotation of edges. Engineering Fracture Mechanics, 2013, 110: 113–137

[86] Areias P, Rabczuk T, Camanho P P. Initially rigid cohesive laws and fracture based on edge rotations. Computational Mechanics, 2013, 52(4): 931–947

[87] Areias P, Rabczuk T. Finite strain fracture of plates and shells with congurational forces and edge rotation. International Journal for Numerical Methods in Engineering, 2013, 94(12): 1099–1122

[88] Silani M, Talebi H, Hamouda A S, Rabczuk T. Nonlocal damage modeling in clay/epoxy nanocomposites using a multiscale approach. Journal of Computational Science, 2016, 15: 18–23

[89] Talebi H, Silani M, Rabczuk T. Concurrent multiscale modelling of three dimensional crack and dislocation propagation. Advances in Engineering Software, 2015, 80: 82–92

[90] Silani M, Talebi H, Ziaei-Rad S, Hamouda A M S, Zi G, Rabczuk T. A three dimensional extended Arlequin method for dynamic fracture. Computational Materials Science, 2015, 96: 425–431

[91] Silani M, Ziaei-Rad S, Talebi H, Rabczuk T. A semi-concurrent multiscale approach for modeling damage in nanocomposites. Theoretical and Applied Fracture Mechanics, 2014, 74: 30–38

[92] Talebi H, Silani M, Bordas S, Kerfriden P, Rabczuk T. A computational library for multiscale modelling of material failure. Computational Mechanics, 2014, 53(5): 1047–1071

[93] Talebi H, Silani M, Bordas S P A, Kerfriden P, Rabczuk T. Molecular dynamics/XFEM coupling by a three-dimensional extended bridging domain with applications to dynamic brittle fracture. International Journal for Multiscale Computational Engineering, 2013, 11(6): 527–541

[94] Yang S W, Budarapu P R, Mahapatra D R, Bordas S P A, Zi G, Rabczuk T. A meshless adaptive multiscale method for fracture. Computational Materials Science, 2015, 96: 382–395

[95] Budarapu P, Gracie R, Bordas S, Rabczuk T. An adaptive multiscale method for quasi-static crack growth. Computational Mechanics, 2014, 53(6): 1129–1148

[96] Budarapu P, Gracie R, Yang S W, Zhuang X, Rabczuk T. Ecient coarse graining in multiscale modeling of fracture. Theoretical and Applied Fracture Mechanics, 2014, 69: 126–143

[97] Zhuang X, Wang Q, Zhu H. Multiscale modelling of hydro-mechanical couplings in quasi-brittle materials. International Journal of Fracture, 2017, 204(1): 1–27

[98] Zhu H, Wang Q, Zhuang X. A nonlinear semi-concurrent multiscale method for fractures. International Journal of Impact Engineering, 2016, 87: 65–82

[99] Zhuang X, Wang Q, Zhu H. A 3D computational homogenization model for porous material and parameters identication. Computational Materials Science, 2015, 96: 536–548

[100] Kouznetsova V, Geers M G D, Brekelmans W A M. Multi-scale constitutive modelling of heterogeneous materials with a gradient-enhanced computational homogenization scheme. International Journal for Numerical Methods in Engineering, 2002, 54(8): 1235–1260

[101] Rabczuk T, Ren H. Peridynamic formulation for the modelling of quasi-static fractures and contacts in brittle rocks. Engineering Geology, 2017, 225: 42–48

[102] Amani J, Oterkus E, Areias P M A, Zi G, Nguyen-Thoi T, Rabczuk T. A non-ordinary state-based peridynamics formulation for thermoplastic fracture. International Journal of Impact Engineering, 2016, 87: 83–94

[103] Ren H, Zhuang X, Rabczuk T. A new peridynamic formulation with shear deformation for elastic solid. Journal of Micromechanics and Molecular Physics, 2016, 1(2): 1650009

[104] Ren H, Zhuang X, Cai Y, Rabczuk T. Dual-horizon peridynamics. International Journal for Numerical Methods in Engineering, 2016, 108(12): 1451–1476

[105] Ren H, Zhuang X, Rabczuk T. Dual-horizon peridynamics: a stable solution to varying horizons. Computer Methods in Applied Mechanics and Engineering, 2017, 318: 762–782

[106] Gaun F, Moore I D, Lin G. Seismic Analysis of Reservoir-Dam-Soil Systems in the Time Domain. The 8th international conference on Computer Methods and Advances in Geomechanics, Siriwardane & Zaman (Eds), 1994, 2: 917–922

[107] Zeidan B. State of Art in Design and Analysis of Concrete Gravity Dams. Tanta University, 2014, 1–56

[108] Wieland M, Ahlehagh S. Dynamic stability analysis of a gravity dam subject to the safety evaluation earthquake. 9th European club symposium IECS2013, Venice, Italy, 10–12 April, 2013

[109] Benerjee A, Paul D K, Dubey R N. Estimation of cumulative damage in concrete gravity dam considering foundatiom-reservoir interaction. Tenth U.S. National Conference on Earthquake Engineering Frontiers of Earthquake Engineering, Anchorage, Alaska, USA, 21–25 July, 2014

[110] Shariatmadar H, Mirhaj A. Modal Response of Dam-Reservoir-Foundation Interaction. 8th International Congress on Civil Engineering, Shiraz University, Shiraz, Iran, 11–13 May, 2009

[111] Berrabah A T, Belharizi M, Bekkouche A. Modal behavior of dam reservoir foundation system. Electronic Journal of Geotechnical Engineering, 2011, 16: 1593–1605

[112] Khosravi S, Heydari M. Modelling of concrete gravity dam including dam-water-foundation rock interaction. World Applied Sciences Journal, 2013, 22(4): 538–546

[113] Benerjee A, Paul D K, Dubey R N. Modelling issues in the seismic analysis of concrete gravity dams. Dam Engineering, 2014, 24(2): 18–38

[114] Joghataie A, Dizaji M S, Dizaji F S. Neural Network Software for Dam-Reservoir-Foundation Interaction. International Conference on Intelligent Computational Systems (ICICS'2012), Dubai, UAE, 7–8 January, 2012

[115] Hung J, Zerva A. Nonlinear analysis of a dam-reservoir-foundation system under spatially variable seismic excitations. The 14th World Conference on Earthquake Engineering, Beijing, China, 12–17 October, 2008

[116] Heirany Z, Ghaemian M. The effect of foundation’s modulus of elasticity on concrete gravity dam’s behavior. Indian Journal of Science and Technology, 2012, 5(5): 2738–2740

[117] Berrabah A T. Dynamic Soil-Fluid –Structure Interaction Applied for Concrete Dam. Dissertation for PhD degree. Chetouane, Algeria, Aboubekr Belkaid University, 2012

[118] Zeidan B. Effect of foundation flexibility on dam-reservoir-foundation interaction. Eighteenth International Water Technology Conference (IWTC18), Sharm El-Sheikh, Egypt, 12–14 March, 2015

[119] Burman A, Maity D, Sreedeep S. Iterative analysis of concrete gravity dam-nonlinear foundation interaction. International Journal of Engineering Science and Technology, 2010, 2(4): 85–99

[120] Chopra A K. Hydrodynamic pressures on dams during earthquakes. Journal of Engineering Mechanics, ASCE, 1967, 93(6): 205–223

[121] Fereydooni M, Jahanbakhsh M. An analysis of the hydrodynamic pressure on concrete gravity dams under earthquake forces using ANSYS software. Indian Journal of Fundamental and Applied Life Sciences, 2015, 5(S4): 984–988

[122] Berrabah AT, Armouti N, Belharizi M, Bekkouche A. Dynamic soil structure interaction study. Jordan Journal of Civil Engineering, 2012, 6(2): 161–173

[123] Pekan O A, Cui Y Z. Failure analysis of fractured dams during earthquakes by DEM. Eng Struct, 2004, 26: 1483–1502.

[124] Patil S V, Awari U R. Effect of soil structure interaction on gravity dam. International Journal of Science, Engineering and Technology Research (IJSETR), 2015, 4(4): 1046– 4053

[125] Chopra A K, Chakrabarti P. The earthquake experience at Koyna dam and stress in concrete gravity dam. Earthquake Engineering and Structural Dynamic, 1972, 1(2): 151–164

[126] Pal W. Seismic cracking of concrete gravity dam. Journal of Structural Division, 1976, 102 (ST9): 1827–1844

[127] Pasma S A, Daik R, Maskat M Y, Hassan O. Application of box-behnken design in optimization of glucose production from oil palm empty fruit bunch cellulose. International Journal of Polymer Science, 2013, 2013: 104502

[128] Qiu P, Cui M, Kang K, Park B, Son Y, Khim E, Jang M, Khim J. Application of Box–Behnken design with response surface methodology for modeling and optimizing ultrasonic oxidation of arsenite with H₂O₂. Central European Journal of Chemistry, 2014, 12(2): 164–172

[129] Tekindal M A, Bayrak H, Ozkaya B, Genc Y. Box Behnken experimental design in factorial experiments: the importance of bread for nutrition and health. Turkish Journal of Field Crops, 2012, 17(2): 115–123

[130] Amenaghawon N A, Nwaru K I, Aisien F A, Ogbeide S E, Okieimen C O. Application of Box-Behnken design for the optimization of citric acid production from corn starch using aspergillus niger. British Biotechnology Journal, 2013, 3(3): 236–245

[131] Souza A S, dos Santos W N L, Ferreira S L C. Application of Box–Behnken design in the optimization of an on-line pre-concentration system using knotted reactor for cadmium determination by flame atomic absorption spectrometry. Spectrochimica Acta Part B: Atomic Spectroscopy, 2005, 60(5): 737–742

[132] Aslan N, Cebeci Y. Application of Box-Behnken design and response surface methodology for modeling of some Turkish coals. Fuel, 2007, 86(1–2): 90–97

RIGHTS & PERMISSIONS

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