Dynamic failure analysis of concrete dams under air blast using coupled Euler-Lagrange finite element method

Abstract

In this study, the air blast response of the concrete dams including dam-reservoir interaction and acoustic cavitation in the reservoir is investigated. The finite element (FE) developed code are used to build three-dimensional (3D) finite element models of concrete dams. A fully coupled Euler-Lagrange formulation has been adopted herein. A previous developed model including the strain rate effects is employed to model the concrete material behavior subjected to blast loading. In addition, a one-fluid cavitating model is employed for the simulation of acoustic cavitation in the fluid domain. A parametric study is conducted to evaluate the effects of the air blast loading on the response of concrete dam systems. Hence, the analyses are performed for different heights of dam and different values of the charge distance from the charge center. Numerical results revealed that 1) concrete arch dams are more vulnerable to air blast loading than concrete gravity dams; 2) reservoir has mitigation effect on the response of concrete dams; 3) acoustic cavitation intensify crest displacement of concrete dams.

Keywords

air blast loading / concrete dams / finite element / dam-reservoir interaction / cavitation / concrete damage model

Introduction

Over the last two decades, the use of explosives by terrorist groups around the world that target civilian buildings and other structures is becoming a growing problem and has resulted in considerable attention which has been raised on the behavior of engineering structures under blast or impact loading.

However, blast response of many important civil infrastructures has not yet been well understood as a result of the complexities in their material behavior, loading and higher nonlinearities. Concrete dams which are indivisible parts of any society are an example of such important civil infrastructure used for storage of water, generation of electricity, etc. Blast loading on concrete dams may result in disaster due to water crisis, consequent flood flows and related damage events. Hence, understanding the dynamic behavior of concrete dams under blast loading through numerical simulations is of utmost importance. Therefore, the study of the effects of blast loading on structures has been a field of active research over the last decade [¹,²]. In the present research, different parameters affecting the nonlinear response of concrete dams and their mathematical representation are explained. The major contribution in this research are considering development of cavitating region in the reservoir that induced nonlinear behavior of water and implementing a strain rate dependent material law for concrete in the material nonlinearity analysis of dam body. Additionally, a detailed study is performed for the evaluation of performance of gravity and arch concrete dams. The dynamic response of the dams subjected to air blast loading is carried out for different dam height ranging from 50 to 150 m and 100 to 250 m for gravity and arch dams, respectively. The influence of dam height, standoff distance and reservoir cavitating on the performance of the concrete dams is also investigated. All numerical simulations are carried out based on the developed FORTRAN finite element code. A fully coupled numerical approach with combined Lagrangian and Eulerian methods, in which the reservoir are modeled using an Eulerian mesh, while the dam concrete is modeled using a Lagrangian mesh, is adopted to permit for the incorporation of the essential processes, namely the air blast loading, fluid-structure interaction, acoustic cavitation in the reservoir and nonlinear structural response. In the developed code, an iterative partitioned implicit scheme which is used to time integration of dynamic nonlinear equilibrium equations of fluid and structure domains and element-by-element preconditioned conjugate gradient (PCG) solver together with diagonal preconditioning are used to solve the large equation system resulting from the finite element discretization of the governing equations of fluid and structure domains. Figure 1 provides the sketch of the involved domains and boundary conditions of the finite element model. It is observed in Fig. 1 that the outer surface of the fluid domain is set to be a non-reflecting boundary.

Literature review

Concrete dam structures may experience air blast loading due to the attack of rockets and missiles. Study on the failure modes and performance of concrete dams subjected to air blast loading is essential to assess their safety. While the physical processes during an explosive detonated in air and shock wave propagation are extremely complex, the subsequent response of the dam subjected to explosion shock loading is much more complicated than that under other loading such as static and earthquake loadings. Although many researchers have conducted comprehensive experimental and numerical investigations related to model the blast effects on building structures [³–⁶], bridge structures [⁷–¹⁰], underground structures [¹¹–¹⁴], marine structures [¹⁵–¹⁷], and gravity dam structures [¹⁸,¹⁹], relatively little attention has been paid to dam structures under blast loads especially arch dams. This is possibly due to large size and the interactions with the reservoir and dam structure and the foundation which is responsible for the exorbitant price of both numerical modelling and experimental tests. In addition, the experimental tests require the use of relatively large amounts of charges which involves potential risks and need careful handling, which is typically not feasible in civilian research [²⁰]. In October 1940, the first in a series of experiments was conducted on a scale model of Mohne dam to determine if it could be destroyed by a huge conventional bomb. Further experiments involving a one-fiftieth scale model of Mohne dam and a full-size dam in mid-Wales demonstrated the Mohne dam could be destroyed if 6500 ib of high explosive detonated against the inner wall of the dam [²¹]. Xue et al. [²²] built a three-dimensional anisotropic dynamic damage model for Dagangshan arch dam subjected to blast loading. They concluded that the nonlinear behavior of concrete dams can be satisfactorily predicted using the proposed model. Some efficient numerical models for crack propagation under dynamic loading are developed by Rabczuk et al. [²³,²⁴]. In these studies, a meshfree scheme are used and also proposed a dynamic cohesive law for concrete that take into account the change of fracture energy and strain rate behavior under high dynamic loading conditions. An immersed particle method for fluid-structure interaction (FSI) problems under high velocity loading is presented in Ref. [²⁵], the proposed method is applied to predict the dynamic behavior of cracked structures that interact with adjacent fluid under explosive loading. However, few studies have focused their attention on cavitation in reservoir during blast loading and dynamic behavior of concrete dams especially arch dams under air blast loading. Acoustic cavitation can be formed in the reservoirs of dams during strong shaking of dam structures [²⁶]. The importance of cavitation in dam-water systems is as a result of its probable effects on the hydrodynamic forces acting on the dam which play an important role in the dam response. Previously developed cavitation model [²⁷] was used to study the effects of reservoir cavitation on the dam structures under air blast loading in the present research.

Air blast

Air blast which involves the detonation of explosive material in air can be modeled with a decaying exponential equation form that uses a series of parameters, and rely on explosive charge size, type and standoff point. Many different sets of air blast parameter data in both graphical and equation form are available in Refs. [²⁸–³³]. Thus, the objectives of this section are to describe blast parameters definitions from the open literature and implement these in a load generation code to produce air blast loading for finite element simulation of dam’s structures subjected to explosive air burst. The evolution of an air blast involves several stages: Detonation, shock wave formation and propagation and decay in the shock wave strength that terminate with a return to ambient conditions. Figure 2 illustrates a schematic diagram of a typical time history of overpressure at a stationary location affected by an explosion in air. This air blast profile comprises a positive phase and a negative phase. The negative phase of the air blast time history can be disregarded for structures in several situations [³⁰]. Any structure in the path of the shock wave will reflect the wave. The obstruction of the air velocity at the structure surface induces a significant increase in load on the structure. Thus, when describing the loading on a structure induced by explosive air burst, reflected and side-on cases are treated separately, with different values for reflected peak overpressure and side-on peak overpressure. At a stationary point in space, the effects of air blast have frequently been modeled with a modified Friedlander’s equation (Eq. (1)) [³⁴], as

(1)

where P_i\left(t\right) is the overpressure at time t after detonation, P_{\max } is the peak overpressure, t_{\mathrm{a}} is the arrival time of the shock wave, t_{\mathrm{d}} is the duration of the shock wave, and b is the decay constant. P_{\max } is either equal to P_{\mathrm{s}}, the peak side-on overpressure or P_r the maximum reflected overpressure depending on if the point of interest is located on the surface of an object or not. In empirical approach, air blast parameters are most often presented for a reference explosion and some type of scaling is subsequently used to obtain the parameter values for the actual charge weight of interest [²⁸–³¹,³³]. The Hopkinson scaling law states that when two charges of the same explosive material are detonated in the same atmospheric conditions, similar shock wave effects are experienced at equivalent scaled distances, Z, which is defined as

(2)Z={\displaystyle \frac{R}{\sqrt[\mathrm{3}]{W}}},

where R is the standoff distance, and W is the charge weight. Values for the parameters describing air blast in Eq. (1) can be found in a few different sources. The equations are taken from Ref. [³¹] and time duration is given as

(3){\displaystyle \frac{t_{\mathrm{d}}}{W^{\frac{\mathrm{1}}{\mathrm{3}}}}}={\displaystyle \frac{\mathrm{980}\left[\mathrm{1}+{\left({\displaystyle \frac{Z}{\mathrm{0.54}}}\right)}^{\mathrm{10}}\right]}{\left[\mathrm{1}+{\left({\displaystyle \frac{Z}{\mathrm{0.02}}}\right)}^{\mathrm{3}}\right]\left[\mathrm{1}+{\left({\displaystyle \frac{Z}{\mathrm{0.74}}}\right)}^{\mathrm{6}}\right]\sqrt{\mathrm{1}+{\left({\displaystyle \frac{Z}{\mathrm{6.9}}}\right)}^{\mathrm{2}}}}},

where t_{\mathrm{d}} is the duration of the positive phase of the blast profile (unit: s). Information about the peak overpressure in free air is also taken directly from Ref. [²⁸] and is defined as

P_{\mathrm{s}}=

(4){\displaystyle \frac{\mathrm{808}{P}_{\mathrm{atm}}\mathrm{1}+{\left({\displaystyle \frac{Z}{\mathrm{4.5}}}\right)}^{\mathrm{2}}}{\sqrt{\mathrm{1}+{\left({\displaystyle \frac{Z}{\mathrm{0.048}}}\right)}^{\mathrm{2}}\left[\mathrm{1}+{\left({\displaystyle \frac{Z}{\mathrm{0.32}}}\right)}^{\mathrm{2}}\right]\left[\mathrm{1}+{\left({\displaystyle \frac{Z}{\mathrm{1.35}}}\right)}^{\mathrm{2}}\right]}}},

where P_{\mathrm{s}} is the peak side-on overpressure in units of bars, and P_{\mathrm{atm}} is the atmospheric pressure (unit: bar, 1 bar=100 kPa). Information regarding peak reflected overpressure P_{\mathrm{r}} is much difficult to obtain than for incident overpressure in the open literature. All sources that do include parameter information for reflected overpressure, present data for the normally reflected case with angle of incident effects separately treated. In the far field limit for explosions of any size, or for small explosions, the air can be treated as an ideal gas in order to establish a relation between the peak side-on overpressure and peak reflected overpressure at a surface. According to Ref. [³⁵], this relation is

(5)

P_{\mathrm{r}}=P_{\mathrm{s}}(\mathrm{2}+{\displaystyle \frac{\mathrm{0.03851}{P}_{\mathrm{s}}}{\mathrm{1}+\mathrm{0.0025061}{P}_{\mathrm{s}}+\mathrm{4.041}\times {\mathrm{10}}^{-\mathrm{7}}{P}_{\mathrm{s}}^{\mathrm{2}}}}

(6)+{\displaystyle \frac{\mathrm{0.004218}+\mathrm{0.7011}{P}_{\mathrm{s}}+\mathrm{0.001442}{P}_{\mathrm{s}}^{\mathrm{2}}}{\mathrm{1}+\mathrm{0.1160}{P}_{\mathrm{s}}+\mathrm{8.086}\times {\mathrm{10}}^{-\mathrm{4}}{P}_{\mathrm{s}}^{\mathrm{2}}}},P_{\mathrm{s}}\ge \mathrm{6.9}\mathrm{bar},

where P_{\mathrm{r}} is the maximum overpressure for normal reflection, P_{\mathrm{s}} is the peak side-on overpressure, and P_{\mathrm{atm}} is the ambient air pressure. An implicit assumption in this equation is that \mathit{\gamma }=1.4, where \mathit{\gamma } is the heat capacity ratio of the air. When overpressure values exceed 6.9 bar, the molecules in the air interact with each other and the ideal gas assumption becomes invalid [³⁵]. Suitable expressions for the decay constant and arrival time are excluded from Refs. [³¹,³⁵]. Thus, in the present research is employed the proposed method in Ref. [²⁸] that fits piecewise polynomials to data for a 1 kg trinitrotoluene (TNT) reference explosion in Ref. [³¹]. The data includes arrival times and decay coefficients over a range of scaled distances. The resulting expression for arrival time is

(7){\displaystyle \frac{t_{\mathrm{a}}}{W^{\frac{\mathrm{1}}{\mathrm{3}}}}}={\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{4}}{a}_i{Z}^{i-\mathrm{1}},\mathrm{0.3}{\mathrm{m/kg}}^{\frac{\mathrm{1}}{\mathrm{3}}}\le Z\le \mathrm{500}{\mathrm{m/kg}}^{\frac{\mathrm{1}}{\mathrm{3}}}}.

Values for the fitted polynomial coefficients, a_i, are included in Table 1 [²⁸] for various ranges of Z.

Same as arrival time for the decay constant form in Ref. [²⁸] is given as

(8)b={\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{6}}{c}_i{Z}^{i-\mathrm{1}},\mathrm{0.3}{\mathrm{m/kg}}^{\frac{\mathrm{1}}{\mathrm{3}}}\le Z\le \mathrm{500}{\mathrm{m/kg}}^{\frac{\mathrm{1}}{\mathrm{3}}}},

where b is the dimensionless decay constant for side-on air blast. Values for the fitted polynomial coefficients, c_i, are illustrated for different range of Z in Table 2 [²⁸].

Structural equations

The second order equations of motion for solid can be written as

(9)\nabla ·\mathit{\sigma }+{\mathit{f}}_{\mathrm{s}}^{\mathrm{ext}}={\rho }_{\mathrm{s}}{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{u}}{\partial t^{\mathrm{2}}}},

where \mathit{u} corresponds to the displacement of the structure, \mathit{\sigma } is the structural stress tensor, {\mathit{\rho }}_{\mathrm{s}} refers to solid density, and {\mathit{f}}_{\mathrm{s}}^{\mathrm{ext}} is the load vector due to the external structural loads. The equations of motion for solid (Eq. (9)) is written in the most general form, which could include both material and geometric nonlinearities. Using standard procedures for finite element discretization of the structural domain, the equations of motion for the structure subjected to external blast forces may be written in standard finite element form as

(10)\mathit{M}\ddot{\mathit{u}}+\mathit{C}\dot{\mathit{u}}+{\displaystyle \int {\mathit{B}}^{\mathrm{T}}\mathit{\sigma }}\mathrm{d}\mathrm{\Omega }={\mathit{F}}_{\mathrm{a}}+{\mathit{F}}_{\mathrm{s}},

where \mathit{M} and \mathit{C} are the structural mass matrix and damping matrix, respectively, \mathit{B} is displacement-strain related matrix, \mathit{\sigma } is the tensor of internal stresses of structure, \mathit{u} is the vector of nodal displacements relative to the ground, and and {\mathit{F}}_{\mathrm{s}} is the vector of forces associated with the air blast loading and hydrodynamic pressure produced by fluid domain, respectively. Rayleigh damping is assumed in the present work and thus the global damping matrix is computed using the expression:

(11)\mathit{C}=\alpha \mathit{M}+\beta \mathit{K},

where K is the initial stiffness matrix of dam structure, \alpha and \beta are proportionality constants selected to control the damping ratios of the lowest and highest modes which are expected to contribute significantly to the response. These can be calculated from the relation

(12)\alpha +\beta {\omega }_i^{\mathrm{2}}=\mathrm{2}{\omega }_i{\xi }_i,

where {\xi }_i are the damping ratios and {\omega }_iare the ith natural frequency of the system. The main disadvantage of the Rayleigh damping method is that the higher modes are considerably more damped than lower modes, and that the damping can be controlled at only two modes of vibration. In various practical structural problems, the mass damping may be ignored and then the structural damping can be calculated as [³⁶]:

(13)\mathit{C}={\displaystyle \frac{\xi }{\mathrm{\pi }\omega }}\mathit{K},

where \omega is the main frequency of the structure, and \xi is a damping ratio. In the current study \xi =\mathrm{0.1} is considered.

The effect of strain rate on the concrete compressive and tensile strength is typically represented by a parameter, known as the dynamic increase factor (DIF). It is a ratio of the dynamic to static material constants versus strain rate. Usually, the strain rate effect is also response time history dependent [³⁷,³⁸]. However, in practice it is always assumed to depend solely on the strain rate. Various empirical relations use to evaluate strain rate effect on concrete material properties are obtainable in Refs. [³⁹,⁴⁰]. Recently, meshfree methods due to their ability and accuracy in the numerical simulation of lager deformation and crack path tracking in the concrete material is gain considerable attention and efficient numerical meshfree methods is developed [⁴¹,⁴²]. In the present study the same procedure used in Ref. [⁴³] is implemented. To define failure surface of concrete material, a piece-wise Drucker-Prager strength criterion is employed in the present study (similar to the criterion used to model concrete slab under explosion by Zhou et al. [⁴⁴]). Heterogeneity and porosity of concrete resulted in its complex nonlinear compression behavior, for this reason, the porous equation of state proposed by Herrmann [⁴⁵] that is well adapted to capture the major macroscopic phenomena of concrete and also describe the relationship between hydrostatic pressure and volume is applied in the present research as the equation of state (EOS) of concrete material. The material parameters for ordinary concrete adopted in the present work are based on Ref. [⁴⁴], as shown in Table 3 [⁴⁶].

Fluid equations

The set of governing equations, which describes the fluid domain in Cartesian coordinates, are the Euler equations where the viscosity, thermal conductivity, surface tension and turbulence are generally ignored. As the liquid and cavitating fluid are assumed to be compressible and barotropic in this work, total energy equation is not required to be solved directly. These equations can be stated as follow.

The Continuity equation:

(14){\displaystyle \frac{\partial \rho }{\partial t}}+{\rho }_{\mathrm{0}}\nabla ·\mathit{v}=q.

The momentum equation:

(15){\rho }_{\mathrm{0}}\left({\displaystyle \frac{\partial \mathit{v}}{\partial t}}+\mathit{v}·\nabla \mathit{v}\right)=-\nabla P+{\rho }_{\mathrm{0}}b.

In Eqs. (14) and (15), \rho and {\rho }_{\mathrm{0}} are the fluid density and the reference density of the fluid, respectively [⁴⁷], pis the pressure, v is the fluid velocity vector with three components, b is the body force, and q is the added fluid mass per unit volume and time. The nonlinear convection term \mathit{v}\cdot \nabla \mathit{v} in the Euler equation can be neglected for the acoustic fluid when the fluid velocity is small compared to the dimensions of the model. Therefore, the linearized Euler equation takes the form:

(16){\displaystyle \frac{\partial \mathit{v}}{\partial \mathit{t}}}=-{\displaystyle \frac{\nabla P}{{\rho }_{\mathrm{0}}}}+b.

The equation of state will be described in Section 5.1:

(17)P=H\left(\rho \right).

Density of vapor and fluid mixture can be defined as

(18)\rho =\lambda {\rho }_{\mathrm{vapor}}+\left(\mathrm{1}-\lambda \right){\rho }_{\mathrm{liquid}},

where {\rho }_{\mathrm{vapor}} and {\rho }_{\mathrm{liquid}} are the vapor and liquid densities, respectively. The volume ratio \lambda characterizes the volume of vapor in each element: \lambda =\mathrm{1} means that the element is completely filled with vapor; conversely, a complete liquid element is represented by \lambda =\mathrm{0}. Therefore, fluid is considered as single-fluid, whose density varies from liquid density to vapor density according to a barotropic equation of state, based on the expression of the local sound speed depending on the local volume fraction \lambda. By differentiating Eq. (14) with respect to time and eliminating the velocity field we have [⁴⁶]

(19){\nabla }^{\mathrm{2}}P={\displaystyle \frac{{\partial }^{\mathrm{2}}\rho }{\partial t^{\mathrm{2}}}}+{\rho }_{\mathrm{0}}\nabla b-{\displaystyle \frac{\partial q}{\partial t}}.

For homogenous and thermal equilibrium of the liquid and vapor phases with negligible viscous dissipation in isentropic isothermal cavitating water, an energy statement can be expressed as

(20)\mathrm{d}P=c^{\mathrm{2}}\mathrm{d}\rho ,

where c is the speed of sound. Equation (20) represents cavitation as an inertially dominated process with negligible thermodynamic effects. This assumption is valid since the heat transfer from the vapor bubble to the surrounding liquid field is almost instantaneous [⁴⁸]. When cavitation occurs, the density in the flow equations is treated as a pseudo density, which is related to the vapor and liquid density by Eq. (18). Substituting Eq. (20) into Eq. (19) we have

(21){\nabla }^{\mathrm{2}}P={\displaystyle \frac{\mathrm{1}}{c{\left(\lambda P_{\mathrm{absolute}}\right)}^{\mathrm{2}}}}\hspace{0.17em}{\displaystyle \frac{{\partial }^{\mathrm{2}}P}{\partial t^{\mathrm{2}}}}+\rho \nabla b-{\displaystyle \frac{\partial q}{\partial t}},

where P_{\mathrm{absoulte}} is the absolute pressure of fluid.

Equation of state for cavitation mixture

The acoustic cavitation modeling method used in this study is the same method that developed in our previous study [²⁷,⁴⁹]. Since the formulation and derivation are described in detail in Ref. [²⁷], only the outline is given here. The one-fluid cavitation model is found to be a simple yet efficient EOS for the acoustic cavitation [⁵⁰]. The major characteristic of the one-fluid model is that all flow phases are treated as compressible. The model is mathematically and physically consistent, and the cavitation dynamics are captured automatically. For the simulation of cavitation which occurs in cold water and under high pressure condition, the effects of thermal conductivity and viscosity on cavitation evolution are usually very much smaller than the influence of the pressure driving force of the system [⁴⁶]. Therefore, flow viscosity, turbulence, thermal non-equilibrium and cavitation surface tension are usually ignored in the application of the one-fluid cavitation model. The one-fluid model is considered to be appropriate for the simulation of unsteady cavitation where the explicit phase exchange rates are replaced by Eq. (29) and the mass transfer is evaluated by updating the void fraction at each time step. In the present study, it is used an EOS based on the expression of acoustic speed of the two-phase flow formulated by Wallis [⁵¹] as

c(\lambda p_{\mathrm{absolute}})={\left[(\lambda {\rho }_{\mathrm{vapor}}+(\mathrm{1}-\lambda ){\rho }_{\mathrm{liquid}})\left({\displaystyle \frac{\lambda }{{\rho }_{\mathrm{vapor}}{c}_{\mathrm{vapor}}^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{1}-\lambda }{{\rho }_{\mathrm{liquid}}{c}_{\mathrm{liquid}}^{\mathrm{2}}}}\right)\right]}^{-\frac{\mathrm{1}}{\mathrm{2}}},where {\rho }_{\mathrm{vapor}} and {\rho }_{\mathrm{liquid}} in Eq. (22) are the density of vapor and liquid phases, respectively, and the values of c_{\mathrm{vapor}}^{} and c_{\mathrm{liquid}} are the speed of sound in pure vapor and liquid, respectively.

Finite element implementation of fluid domain

A weak formulation of Eq. (21) can be constructed by multiplying with a weight function, w, and integrating by part. The weak formulation can be written as follows:

{\displaystyle \underset{{\mathrm{\Omega }}_{\mathrm{f}}}{\int }{\displaystyle \frac{\mathrm{1}}{c^{\mathrm{2}}P}}}\ddot{P}w\mathrm{d}\mathrm{\Omega }+{\displaystyle \underset{{\Omega }_{\mathrm{f}}}{\int }\nabla P·\nabla w}\mathrm{d}\mathrm{\Omega }

(23)={\displaystyle \underset{{\mathrm{\Gamma }}_{\mathrm{f}}}{\int }{\displaystyle \frac{\partial P}{\partial n}}}w\mathrm{d}\mathrm{\Gamma }+{\displaystyle \underset{{\Omega }_{\mathrm{f}}}{\int }{\rho }_{\mathrm{f}}\nabla bw\mathrm{d}\mathrm{\Omega }}-{\displaystyle \underset{{\Omega }_{\mathrm{f}}}{\int }\dot{q}w\mathrm{d}\mathrm{\Omega }},

where, {\mathrm{\Omega }}_{\mathrm{f}} and {\mathrm{\Gamma }}_{\mathrm{f}} denote the fluid domain and its boundaries. With substituting fluid domain boundary conditions into Eq. (23) and representing the unknown by a finite element summation, P={\displaystyle \sum _{i=\mathrm{1}}^m{P}_i{N}_{\mathrm{f}i}=P^{\mathrm{T}}{N}_{\mathrm{f}}} and assuming w_i=N_i lead to the following set of nonlinear ordinary differential equations in time:

(24){\mathit{F}}_{\mathrm{f}}^{\mathrm{int}}\left(\ddot{P},\dot{P},P\right)={\mathit{F}}_{\mathrm{f}}^{\mathrm{e}\mathrm{x}\mathrm{t}}\left(\ddot{u},b,\dot{q}\right).

For an individual fluid element (e), these are given by the following expressions:

(25)\begin{array}{l}{\mathit{F}}_{\mathrm{f}}^{\mathrm{int}}{\left(\ddot{P},\dot{P},P\right)}_e=\left({\displaystyle \underset{\mathrm{\Omega }}{\int }{\displaystyle \frac{\mathrm{1}}{c^{\mathrm{2}}}}{N}_{\mathrm{f}}^{\mathrm{T}}{N}_{\mathrm{f}}\mathrm{d}\mathrm{\Omega }}+{\displaystyle \frac{\mathrm{1}}{g}}{\displaystyle \underset{{\mathrm{\Gamma }}_{\mathrm{6}}}{\int }{N}_{\mathrm{f}}^{\mathrm{T}}{N}_{\mathrm{f}}\mathrm{d}\mathrm{\Gamma }}\right){\ddot{\mathit{p}}}_e\hfill \\ \hfill \\ +({\displaystyle \underset{{\mathrm{\Gamma }}_1}{\int }{\displaystyle \frac{\mathrm{1}}{c{\beta }_{\mathrm{1}}}}{N}_{\mathrm{f}}^{\mathrm{T}}{N}_{\mathrm{f}}\mathrm{d}\mathrm{\Gamma }+{\displaystyle \underset{{\mathrm{\Gamma }}_{\mathrm{2}}}{\int }{\displaystyle \frac{\mathrm{1}}{c{\beta }_{\mathrm{2}}}}{N}_{\mathrm{f}}^{\mathrm{T}}{N}_{\mathrm{f}}\mathrm{d}\mathrm{\Gamma }}}\hfill \\ \hfill \\ {\displaystyle +{\displaystyle \underset{{\mathrm{\Gamma }}_{\mathrm{3}}}{\int }{\displaystyle \frac{\mathrm{1}}{c{\beta }_{\mathrm{3}}}}{N}_{\mathrm{f}}^{\mathrm{T}}{N}_{\mathrm{f}}\mathrm{d}\mathrm{\Gamma }}{\displaystyle +\underset{{\mathrm{\Gamma }}_{\mathrm{4}}}{\int }{\displaystyle \frac{\mathrm{1}}{c}}{N}_{\mathrm{f}}^{\mathrm{T}}{N}_{\mathrm{f}}\mathrm{d}\mathrm{\Gamma }}}{\dot{\mathit{P}}}_e\hfill \\ \hfill \\ +\left({\displaystyle \underset{\mathrm{\Omega }}{\int }\nabla N_{\mathrm{f}}^{\mathrm{T}}\nabla N_{\mathrm{f}}\mathrm{d}\mathrm{\Omega }+{\displaystyle \underset{{\mathrm{\Gamma }}_{\mathrm{4}}}{\int }{\displaystyle \frac{\mathrm{\pi }}{\mathrm{2}h}}{N}_{\mathrm{f}}^{\mathrm{T}}{N}_{\mathrm{f}}\mathrm{d}\mathrm{\Gamma }}}\right){\mathit{P}}_e\hfill & ,\end{array}

and

(26){\mathit{F}}_{\mathrm{ext}}^{\mathrm{int}}{\left(\ddot{\mathit{u}},b,\dot{q}\right)}_e=-{\rho }_{\mathrm{f}}{\displaystyle \underset{{\mathrm{\Gamma }}_{\mathrm{1}}}{\int }{N}_{\mathrm{f}}n{N}_{\mathrm{s}}^{\mathrm{T}}\mathrm{d}\mathrm{\Gamma }}{\ddot{\mathit{u}}}_e+{\mathit{F}}_{\mathrm{b}}{\left(b,\dot{q}\right)}_e,

where, F_{\mathrm{f}}^{\mathrm{int}} is the internal force of fluid element, which depends on P and its first two derivatives, F_{\mathrm{f}}^{\mathrm{ext}} is the external force, {\mathrm{\Gamma }}_{\mathrm{1}},{\mathrm{\Gamma }}_{\mathrm{2}},\mathrm{...},{\mathrm{\Gamma }}_{\mathrm{6}} are reservoir boundary surfaces, N_{\mathrm{s}} and N_{\mathrm{f}} are shape functions of structure and fluid domain, respectively, n is normal to the interface and its direction as well as its value at a point are constant, \ddot{\mathit{u}} is the nodal acceleration produced by the flexible structure, {\mathit{F}}_{\mathrm{b}}{\left(b,\dot{q}\right)}_e is body force that acts on eth fluid element, {\beta }_{\mathrm{1}},{\beta }_{\mathrm{2}} and {\beta }_{\mathrm{3}} are acoustic impedance coefficients of material in the bottom, right bank and left bank of reservoir, respectively. The matrix form of Eq. (24) can be written as follow

(27)\mathit{G}\left(\mathit{p}\right)·\ddot{\mathit{P}}+\mathit{D}\left(\mathit{P}\right)·\dot{\mathit{P}}+\mathit{H}·\mathit{P}=-\rho {\mathit{Q}}^{\mathrm{T}}\ddot{\mathit{u}}+{\mathit{F}}_{\mathrm{f}}.

where {\mathit{G}}_{ij}={\displaystyle \sum {\mathit{G}}_{ij}^e}, {\mathit{D}}_{ij}={\displaystyle \sum {\mathit{D}}_{ij}^e}, and {\mathit{H}}_{ij}={\displaystyle \sum {\mathit{H}}_{ij}^e} are matrices representing the mass, damping and stiffness matrix of the fluid domain, respectively. {\mathit{F}}_{\mathrm{f}}={\displaystyle \sum {\left\{{\mathit{F}}_{\mathrm{b}}\right\}}_i^e} is external body force of fluid domain and Q_{ij}={\displaystyle \sum {\mathit{Q}}_{ij}^e} is the total coupling matrix. {\ddot{\mathit{u}}}_{\mathrm{total}} is total acceleration along interface. The coefficient G_{ij}^e, D_{ij}^e, H_{ij}^e and F_i^e for eth fluid element may be defined as

(28a){\mathit{G}}_{ij}^e\left(\mathit{p}\right)={\displaystyle \underset{\mathrm{\Omega }}{\int }{\displaystyle \frac{\mathrm{1}}{c{\left(\mathit{P}\right)}^{\mathrm{2}}}}{\left\{N_i\right\}}_{\mathrm{f}}^{\mathrm{T}}{\left\{N_j\right\}}_{\mathrm{f}}\mathrm{d}\mathrm{\Omega }}+{\displaystyle \frac{\mathrm{1}}{g}}{\displaystyle \underset{{\mathrm{\Gamma }}_{\mathrm{6}}}{\int }{\left\{N_i\right\}}_{\mathrm{f}}^{\mathrm{T}}{\left\{N_j\right\}}_{\mathrm{f}}\mathrm{d}\mathrm{\Gamma }},

{\mathit{D}}_{ij}^e\left(\mathit{P}\right)=[{\displaystyle \underset{\mathrm{\Gamma }1}{\int }{\displaystyle \frac{\mathrm{1}}{\mathit{c}\left(\mathit{P}\right){\beta }_{\mathrm{1}}}}{\left\{{\mathit{N}}_i\right\}}_{\mathrm{f}}^{\mathrm{T}}{\left\{{\mathit{N}}_j\right\}}_{\mathrm{f}}\mathrm{d}\mathrm{\Gamma }}

{\displaystyle +{\displaystyle \underset{\mathrm{\Gamma }\mathrm{2}}{\int }{\displaystyle \frac{\mathrm{1}}{\mathit{c}\left(\mathit{P}\right){\beta }_{\mathrm{2}}}}{\left\{{\mathit{N}}_i\right\}}_{\mathrm{f}}^{\mathrm{T}}{\left\{{\mathit{N}}_j\right\}}_{\mathrm{f}}\mathrm{d}\mathrm{\Gamma }}}

{\displaystyle +\underset{\mathrm{\Gamma }\mathrm{3}}{\int }{\displaystyle \frac{\mathrm{1}}{\mathit{c}\left(\mathit{P}\right){\beta }_{\mathrm{3}}}}{\left\{{\mathit{N}}_i\right\}}_{\mathrm{f}}^{\mathrm{T}}{\left\{{\mathit{N}}_j\right\}}_{\mathrm{f}}\mathrm{d}\mathrm{\Gamma }}

(28b){\displaystyle +\underset{\mathrm{\Gamma }\mathrm{4}}{\int }{\displaystyle \frac{\mathrm{1}}{c\left(\mathit{P}\right)}}{\left\{{\mathit{N}}_i\right\}}_{\mathrm{f}}^{\mathrm{T}}{\left\{{\mathit{N}}_j\right\}}_{\mathrm{f}}\mathrm{d}\mathrm{\Gamma }},

(28{{\rm\tf="fp1"\char"0063}})H_{ij}^e={\displaystyle \underset{\mathrm{\Omega }}{\int }\nabla {\left\{N_i\right\}}_{\mathrm{f}}^{\mathrm{T}}\nabla {\left\{N_j\right\}}_{\mathrm{f}}\mathrm{d}\mathrm{\Omega }+{\displaystyle \underset{{\mathrm{\Gamma }}_{\mathrm{4}}}{\int }{\displaystyle \frac{\mathrm{\pi }}{\mathrm{2}h}}{\left\{N_i\right\}}_{\mathrm{f}}^{\mathrm{T}}{\left\{N_j\right\}}_{\mathrm{f}}\mathrm{d}\mathrm{\Gamma }}},

and finally coupling matrix is defined as

(29){{\mathit{Q}}_{ij}}^e={\displaystyle \int {\mathit{N}}_{\mathrm{s}}\mathit{n}{\mathit{N}}_{\mathrm{f}}^{\mathrm{T}}\mathbf{d}\mathrm{\Gamma }}.

In the integration on the wet interface, the normal vector \mathit{n}is defined to be positive going from solid into the fluid.

Nonlinear fluid-structure coupling scheme

Coupling between domains

The coupling between the fluid and structure domains is applied through the forcing terms. Since have been already assumed the fluid to be inviscid, the coupling occurs only in the normal direction. It is clear that the surface traction acting on the fluid due to the interaction with the structure is equal and opposite of the pressure loading on the structure by fluid. The work done by hydrodynamic pressure on the interaction surface of the structure must be equal to the work of the equivalent nodal forces on the interface boundary of the fluid element. The coupling matrix (Eq. (29)) relates the pressure of the reservoir and the forces on the dam-reservoir interface as

(30)\mathit{Q}\mathit{P}=\mathit{F},

where \mathit{F} is the force vector acting on the structure due to the pressure loading.

Governing equations

Equations (10) and (27) describe the finite element discretized equations for the dam-reservoir interaction problem and can be written in an assembled form as

\left[\begin{array}{cc}\hfill \mathit{M}\hfill & \hfill \mathbf{0}\hfill \\ \hfill -{\rho }_{\mathrm{f}}\mathit{Q}\hfill & \hfill \mathit{G}\left(\mathit{p}\right)\hfill \end{array}\right]\left\{\begin{array}{c}\hfill {\ddot{\mathit{u}}}^{n+\mathrm{1}}\hfill \\ \hfill {\ddot{\mathit{P}}}^{n+\mathrm{1}}\hfill \end{array}\right\}+\left[\begin{array}{cc}\hfill \mathit{C}\hfill & \hfill \mathbf{0}\hfill \\ \hfill \mathbf{0}\hfill & \hfill \mathit{D}\left(\mathit{P}\right)\hfill \end{array}\right]\left\{\begin{array}{c}\hfill {\dot{\mathit{u}}}^{n+\mathrm{1}}\hfill \\ \hfill {\dot{\mathit{P}}}^{n+\mathrm{1}}\hfill \end{array}\right\}

(31)+\left[\begin{array}{cc}\hfill {\mathit{K}}_{\mathrm{t}}\hfill & \hfill \mathit{Q}\hfill \\ \hfill \mathbf{0}\hfill & \hfill \mathit{H}\hfill \end{array}\right]\left\{\begin{array}{c}\hfill {\mathit{u}}^{n+\mathrm{1}}\hfill \\ \hfill {\mathit{P}}^{n+\mathrm{1}}\hfill \end{array}\right\}=\left\{\begin{array}{c}\hfill \mathit{Q}\left({\mathit{P}}_{\mathrm{i}}^{n+\mathrm{1}}+{\mathit{P}}_{\mathrm{r}}^{n+\mathrm{1}}\right)\hfill \\ \hfill {\mathit{F}}_{\mathrm{f}}\hfill \end{array}\right\},

where \mathit{u},\dot{\mathit{u}} and\ddot{\mathit{u}} are structure nodal displacements vector and its first and second time derivatives, respectively, and \mathit{p},\dot{\mathit{p}} and \ddot{\mathit{p}} are hydrodynamic pressure vector of reservoir domain and its first and second time derivatives, respectively. p_iis the impinging shock pressure wave that computed using Eq. (1) and p_r is the maximum reflected overpressure calculated by Eqs. (5) and (6). K_t, the tangent stiffness matrix of structural domain, is computed with considering material nonlinearity behavior of concrete and strain rate effect.

Validation of simulation

Validation for dynamic analysis of simple coupled system

The objective of this example is to assess the accuracy of the applied fluid-structure interaction solution method. Hence, a benchmark problem has been solved and compared with the existing literature. Dynamic interaction response of a simple 3D fluid-structure system that is illustrated in Fig. 3(a) is investigated and subjected to a ramp excitation (see Fig. 3(b)). The data assumed for the problem are: depth of the fluid domain is, H=\mathrm{180} m, the width of the fluid domain is, B=\mathrm{20} m, the reservoir behind the structure extends a distance of L=\mathrm{3}H=\mathrm{540} m, and the thickness of wall structure is,t=\mathrm{15} m. The speed of sound and mass density of water are taken as, c_w=\mathrm{1440} m/s and {\rho }_f=\mathrm{1000}{\mathrm{kg/m}}^3, respectively. The material of structure is assumed to be linearly elastic. Its elasticity modulus, mass density and Poisson’s ratio are taken as \mathrm{3.43}\times {\mathrm{10}}^{\mathrm{10}}{\mathrm{kg/m}}^2, \mathrm{2440}{\mathrm{kg/m}}^3and 0\mathrm{.2}, respectively. Figure 3(c) shows a typical finite element discretization of the fluid- structure system. The fluid domain is discretized by 576 eight-node brick elements with linear pressure interpolation functions and the structure domain is modeled using 24 twenty-node brick elements with quadratic displacement interpolation function. This simple example was previously studied by other researchers analytically [⁵²,⁵³] and numerically [⁵⁴]. The computed results of the analysis of dynamic interaction of flexible wall with adjacent reservoir are illustrated in Figs. 4(a) and 4(b). The results are compared with those presented by Tsai et al. [⁵²], which validates the proposed algorithm at the case of linear fluid behavior. Fundamental frequencies of wall structure for empty and full reservoir are calculated to be 11.42 and 6.68 rad/s, respectively. Same wall structure-fluid system assuming the wall to be rigid is also analyzed. The time history of dynamic pressure is shown in Fig. 5 and compared with result obtained analytically by Tsia et. al. [⁵³], and good agreement is shown between two results.

Validation for shock wave loading and cavitation effects

This example consists of a floating structure that is subjected to a step exponential plan wave and is one-dimensional Bleich- Sandler problem [⁵⁵]. This test case is often referred to in literature, because it has been solved analytically and validated by many other researchers [⁵⁶–⁵⁹]. The structure is modeled using five 20-node brick elements while water under the structure is modeled by 100, 8-node fluid volume elements as is shown in Fig. 6. The depth of the fluid is 3.81 m. The peak pressure, P_{\mathrm{0}}, is 0.710 MPa and the decay time,\theta is 0.9958 ms. The density and sound speed of the fluid in this case are 1000 kg/m³ and 1450 m/s, respectively. The atmospheric pressure is 1 Pa. These material properties and parameters are identical to those used in the papers published by Sprague and Geers [⁶⁰]. Figure 7 illustrates the velocity history of the floating structure center for both cases of with and without cavitation, which is well fitted with Bleich-Sandler analytic solution [⁵⁵].

Examples and numerical simulation

In this study, two individual concrete dams are selected as representative concrete dams in order to study the effect of air blast on the dynamic performance of concrete dams. The first is Konya Dam which that is a large concrete gravity dam located in Maharashtra, India. Built over 1954‒1963, it has a height of 103 m, width of 808 m, volume of 1555 m³ and water storage capacity of 2797 km³. Figure 8 shows a representative coarse and fine finite element model of the concrete gravity dam along with relevant dimensions. The second benchmark dam is the Morrow point arch dam on the Gunnison River located in Colorado with a height of 143 m. The finite element model of Morrow point arch dam is presented in Fig. 9. Material of the dams is assumed to be nonlinear viscoplastic, isotropic and homogeneous. The water is assumed to be nonlinear acoustic due to cavitation. Fluid and structure are modeled using eight-node and twenty-node hexahedron elements, respectively. In fluid domain, mesh is truncated at non-reflecting boundary. Physical parameters chosen for the system are as follows. For the dam concrete:E_c=\mathrm{34.3}\times {\mathrm{10}}^{\mathrm{9}}\mathrm{Pa}, {\nu }_c=\mathrm{0.2}. For the reservoir water: c=\mathrm{1439}\mathrm{m}/\mathrm{s} and {\rho }_w=\mathrm{1019}\mathrm{kg}/{\mathrm{m}}^3.

Choosing the mesh size of dam finite element models

In Refs. [¹⁷,¹⁸,⁶¹] discussed the effect of mesh size on the simulation and prediction the response of structures subjected to explosion load. They believed that if the grid size is taken fine, the spread way of explosion load can be simulated accurately. This is mainly because of accurately capturing the transition of the shock wave energy on the structure in the explosion. In the present research, to overcome to this need, the blast loading is assumed as conventional weapon blast loading, therefore, blast loading is applied on the dam structure through an equivalent pressure time history generated using the manual and modified Friedlander’s equation (Eq. (1)), namely, the effective blast induced pressure load is computed with summation of impinging pressure p_i(t) (Eq. (1)) and reflected overpressure p_r(t) (Eq. (5) or Eq. (6)) that applied on the structure on the assumed affected zone as a distributed pressure loading, therefore, only nodal points in this zone are loaded due to blast loading. In the coarse and refined finite element (FE) meshes the intensity of distributed pressure load and area of affected zone are identical only the number of loaded nodal points is varied but with different magnitude of equivalent nodal loads in two meshes. With uses this blast load applying method, the responses of coarse and refined FE model of dams do not alter significantly and can be used relatively coarse FE mesh in the simulation as shown in the following sections.

(32)P_{\mathrm{blast}}=P_{\mathrm{i}}(t)+P_{\mathrm{r}}(t)\Rightarrow F_{\mathrm{blast}}\mathrm{(t)}=P_{\mathrm{blast}}\times A_{\mathrm{area\; of\; affected\; zone}}=\mathrm{Equivalent\; nodal\; loads}.

Concrete gravity dam: Konya dam

Four fully coupled numerical models with heights of 50, 80, 110, and 150 m, respectively, are built in order to study the effect of the dam height on the dynamic behavior of concrete gravity dam under air blast loading. As depicted in Fig. 10, two different affected zones are assumed in downstream face of dam, one close to crest in the upper part of dam and later close to middle part of downstream face of dam. For each case, explosive charge weight selected depends on the tolerated limit of each dam. It should be noted that the tolerated limit of each dam model is computed through several incremental analysis.

In the FE model with coarse mesh, is used 1728, 8-node fluid elements and 768, 20-node brick elements to model fluid-structure system and refined mesh is included 7200, 8-node fluid element and 3600, 20-node brick elements to discretized coupled dam-reservoir system. The location of far end boundary of the reservoir domain is chosen through the sensitivity analysis which indicated that assuming reservoir length equal to 210 m from upstream face of dam is sufficient and utilizes larger length has no effect on the dam's response.

To obtain accurate results and capture shock wave propagation time step, \Delta t, is assumed \mathrm{2}\times {\mathrm{10}}^{-\mathrm{5}} in the analysis. The response of the dam to the applied air blast loading that is computed using coarse and refined FE meshes is illustrated in Fig. 11. Obtained results indicate that the response of dam is independent of the numerical mesh chosen in the simulation. Fig. 12 depicted time history graphs of crest displacement of the concrete gravity dam for the first 3 s of air blast loading close to the upper part of dam. As shown, the displacement of dam crest is observed to fluctuate with time. Three different conditions, namely, full reservoir without considering cavitation, full reservoir with cavitation and empty reservoir are considered in Fig. 12. Cavitation has affected the response of dam and increases the displacement, but this effect decreases with decrease of dam height, that is, cavitation increases displacement by14.4% for dam height which is equal to 150 m and 5.3% for dam height which is equal to 50 m. The reason behind this is that the cavitating region which is adjacent to the upstream face of dam can create and expand in deep reservoirs and collapsing of this region resulting in a reloading of dam structure. In principle, the displacement of dam crest increases with the increase of the dam height. As a result, the maximum displacement of the dam crest for the dam height of 150 m is 9.36 cm toward the upstream direction. With decrease of dam height, the fluctuation nature of dam response decrease and collapsing instant displacement also decrease which is due to the increase in the rigidity of dam structure and reduction of flexibility of dam with decrease of dam height. The analysis of dam response which was carried out with full and empty reservoir, as shown in Fig. 12 illustrates that the interaction of reservoir and dam structure decreases the response of dam subjected to air blast loading, as the dam crest displacement decreases from 17.5%, 10.9%, 7.9%, and 5% for dam height of 150, 110, 80, and 50 m, respectively.

Obtained results reveal that the maximum weight of charge that is located close to the upper part of downstream face of gravity dam that collapse gravity dam decreases with decrease of dam height, thus for dam heights 150 m (R=\mathrm{10}\mathrm{m},h=\mathrm{135}\mathrm{m}), 110 m (R=\mathrm{10}\mathrm{m},h=\mathrm{96}\mathrm{m}), 80 m (R=\mathrm{10}\mathrm{m},h=\mathrm{66}\mathrm{m}), 50 m (R=\mathrm{10}\mathrm{m},h=\mathrm{36}\mathrm{m}), the maximum weight of charges obtained are 660, 400, 320, and 270 kg, respectively. Figure 13 shows the number of cavitating nodes in the reservoir during analysis. It can be observed from Fig. 13 that the cavitating region continues to expand beyond the point of hydrostatic equilibrium until a point of dynamic equilibrium is attained. The cavitating region then reverses its motion and continues to contract until dynamic equilibrium is again attained, where it quickly rebounds and again begins to expand. This oscillating cavitating region expansion and contraction continues until the energy of the reaction is fully dissipated or the cavitating region finally reaches the surface, venting the gaseous by products of the explosion. As the cavitating region rebounds, it greatly accelerates the surrounding water, generating a substantial pressure pulse, known as the bubble pulse. This bubble pulse can impart significant loads on structures in the vicinity. Also Fig. 13 illustrates that the number of cavitating nodes in the reservoir in the first duration shock wave induced by air blast reduces with decrease of dam height. The time history graphs of the hydrodynamic pressure at the point C which is located at the upper part of reservoir in the vicinity of upstream face of dam is shown in Fig. 14. As illustrated in Fig. 14(a), that is computed for air blast loading close to the upper part of dam, hydrodynamic pressure attains its maximum value in the early stage of analysis and then decrease rapidly and tend to zero, as hydrodynamic pressure maximum value obtained for dam with height 150 m equal to 8057.74 kPa but its value reduces to 4293.84 kPa for height decrease to 110 m, a 26.7% decrease in the height of dam induced 46.7% reduction in hydrodynamic pressure value, but with further decrease in the height of dam, it has lesser effect on hydrodynamic pressure value. The reason is that the deflection of the dam body toward the lateral direction is larger with the increase of the dam height and this result in a relatively strong coupling between dam structure and reservoir that induced larger hydrodynamic pressure with larger height. Comparisons between Figs. 14(a) and 14(b) reveal that change of standoff point of charge has significant effect on the hydrodynamic pressure distribution in the reservoir. For the two different standoff positions for dam height 150 m; one is (R=\mathrm{10}\mathrm{m},h=\mathrm{135}\mathrm{m}) and the other is (R=\mathrm{10}\mathrm{m},h=\mathrm{72}\mathrm{m}), the obtained maximum value of hydrodynamic pressures are 8057.74 kPa and 1322.68 kPa, respectively, indicating that the maximum value of hydrodynamic pressure is almost six times reduced with change standoff location from the upper to middle part of the dam. The failure elements of gravity dams under 660, 400, 320, and 270 kg TNT charge for dam height 150, 110, 80, and 50 m, respectively, are shown in Fig. 15, with standoff of charge close to upper part of gravity dam, at all cases the standoff distance is fix to 10 m. It can be seen from Fig. 15 that localized damage of concrete is observed on the downstream face of the dam facing the explosion center and the detonation pressure acts on the downstream surface of the dam. The failure modes of the dam are similar with each other in different cases and changing the dam height not affected the failure mode of gravity dam. As illustrated in Fig. 16, the time history diagrams of displacement of gravity dam crest at Point A, is depicted at different detonation height ranging from 24 to 72 m (h=24, 42, 54, and 72 m, for dam heights 50, 80, 110, and 150 m, respectively), The standoff distance is 10 m (R=\mathrm{10}\mathrm{m}). Three different situations namely full reservoir without considering cavitation, full reservoir with cavitation forming and empty reservoir are considered in Fig. 16. Cavitation changes the response of dam and intensifies the maximum displacement of dam crest. The influence of cavitation increases with increase of dam height. As shown in Fig. 16, cavitation has moderately little effect in the increase of dam crest displacement, that is, cavitation increases displacement of dam crest by 12% for dam height of 150 m and 9.5% for dam height of 50 m. It can be seen in Fig. 16 that the displacement of dam crest increases with increase of the dam height, and therefore, the maximum displacement of the dam crest for the dam height of 150 m is 4.9 cm toward the upstream direction. As shown in Fig. 16, interaction of reservoir and dam structure decreases the response of dam subjected to air blast loading, as the dam crest displacement decreases by 14.2%, 9.7%, 7.7%, and 2.3% for dam heights of 150, 110, 80, and 50 m, respectively. As can be seen from Fig. 16, the maximum weight of charge that is located close to the middle part of downstream face of gravity dam that resulted in the collapse of the dam decreases with decrease of dam height, thus for dam heights of 150, 110, 80, and 50 m, the maximum weight of charges obtained are 925, 535, 475, and 420 kg, respectively. The failure elements of gravity dams under 925, 535, 475, and 420 kg TNT charge for dam height H=150, 110, 80, and 50 m, respectively, are shown in Fig. 17, with standoff of charge close to middle part of gravity dam, at all cases the standoff distance is fix to 10 m. It can be seen form Fig. 17 that the crack propagates in the dam body and penetrate in the whole section of the dam with decrease the dam height, as in the case of dam height is 50 m the crack around the elevation at which the slope of downstream face changes abruptly extends completely to the upstream face. However, the failure modes of the dam are similar with each other in different cases and changing the dam height not affected significantly the failure mode of gravity dam.

Concrete arch dam: Morrow point dam

The interest in this example is directed towards application of the proposed method for the evaluation of response of arch dams under blast loading. Morrow point dam depicted in Fig. 9, is a double curved concrete arch dam chosen for bench mark problem in the present study. Although it has a height of 142 m, but it is assumed four different heights of 100, 150, 200, and 250 m in order to study the effect of the dam height on the blast loading behavior of concrete arch dams. As shown in Fig. 18, two different locations for charge are assumed, one is in the upper part of dam close to crest and the other is in the middle of dam. Figure 9 shows coarse and refined 3D FE model of coupled system along with its relevant dimensions. It is employed 13360, 8-node fluid element in the model of reservoir and 800, 20-node brick element for discretization of dam structure in the coarse FE model and 49600, 8-node fluid element in the model of reservoir and 3200, 20-node brick element for discretization of arch dam structure in the refined mesh. Two element layers are used through cross direction of the dam thickness to determine the damage of upstream and downstream faces. A time step, \Delta t=\mathrm{2}\times {\mathrm{10}}^{-\mathrm{5}}, is used in the dynamic analysis. Figure 19 illustrates the response of morrow point arch dam in the first 2.5 s of applied air blast loading at the upper part of the arch dam that is evaluated using coarse and refined FE meshes. Time history graphs of the displacement of Point A, which is located at the crest of dam, were depicted in Fig. 20 for different standoff points. For height of 250 m, standoff point coordinate is(R=\mathrm{10}\mathrm{m},h=\mathrm{235}\mathrm{m}), for height of 200 m, standoff coordinate is (R=10 m, h=185 m), for height of 150 m standoff point coordinate is (R=10 m, h=135 m) and for height of 100 m, standoff point coordinate is (R=10 m, h=86 m). It can be observed from Fig. 20 that the maximum displacement of the crest of arch dam for height of 250 m is 12.4 cm toward upstream face of dam with considering cavitation, but crest displacement attained 10.5 cm without considering cavitation which shows an 18.1% increase of displacement due to bulk cavitation forming in the reservoir. Cavitation affects the response of high arch dams, as for heights of 200, 150, and 100 m, crest displacement increases were 15.5%, 11.5% and 15%, respectively. But in all, cavitation effects decrease with decrease of dam height. The reason is that the cavitating region adjacent to upstream face of high arch dam expanded violently and as a result, the collapsing of this region induced larger reloading on dam structure. Other aspect that is studied in Fig. 20 is the reservoir effects on the coupled system performance that arch dam reservoir systems analyses with full reservoir and empty reservoir. Obtained results revealed that reservoir reduces the maximum crest displacement of dam, and as for arch dam heights of 250, 200, 150, and 100 m, the crest displacement decreases are 7.9%, 10.9%, 12.5%, and 18.3%, respectively. Reduction effect of reservoir on arch dam dynamic response increases with decreasing of arch dam height which is contrary to concrete gravity dam behavior. As shown in Fig. 20, the maximum charge weight located close to crest of arch dam that resulted in the collapse of the dam decrease with reduction of dam height, and as for dam heights of 250, 200, 150, and 100 m, the maximum charge weight is computed as 150, 140, 140, and 135 kg, respectively. In contrast to gravity dam, the maximum charge weight for arch dam decreases very smoothly with decreasing of height of dam.

The damaged element of arch dam in downstream and upstream face of dam body for 150, 140, 140, and 135 kg TNT charge for dam heights H=250, 200, 150, and 100 m, respectively and with standoff distance 10 m (R=10 m) and standoff height (h=235, 185, 135, and 86 m, respectively) are shown in Fig. 21. As detonation occur close to upper part of dam crest, damaged elements propagate at the center of arch dam body close to dam crest. With decrease of arch dam height, the pattern of damage is changed, and failure elements gradually propagate along the dam crest and damaged zone is extended.

Figure 22 shows time history diagrams of the crest displacement of the arch dam corresponding to the point A for different dam heights and for charge located at the middle of arch dam. Consequently, for dam height of 250 m, standoff point coordinate is (R=10 m, h=130 m), for dam height of 200 m, standoff point coordinate is (R=10 m, h=108 m), for dam height of 150 m, standoff point coordinate is (R=10 m, h=80 m) and for dam height of 100 m, standoff point coordinate is (R=10 m, h=48 m). As illustrated in Fig. 22, the displacement of dam crest pattern obtains different from in Fig. 20 and cavitation totally changes the crest displacement manner and intensifies the crest displacement. This is for the complicated nature of arch dam —reservoir system and interaction between different domains. It can be seen from Fig. 22 that the maximum displacement of the crest of arch dam for dam height 250 m is 4.9 cm, for dam height 200 m is 4.8 cm and for dam height 150 m is 4 cm toward downstream face of dam with considering cavitation and for dam height 100 m the maximum crest displacement obtain 3.3 cm toward upstream face of dam with considering cavitation. As presented in Fig. 22, the maximum charge weight located close to the middle of arch dam that resulted in collapse of dam decreases with reduction of dam height, and as for dam heights of 250, 200, 150, and 100 m, the maximum charge weights computed are 200, 200, 185, and 175 kg, respectively.

The damaged element of arch dam in downstream and upstream face of dam body for 200, 200, 185, and 175 kg TNT charge for dam heights H=250, 200, 150, and 100 m, respectively and with standoff distance 10 m (R=10 m) and standoff height (h=130, 108, 80, and 48 m, respectively) are shown in Fig. 23. As shown in Fig. 23, the blast pressure generated from air blast loading causes a highly localized damage close to the middle part of downstream face of dam body facing the explosion center and causes a confined damage zone in the downstream face of dam close to abutments. The crack damage in the upstream face, near the bank of dam, extends due to tensile stress arise from blast loading and the lower tensile strength of concrete material than its compressive strength. Figure 24 shows hydrodynamic pressure time history graphs at Point C that is located close to the crest of arch dam for two different charges position. As shown in Fig. 24(a), that when charge is located close to the upper part of arch dam, hydrodynamic pressure rises rapidly and attains its maximum value in the early time of blast and then decreases and tend to zero, but as shown in Fig. 20(b), hydrodynamic pressure has a fluctuated nature when charge is located close to the middle of arch dam. It can be seen from Fig. 24(a) that the maximum value of hydrodynamic pressure at Point C for heights of 250, 200, 150, and 100 m are 478.32, 462.73, 428.47, and 414.90 t/m², respectively. As shown in Fig. 24(b), the maximum value of hydrodynamic pressure at Point C for heights of 250, 200, 150, and 100 m are 148.2, 140.12, 111.0, and 96.0 t/m², respectively. In comparison with hydrodynamic pressure variation that is depicted in Fig. 14, hydrodynamic pressure reduces smoothly with decrease of height of arch dam. Maximum principal stress contour of arch dam is shown in Fig. 25 for standoff point of charge located close to the upper part of arch dam and for different heights of dam. It can be seen from Fig. 25 that region with tensile stresses in dam body is created in the middle part of dam close to the upper part confronting charge standoff point and in the bottom of dam body but region with compressive stress is formed close to the banks of arch dam. Figure 26 shows maximum principal stress contour of arch dam for standoff point of charge located close to the middle of arch dam and for different heights of dam. As shown in Fig. 26, region of dam body with tensile stress is bounded to the lower part of dam body close to the bottom of dam and compressive stress region is created on the banks.

Conclusions

In the present study, the three-dimensional finite element developed Fortran code has been used to investigate the effect of air blast loading on concrete dams and analysis procedure is done for typical concrete gravity dam and concrete arch dam with different heights and different standoff point of charge with empty and full reservoirs. In addition, bulk cavitation effects on response of concrete dams were studied. The following is a brief outline of the key results obtained from this study:

1) For charge located adjacent to the upper part of the gravity dam, displacement of crest attains its maximum value within very short lag time after air blast shock loading and with changing charge location which is adjacent to the middle part of the dam, the response scheme of dam does not vary. But for arch dam with changing charge location close to the middle part of the dam, the response scheme of the arch dam changed completely.

2) Charge weight that led to collapse of both types of concrete dams increases with increasing height of concrete dams, but collapsing charge weight of gravity dam is approximately four times of arch dam for examples considered in this research. This indicated that concrete arch dams are more vulnerable to air blast loading than concrete gravity dams.

3) The inclusion of reservoir in the numerical model decreases the response of concrete dams under air blast. Increase in the height of concrete gravity dam resulted in the increase in reservoir reduction effects but for arch dam, increase in the height of dam decreases reduction effect of reservoir. Maximum reduction effect of reservoir is approximately 18% for both types of concrete dams.

4) Acoustic cavitation increases the response of con crete dams under blast loading. As cavitation intensifies crest displacement of concrete gravity dam and concrete arch dam at the maximum of 14.4% and 15.5%, respectively.

5) Compressive area generally appeared close to the sides of concrete arch dam, but tensile area is created in the central part and bottom of arch dam body for standoff point of charge located close to the upper part of arch dam and with changing of standoff point location close to the middle part of arch dam, tensile area was confined to the lower part of dam body.

References

[1] Ngo T, Mendis P, Gupta A, Ramsay J. Blast loading and blast effects on structures—An overview. Electronic Journal of Structural Engineering, 2007, 7: 76–91

[2] Remennikov A M. A review of methods for predicting bomb blast effects on buildings. Journal of Battlefield Technology, 2003, 6(3): 5–10

[3] Lu Y, Wang Z. Characterization of structural effects from above-ground explosion using coupled numerical simulation. Computers & Structures, 2006, 84(28): 1729–1742

[4] Tian L, Li Z X. Dynamic response analysis of a building structure subjected to ground shock from a tunnel explosion. International Journal of Impact Engineering, 2008, 35(10): 1164–1178

[5] Jayasooriya R, Thambiratnam D P, Perera N J, Kosse V. Blast and residual capacity analysis of reinforced concrete framed buildings. Engineering Structures, 2011, 33(12): 3483–3492

[6] Parisi F, Augenti N. Influence of seismic design criteria on blast resistance of RC framed buildings: A case study. Engineering Structures, 2012, 44: 78–93

[7] Tang E K C, Hao H. Numerical simulation of a cable-stayed bridge response to blast loads. Part I: Model development and response calculations. Engineering Structures, 2010, 32(10): 3180–3192

[8] Hao H, Tang E K C. Numerical simulation of a cable-stayed bridge response to blast loads. Part II: Damage prediction and FRP strengthening. Engineering Structures, 2010, 32(10): 3193–3205

[9] Son J, Lee H J. Performance of cable-stayed bridge pylons subjected to blast loading. Engineering Structures, 2011, 33(4): 1133–1148

[10] Anwarul Islam A K M, Yazdani N. Performance of AASHTO girder bridges under blast loading. Engineering Structures, 2008, 30(7): 1922–1937

[11] Lu Y, Wang Z, Chong K. A comparative study of buried structure in soil subjected to blast load using 2D and 3D numerical simulations. Soil Dynamics and Earthquake Engineering, 2005, 25(4): 275–288

[12] Wang Z, Lu Y, Hao H, Chong K. A full coupled numerical analysis approach for buried structures subjected to subsurface blast. Computers & Structures, 2005, 83(4–5): 339–356

[13] Ma G, Zhou H, Chong K. In-structure shock assessment of underground structures with consideration of rigid body motion. Journal of Engineering Mechanics, 2011, 137(12): 797–806

[14] Li J C, Li H B, Ma G W, Zhou Y X. Assessment of underground tunnel stability to adjacent tunnel explosion. Tunnelling and Underground Space Technology, 2013, 35: 227–234

[15] Zhang A, Zeng L, Cheng X, Wang S, Chen Y. The evaluation method of total damage to ship in underwater explosion. Applied Ocean Research, 2011, 33(4): 240–251

[16] Jin Q, Ding G. A finite element analysis of ship sections subjected to underwater explosion. International Journal of Impact Engineering, 2011, 38(7): 558–566

[17] Zhang A, Zhou W, Wang S, Feng L. Dynamic response of the non-contact underwater explosions on naval equipment. Marine Structures, 2011, 24(4): 396–411

[18] Wang G, Zhang S. Damage prediction of concrete gravity dams subjected to underwater explosion shock loading. Engineering Failure Analysis, 2014, 39: 72–91

[19] Zhang S, Wang G, Wang C, Pang B, Du C. Numerical simulation of failure modes of concrete gravity dams subjected to underwater explosion. Engineering Failure Analysis, 2014, 36: 49–64

[20] De A. Numerical simulation of surface explosions over dry, cohesionless soil. Computers and Geotechnics, 2012, 43: 72–79

[21] Falconer J. The Dam Busters Story. London: Sutton Publishing Ltd., 2007

[22] Xue X, Yang X, Zhang W. Numerical modeling of arch dam under blast loading. Journal of Vibration and Control, 2012, 20(2): 256–265

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

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

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

[26] Hall J F. Study of the earthquake response of pine flat dam. Earthquake Engineering & Structural Dynamics, 1986, 14(2): 281–295

[27] Kalateh F, Attarnejad R. A new cavitation simulation method: Dam-reservoir systems. International Journal for Computational Methods in Engineering Science and Mechanics, 2012, 13(3): 161–183

[28] Guzas E L, Earls C J. Air blast load generation for simulating structural response. Steel and Composite Structures, 2010, 10(5): 429–455

[29] Baker W E. Explosions in Air. Austin: University of Texas Press, 1973

[30] Baker W E, Cox P, Westine P S, Kulesz J J, Strehlow R A. Explosion Hazards and Evaluation. New York: Elsevier, 1983

[31] Kinney G F, Graham K J. Explosive Shock in Air. 2nd ed. New York: Springer, 1985

[32] Kingery C N, Bulmash G. Airblast parameters from TNT Spherical Air Burst and Hemispherical Surface Burst. Army BRL Report ARBRL-TR-02555, 1984

[33] Smith P D, Hetherington J G. Blast and Ballistic Loading of Structures. London: Butterworth-Heinemenn, 1994

[34] Borenstein E, Benaroya H. Sensitivity analysis of blast loading parameters and their trends as uncertainty increases. Journal of Sound and Vibration, 2009, 321(3–5): 762–785

[35] Brode H L. Quick Estimates of Peak Overpressure from Two Simultaneous Blast Waves. Topical Report Oct-Dec 77, 1977

[36] ANSYS Inc. Theory Reference, Release 10.0 Documentation for ANSYS software ANSYS Inc., 2006

[37] Eibl J, Schmidt-Hurtienne B. Strain-rate sensitive constitutive law for concrete. Journal of Engineering Mechanics, 1999, 125(12): 1411–1420

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

[39] Bischoff P H, Perry S H. Compressive behavior of concrete at high strain rate. Materials and Structures, 1991, 24(6): 425–450

[40] Malvar L J, Ross C A. Review of strain rate effects for concrete in tension. ACI Materials Journal, 1998, 95(M73): 735–739

[41] Wu Y, Wang D, Wu C T. Three dimensional fragmentation simulation of concrete structures with a nodally regularized meshfree method. Theoretical and Applied Fracture Mechanics, 2014, 72: 89–99

[42] Wu Y, Wang D, Wu C T, Zhang H. A direct displacement smoothing meshfree particle formulation for impact failure modeling. International Journal of Impact Engineering, 2016, 87: 169–185

[43] Zhou X Q, Hao H, Deeks A J. Modelling dynamic damage of concrete slab under blast loading. In: Proceedings of the 6th International Conference on Shock and Impact Loads on Structures. Perth: CI-Premier, 2005, 703–710

[44] Zhou X Q, Kuznetsov V A, Hao H, Waschl J. Numerical prediction of concrete slab response to blast loading. International Journal of Impact Engineering, 2008, 35(10): 1186–1200

[45] Herrmann W. Constitutive equation for the dynamic compaction of ductile porous materials. Journal of Applied Physics, 1969, 40(6): 2490–2499

[46] Xie W F, Liu T G, Khoo B C. Application of a one-fluid model for large scale homogenous unsteady cavitation: The modified Schmidt model. Computers & Fluids, 2006, 35(10): 1177–1192

[47] Sandberg G. A new finite element formulation of shock-induced hull cavitation. Computer Methods in Applied Mechanics and Engineering, 1995, 120(1–2): 33–44

[48] Brennen C E. Cavitation and Bubble Dynamics. Oxford: Oxford University Press, 1995

[49] Kalateh F, Attarnejad R. Finite element simulation of acoustic cavitation in the reservoir and effects on dynamic response of concrete dams. Finite Elements in Analysis and Design, 2011, 47(5): 543–558

[50] Xie W F, Liu G, Khoo B C. The simulation of cavitation flows induced by underwater shock and free surface interaction. Applied Numerical Mathematics, 2007, 57(5–7): 734–745

[51] Wallis G B. One-dimensional Two-phase Flow. New York: McGraw-Hill, 1969

[52] Tsai C S, Lee G C, Ketter R L. A semi-analytical method for time domain analyses for dam-reservoir interactions. International Journal for Numerical Methods in Engineering, 1990, 29(5): 913–933

[53] Tsai C S, Lee G C, Yeh C S. Time-domain analysis of three-dimensional dam-reservoir interaction by BEM and semi-analytical method. Engineering Analysis with Boundary Elements, 1992, 10(2): 107–118

[54] Gogoi I, Maity D. A non-reflecting boundary condition for the finite element modeling of infinite reservoir with layered sediment. Advances in Water Resources, 2006, 29(10): 1515–1527

[55] Bleich H H, Sandler I S. Interaction between structures and bilinear fluids. International Journal of Solids and Structures, 1970, 6(5): 617–639

[56] Zienkiewicz O C, Paul D K, Hinton E. Cavitation in fluid-structure response (with particular reference to dams under earthquake loading). Earthquake Engineering & Structural Dynamics, 1983, 11(4): 463–481

[57] Newton R E. Finite element study of shock induced cavitation. In: Oñate E, Periaux J, Samuelsson A, eds. The Finite Element Method in the 1990’s. Berlin: Springer, 1991, 389–397

[58] Hamdi M A, Ousset Y, Verchery G. A displacement method for the analysis of vibrations of coupled fluid-structure systems. International Journal for Numerical Methods in Engineering, 1978, 13(1): 139–150

[59] Felippa C A, Deruntz J A. Finite element analysis of shock-induced hull cavitation. Computer Methods in Applied Mechanics and Engineering, 1984, 44(3): 297–337

[60] Sprague M A, Geers T L. Spectral elements and field separation for an acoustic fluid subject to cavitation. Journal of Computational Physics, 2003, 184(1): 149–162

[61] Luccioni B, Ambrosini D, Danesi R. Blast load assessment using hydrocodes. Engineering Structures, 2006, 28(12): 1736–1744

RIGHTS & PERMISSIONS

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