Space-Time Finite Element Method for Seismic Analysis of Concrete Dam

Abstract

Finite element method (FEM) is the most extended approach for analyzing the design of the dams against earthquake motion. In such simulations, time integration schemes are employed to obtain the response of the dam at time tn+1 from the known response at time tn. To this end, it is desirable that such schemes are high-order accurate in time and remain unconditionally stable large time-step size can be employed to decrease the computation cost. Moreover, such schemes should attenuate the high-frequency components from the response of structure being studied. Keeping this in view, this chapter presents the theory of time-discontinuous space-time finite element method (ST/FEM) and its application to obtain the response of dam-reservoir system to seismic loading.

1. Introduction

During an event of earthquake stability of dams is of paramount importance as their failure can cause immense property and environmental damages. When dam-reservoir-foundation system is subjected to the dynamic loading it causes a coupled phenomenon; ground motion and deformations in the dam generate hydrodynamic pressure in the reservoir, which, in turn, can intensify the dynamic response of the dam. Moreover, spatial-temporal variation of stresses in the dam-body depends on the dynamic interactions between the dam, reservoir, and foundation. Therefore, it becomes necessary to use numerical techniques for the safety assessment of a given dam-design against a particular ground motion.

Dynamic finite element method is the most extended approach for computing the seismic response of the dam-reservoir system to the earthquake loading [1]. In this approach finite elements are used for discretization of space domain, and basis functions are locally supported on the spatial domain of these elements and remain independent of time. Furthermore, nodal values of primary unknowns depend only on time. Accordingly, this arrangement yields a system of ordinary differential equations (ODEs) in time which is then solved by employing time-marching schemes based on the finite difference method (FDM), such as Newmark-$\beta$ method, HHT-$\alpha$ method, Houbolt method, and Wilson-$\theta$ method.

In dynamic finite element method (FEM), it is desirable to adopt large time-steps to decrease the computation time while solving a transient problem. Therefore, it is imperative that the time-marching scheme remains unconditionally stable and higher order accurate [1]. In addition, it should filter out the high frequency components from the response of structure. To achieve these goals, Hughes and Hulbert presented space-time finite element method (ST/FEM) for solving the elastodynamics problem [2]. In this method, displacements $\left(\mathbf{u}\right)$ and velocities $\mathbf{v}$ are continuous in space-domain, however, discontinuous in time-domain. Li and Wiberg incorporated the time-discontinuity jump of displacements and velocities in the total energy norm to formulate an adaptive time-stepping ST/FEM [3]. ST/FEM, so far, has been successfully employed for solving linear and nonlinear structural dynamics problems [4, 5], moving-mass problems [6], and dynamical analysis of porous media [7], among other problems.

However, for elastodynamics problem, ST/FEM, yields a larger system of linear equations due to due to the time-discontinuous interpolation of displacement and velocity fields. Several efforts have been made in the past to overcome this issue; both explicit [8] and implicit [9] predictor-multi-corrector iteration schemes have been proposed to solve linear and nonlinear dynamics problems. Recently, to reduce the number of unknowns in ST/FEM, a different approach is taken in which only velocity is included in primary unknowns while displacement and stresses are computed from the velocity in a post-processing step [4, 10]. To this end, the objective of the present chapter is to introduce this method (henceforth, ST/FEM) in a pedagogical manner. The rest of the chapter is organized as follows. Sections 2 and 3 deal with the fundamentals of time-discontinuous Galerkin method. Section 4 describes the dam-reservoir-soil interaction problem, and Section 5 discusses the application of ST/FEM for this problem. Lastly, Section 6 demonstrates the numerical performance of proposed method and in the last section concluding remarks are included.

2. Time-discontinuous Galerkin method (tDGM) for second order ODE

Consider a mass-spring-dashpot system as depicted in Figure 1. The governing equation of motion is described by the following second order initial value problem in time.

Figure 1.

Schematic diagram of the mass-spring-dashpot system.

$$\begin{array}{cc}\frac{d^{2}u}{{\mathit{dt}}^2}+2{\mathit{\zeta \omega }}_n\frac{\mathit{du}}{\mathit{dt}}+{\omega }_n^{2}u& =f\left(t\right)\forall t\in \left[0,T\right]\\ u\left(0\right)& =u_0\\ \frac{\mathit{du}\left(0\right)}{\mathit{dt}}& =v_0\end{array}$$E1

where $u≔u\left(t\right)$ is the unknown displacement, $f\left(t\right)$ is the external force acting on the system. Further, $u_0$ and $v_0$ are the prescribed initial values of the displacement and velocity, respectively. Damping ratio $\zeta$ and the natural frequency of vibration ${\omega }_n$ of the system are related to the mass $m$, stiffness of the spring $k$, and damping coefficient $c$ by:

$${\omega }_n=\sqrt{k/m},\zeta =\frac{c}{2m{\omega }_n}=\frac{c}{2\sqrt{\mathit{mk}}}$$E2

In what follows, this second order ODE will be utilized to discuss the fundamental concepts behind time-discontinuous Galerkin methods (henceforth, tDGM).

2.1 Two-field tDGM

In two-field tDGM (henceforth, uv-tDGM), both displacement ($u$) and velocity ($v$) are treated as independent primary variables and interpolated by using the piecewise polynomials. Both $u$ and $v$ are discontinuous at end-points (i.e., $t_n$ and $t_{n+1}$) of time-slab $I_n=\left(t_n,t_{n+1}\right)$. However, $u$ and $v$ remain continuous inside $I_n$, and approximated by piecewise polynomials (refer, Figure 2). Therefore, discontinuity occurs at discrete times belonging to a set $\left\{t_0,t_1,\cdots ,t_N\right\}$. The jump discontinuity in time for $u$ is denotes by

Figure 2.

Schematic diagram of time discontinuous approximation: (a) piecewise linear interpolation, and (b) piecewise quadratic interpolation.

$${〚u〛}_n=u_n^{+}-u_n^{-}$$E3

where

$$u_n^{+}=\underset{\epsilon \to 0}{lim}u\left(t+\epsilon \right),u_n^{-}=\underset{\epsilon \to 0}{lim}u\left(t-\epsilon \right)$$E4

are the discontinuous values of $u$ at time $t=t_n$. By recasting Eq. (1) into a system of two first-order ODEs one can obtain,

$$\frac{\mathit{dv}}{\mathit{dt}}+2{\mathit{\zeta \omega }}_{n}v+{\omega }_n^{2}u=f\left(t\right)\forall t\in \left[0,T\right]$$E5

$$\frac{\mathit{du}}{\mathit{dt}}-v=0\forall t\in \left[0,T\right]$$E6

$$u\left(0\right)=u_0,v\left(0\right)=v_0$$E7

The weak-form of the uv-tDGM can be stated as: find $u^h\in {\Im }_l^h$ and $v^h\in {\Im }_l^h$, such that for all $\delta u^h\in {\Im }_l^h$ and $\delta v^h\in {\Im }_l^h$, and for all $n=0,\cdots ,N-1$Eq. (8) holds.

$$\begin{array}{l}{\int }_{I_n}\delta v^h\left(\frac{{\mathit{dv}}^h}{\mathit{dt}}+2{\mathit{\zeta \omega }}_n{v}^h+{\omega }_n^2{u}^h-f\left(t\right)\right)\mathit{dt}+\delta v^h\left(t_n\right){\left[\left[v^h\right]\right]}_n\\ +{\int }_{I_n}\delta u^h\left(\frac{{\mathit{du}}^h}{\mathit{dt}}-v^h\right)\mathit{dt}+\delta u^h\left(t_n\right){\left[\left[u^h\right]\right]}_n=0\end{array}$$E8

where ${\Im }_l^h$ denotes the collection of polynomial with order less than or equal to $l$. It is worth noting that the presence ${\left[\left[u^h\right]\right]}_n$ and ${\left[\left[v^h\right]\right]}_n$ correspond to the weakly enforced initial condition for the $u$ and $v$, respectively. Further, since the selection of $\delta u^h$ and $\delta v^h$ is arbitrary one can depict Eq. (8) as,

$${\int }_{I_n}\delta v^h\left(\frac{{\mathit{dv}}^h}{\mathit{dt}}+2{\mathit{\zeta \omega }}_n{v}^h+{\omega }_n^2{u}^h-f\left(t\right)\right)\mathit{dt}+\delta v^h\left(t_n\right){\left[\left[v^h\right]\right]}_n=0$$E9

$${\int }_{I_n}\delta u^h\left(\frac{{\mathit{du}}^h}{\mathit{dt}}-v^h\right)\mathit{dt}+\delta u^h\left(t_n\right){\left[\left[u^h\right]\right]}_n=0$$E10

Eq. (10) denotes that, in uv-tDGM, displacement-velocity compatibility relationship is satisfied in weak form.

2.2 Single field tDGM

To decrease the number of unknowns in comparison to those involved in uv-tDGM, displacement-velocity compatibility condition (cf. Eq. 6) can be explicitly satisfied and velocity can be selected as primary unknown. Henceforth, this strategy will be termed as v-tDGM. In v-tDGM, $v$ is continuous in $I_n$, but discontinuity occurs at the end-points $t_n,t_{n+1}$. Further, $u$ is computed in a post-processing step by integration of $v$, therefore, $u$ remains continuous in time $\left[0,T\right]$.

The weak form of the v-tDGM reads: Find $v^h\in {\Im }_l^h$ such that for all $\delta v^h\in {\Im }_l^h$, and for all $n=0,\cdots ,N-1$Eq. (11) holds.

$${\int }_{I_n}\delta v^h\left(\frac{{\mathit{dv}}^h}{\mathit{dt}}+2{\mathit{\zeta \omega }}_n{v}^h+{\omega }_n^2{u}^h-f\left(t\right)\right)\mathit{dt}+\delta v^h\left(t_n\right){\left[\left[v^h\right]\right]}_n=0$$E11

Note that Eqs. (9) and (11) are identical, however, in former, $u^h$ is an independent variable and, in later, it is a dependent variable which will be computed by using following expression.

$$u^h\left(t\right)=u\left(t_n\right)+{\int }_{t_n}^t{v}^h\left(\tau \right)\mathit{d\tau }$$E12

Let us now focus on the discretization of weak-form (cf. Eq. (11)) by using the locally defined piecewise linear test and trial functions,

$$v^h=T_1{v}_n^{+}+T_2{v}_{n+1}^{-}\delta v^h=T_1\delta v_n^{+}+T_2\delta v_{n+1}^{-}$$E13

$$u^h\left(t\right)=u_n+v_n^{+}\frac{\mathrm{\Delta }{t}_n}{2}\left(1-T_1^2\right)+v_{n+1}^{-}\frac{\mathrm{\Delta }{t}_n}{2}{T}_2^2$$E14

where

$$T_1\left(\theta \right)=\frac{1-\theta }{2}{T}_2\left(\theta \right)=\frac{1+\theta }{2},\theta \in \left[-1,1\right].$$E15

Accordingly, Eq. 11 transforms into following matrix-vector form.

$$\begin{array}{l}\frac{1}{2}\left[\begin{array}{cc}1& 1\\ -1& 1\end{array}\right]\left\{\begin{array}{c}{v}_n^{-}\\ v_{n+1}^{+}\end{array}\right\}+\frac{2{\mathit{\zeta \omega }}_n\mathrm{\Delta }{t}_n}{6}\left[\begin{array}{cc}2& 1\\ 1& 2\end{array}\right]\left\{\begin{array}{c}{v}_n^{-}\\ v_{n+1}^{+}\end{array}\right\}\\ +\frac{{\omega }_n^2{\mathrm{\Delta }}_n^2}{24}\left[\begin{array}{cc}3& 1\\ 5& 3\end{array}\right]\left\{\begin{array}{c}{v}_n^{-}\\ v_{n+1}^{+}\end{array}\right\}=\left\{\begin{array}{c}{J}_{\mathit{ext}}^1\\ J_{\mathit{ext}}^2\end{array}\right\}-\frac{{\omega }_n^2\mathrm{\Delta }{t}_n{u}_n}{2}\left\{\begin{array}{l}1\\ 1\end{array}\right\}+\left\{\begin{array}{c}{v}_n^{-}\\ 0\end{array}\right\}\end{array}$$E16

where $J_{\mathit{ext}}^1$ and $J_{\mathit{ext}}^2$ are given by

$$J_{\mathit{ext}}^1={\int }_{I_n}{T}_{1}f\left(t\right)\mathit{dt}{J}_{\mathit{ext}}^2={\int }_{I_n}{T}_{2}f\left(t\right)\mathit{dt}$$

3. Numerical analysis of tDGM

In this section, numerical analysis of the tDGM schemes, (viz. uv-tDGM and v-tDGM) for the second order ODE will be performed. To assess the stability characteristics and temporal accuracy of these schemes, classical finite difference techniques will be used ([11], Chapter 9). In this context, it is sufficient to consider the following homogeneous and undamped form of Eq. (1):

$$\frac{d^{2}u}{{\mathit{dt}}^2}+{\omega }_n^{2}u=0,u\left(0\right)=u\left(0\right),\frac{\mathit{du}\left(0\right)}{\mathit{dt}}=v_0,t\in \left[0,T\right]$$E17

3.1 Energy decay in v-tDGM

In this section it will be shown that v-tDGM is a true energy-decaying scheme. Consider Eq. (17) which represents the governing equation of a spring-mass system. The total energy (sum of kinetic and potential energy) of the system remains constant because damping and external forces are absent in the system.

$$T_E\left(u,v\right)≔\frac{1}{2}{v}^2+\frac{1}{2}{\omega }_n^2{u}^2=\text{constant}$$E18

Consider the time domain $\left[0,T\right]$ and corresponding $N$ time-slabs; $I_n≔\left(t_n,t_{n+1}\right)$ for $n=0,1\cdots ,N-1$. Let the $u$ and $v$ at time $t_0=0$ be given by $u_0^{+}=u_0^{-}=u_0,$ and $v_0^{+}=v_0^{-}=v_0,$ respectively. Furthermore, the $u$ and $v$ at time $t_N=T$ are denoted by $u_N^{-}$ and $u_N^{-}$, respectively.

Accordingly, it can be shown that

$$T_E\left(u_N,v_N^{-}\right)=T_E\left(u_0,v_0\right)-\frac{1}{2}\sum _{n=0}^{N-1}{\left[\left[v^h\right]\right]}_n^2$$

or

$$T_E\left(u_N,v_N^{-}\right)\le T_E\left(u_0,v_0\right)$$

This shows that v-tDGM is an energy decaying time integration algorithm, in which the total energy during any time step, $T_E\left(u_N^{-},v_N^{-}\right)$, is always bounded from above by the total energy at the first time-step (i.e., $T_E\left(u_0,v_0\right)$).

To assess the energy dissipation characteristics of v-tDGM, Eq. (17) is solved with ${\omega }_n=2\pi$, $u_0=0$, and $v_0=1.0m/s$. The undamped time period $T_0$ of the sinusoidal motion is $1.0$ second, and the total time duration of simulation is $T=50$ seconds. Figure 3a depicts the time history graphs of the normalized total energy (i.e., $T_E\left(u,v\right)/T_E\left(u_0,v_0\right)$) computed by using v-tDGM with different time step sizes. Further, to visualize the effect of energy-dissipation displacement-velocity phase diagram is plotted in Figure 3b. For present problem, phase-diagram should be an ellipse. The presence of energy dissipation in the numerical algorithm, however, decreases the total energy which results in shortening of the radius of ellipse. From these plots it is evident that the dissipation of energy decreases as the time-step size decreases which also indicates that the jump discontinuity in time decreases with time-step size.

Figure 3.

Energy decay characteristics of v-tDGM; (a) temporal variation of normalized total energy and (b) phase diagram obtained with different time-step sizes.

3.2 Stability characteristics of v-tDGM

In this section, to study the stability characteristics of v-tDGM, Eq. (17) is considered. The matrix-vector form corresponding to this problem is given by

$$\begin{array}{l}\frac{1}{2}\left[\begin{array}{cc}1& 1\\ -1& 1\end{array}\right]\left\{\begin{array}{c}{v}_n^{-}\mathrm{\Delta }{t}_n\\ v_{n+1}^{+}\mathrm{\Delta }{t}_n\end{array}\right\}+\frac{{\mathrm{\Omega }}^2}{24}\left[\begin{array}{cc}3& 1\\ 5& 3\end{array}\right]\left\{\begin{array}{c}{v}_n^{-}\mathrm{\Delta }{t}_n\\ v_{n+1}^{+}\mathrm{\Delta }{t}_n\end{array}\right\}\\ =\left\{\begin{array}{c}{v}_n^{-}\mathrm{\Delta }{t}_n\\ 0\end{array}\right\}-\left\{\begin{array}{c}\frac{{\mathrm{\Omega }}^2{u}_n}{2}\\ \frac{{\mathrm{\Omega }}^2{u}_n}{2}\end{array}\right\},\end{array}$$E19

where $\mathrm{\Omega }={\omega }_n\mathrm{\Delta }{t}_n$. Subsequently, eliminating $v_n^{+}$ in Eq. (19),

$$\left\{\begin{array}{c}{u}_{n+1}\\ v_{n+1}^{-}\mathrm{\Delta }{t}_n\end{array}\right\}=\mathbf{A}\left(\mathrm{\Omega }\right)\left\{\begin{array}{c}{u}_n\\ v_n^{-}\mathrm{\Delta }{t}_n\end{array}\right\},$$E20

where $\mathbf{A}$ is the amplification matrix given by,

$$\mathbf{A}\left(\mathrm{\Omega }\right)=\left[\begin{array}{cc}\frac{{\mathrm{\Omega }}^4-30{\mathrm{\Omega }}^2+72}{{\mathrm{\Omega }}^4+6{\mathrm{\Omega }}^2+72}& \frac{-6{\mathrm{\Omega }}^2+72}{{\mathrm{\Omega }}^4+6{\mathrm{\Omega }}^2+72}\\ \frac{6{\mathrm{\Omega }}^4-72{\mathrm{\Omega }}^2}{{\mathrm{\Omega }}^4+6{\mathrm{\Omega }}^2+72}& \frac{-30{\mathrm{\Omega }}^2+72}{{\mathrm{\Omega }}^4+6{\mathrm{\Omega }}^2+72}\end{array}\right].$$E21

To investigate the stability of v-tDGM one should look into the eigenvalues of $\mathbf{A}$ (here, denoted by ${\lambda }_1$ and ${\lambda }_2$). Let the modulus of $\lambda$ be denoted by $\left|\lambda \right|=\sqrt{{\mathit{\lambda \lambda }}^{\ast }}$ with ${\lambda }^{\ast }$ denoting the complex conjugate of $\lambda$. Accordingly, the spectral radius of $\mathbf{A}$ can be described by $\rho \left(\mathbf{A}\right)={max}_{i=1,2}\left|{\lambda }_i\left(\mathbf{A}\right)\right|$. It can be easily shown that v-tDGM satisfies all criteria for the spectral stability [11, Chapter 9]: (a) $\rho \le 1$, (b) eigenvalues of $\mathbf{A}$ of multiplicity greater than one are strictly less than one in modulus. It proves that v-tDGM is an unconditionally stable time-marching scheme (for more details, readers are referred to [4]).

3.3 High-frequency response of TDG/FEM

Figure 4 plots the frequency responses of $\rho \left(\mathbf{A}\right)$ for v-tDGM. It is evident that $\rho \le 1$ which proves that present algorithm is unconditionally stable. The v-TDG/FEM, however, cannot attenuate spurious high-frequency contents since ${\rho }_{\infty }=1$ (see Figure 4). However, v-tDGM provides negligible attenuation in the small frequency regime as $\rho$ is close to one in this regime.

Figure 4.

Frequency response of spectral radius $\rho$ for v-tDGM.

3.4 Accuracy of v-tDGM

In [4], it is shown that $u$ in Eq. (20) satisfies the following finite difference stencil.

$$u_{n+1}-2a_1{u}_n+a_2{u}_{n-1}=0,$$E22

where $a_1=\mathit{\text{Trace}}\mathbf{A}/2$ and $a_2=\mathit{det}\mathbf{A}$. Let us now denote the exact solutions by $u\left(t\right)$ and $v\left(t\right)$. Then the local truncation error $\tau \left(t\right)$ corresponding to Eq. (22) at any time $t$ becomes

$$u\left(t+\mathrm{\Delta }t\right)-2a_{1}u\left(t\right)+a_{2}u\left(t-\mathrm{\Delta }t\right)=\mathrm{\Delta }{t}^2\tau \left(t\right)$$E23

Subsequently, by expanding $u\left(t+\mathrm{\Delta }t\right)$ and $u\left(t-\mathrm{\Delta }t\right)$ about $t$ by using Taylor series, and by using Eq. (17), it can be proved that v-tDGM is consistent and third order accurate, i.e., $\left|\tau \left(t\right)\right|⩽\frac{1}{72}\mathrm{\Delta }{t}^3$ [4]. Accordingly, one can use the Lax equivalence theorem to prove the convergence of the algorithms.

A direct consequence of the convergence is that the solution of Eq. (17) can be given by following expression [1]:

$$u_n=exp\left(-\frac{\overline{\zeta }\overline{\mathrm{\Omega }}{t}_n}{\mathrm{\Delta }t}\right)\left[k_{1}cos\left(\frac{\overline{\mathrm{\Omega }}{t}_n}{\mathrm{\Delta }t}\right)+k_{2}sin\left(\frac{\overline{\mathrm{\Omega }}{t}_n}{\mathrm{\Delta }t}\right)\right]$$E24

with

$$\overline{\mathrm{\Omega }}=arctan\left(\frac{\sqrt{a_2-a_1^2}}{a_1}\right)\overline{\zeta }=-\frac{1}{2\overline{\mathrm{\Omega }}}ln\left(a_2\right)$$E25

where $\overline{\zeta }$ denotes the algorithmic damping ratio, $\overline{\mathrm{\Omega }}$ is the frequency of the discrete solutions, and the coefficients $k_1$ and $k_2$ are determined by the displacement and velocity initial conditions.

Further, to investigate the accuracy of v-tDGM, algorithmic damping ratio, which is a measure of amplitude decay, and relative frequency error $\left(\mathrm{\Omega }-\overline{\mathrm{\Omega }}\right)/\overline{\mathrm{\Omega }}$, which is a measure of relative change in time period, are plotted in Figure 5. From Figure 5a it can be observed that $\zeta$ is comparable with the HHT-$\alpha$ scheme, however, it is significantly smaller than the uv-tDGM. It is evident that the Houbolt and Wilson-$\theta$ methods are too dissipative in the low-frequency range, therefore, these algorithms are not suitable for the long-duration numerical simulations. Furthermore, v-tDGM has smallest frequency error which can be attributed to its third order accuracy (refer, Figure 5b). It can be stated that these characteristics of v-tDGM, such as very low numerical dispersion and dissipation, third-order accuracy, and unconditional stability, make this scheme suitable for long-time simulations. However, at present, the only possible drawback to this method is its incapability to attenuate the spurious high-frequency components.

Figure 5.

Accuracy of v-tDGM: (a) algorithmic damping ratio, and (b) relative frequency error in low frequency regime (after [4]).

4. Statement of problem

A dam-reservoir-soil (DRS) system which is subjected to the spatially uniform horizontal ($a_1^g\left(t\right)$) and vertical ($a_2^g\left(t\right)$) component of ground motion is depicted in Figure 6. Reservoir domain contains linear, inviscid, irrotational, and compressible fluid and solid domain (dam and underlying soil) is treated as isotropic, homogeneous, linear elastic material. Computation domain of soil (${\mathrm{\Omega }}^s$) and fluid (${\mathrm{\Omega }}^f$) are obtained by prescribing the viscous boundary conditions at the artificial boundaries [10]. Let ${\mathrm{\Gamma }}_f^f$ and ${\mathrm{\Gamma }}_{\infty }^f$ be the free surface and upstream artificial boundary of fluid domain. ${\mathrm{\Gamma }}_{\mathit{fs}}^f$, ${\mathrm{\Gamma }}_{\mathit{fd}}^f$, ${\mathrm{\Gamma }}_{\mathit{fd}}^s$ and ${\mathrm{\Gamma }}_{\mathit{fs}}^s$ denote the fluid-soil, fluid-dam, dam-fluid and soil-fluid interfaces, respectively. Further, the outward unit normal vectors to the fluid and solid boundary are given by ${\mathbf{n}}^s$ and ${\mathbf{n}}^f$, respectively.

Figure 6.

Schematic diagram of dam-reservoir-soil (DRS) system subjected to seismic ground motion.

Further, hydrodynamic pressure distribution in the reservoir is modeled by the pressure wave equation,

$$\frac{1}{c^2}\frac{{\partial }^{2}p}{\partial t^2}-{\nabla }^{2}p=0\text{in}{\mathrm{\Omega }}^f\forall t\in \left(0,T\right),$$E26

with following initial and boundary conditions.

$$p\left(\mathbf{x},0\right)=0;\frac{\partial p\left(\mathbf{x},0\right)}{\partial t}=0\text{in}{\mathrm{\Omega }}^f\mathit{at}t=0$$E27

$$p\left(\mathbf{x},t\right)=0\mathrm{on}{\mathrm{\Gamma }}_f^f\forall t\in \left(0,T\right)$$E28

$$\nabla p\cdot {\mathbf{n}}^f=-{\rho }^f\frac{\partial \mathbf{v}}{\partial t}\cdot {\mathbf{n}}^f\mathrm{on}{\mathrm{\Gamma }}_{\mathit{fd}}^f\cup {\mathrm{\Gamma }}_{\mathit{fs}}^f\forall t\in \left(0,T\right)$$E29

$$\nabla p\cdot {\mathbf{n}}^f=-\frac{1}{c}\frac{\partial p}{\partial t}+\frac{1}{c}\frac{\partial p_{\infty }}{\partial t}\mathrm{on}{\mathrm{\Gamma }}_{\infty }^f\forall t\in \left(0,T\right)$$E30

In Eq. (26), ${\nabla }^2$ denotes the Laplace’s operator, $p\left(\mathbf{x},t\right)$ denotes the hydrodynamic pressure in the water (in excess of hydrostatic pressure) and $c$ denotes the speed of sound in water. Eq. (29) denotes the time dependent boundary condition at fluid-dam and fluid-soil interface, respectively, where ${\rho }^f$ is the mass density of fluid, and $\mathbf{v}$ is the velocity of solid domain (i.e., dam or soil). Eq. (30) is due to the viscous boundary condition at the upstream truncated boundary of reservoir. The first term in this equation corresponds to an array of dashpots placed normal to the truncated boundary ${\mathrm{\Gamma }}_{\infty }^f$, and the second term is due to the free-field response of reservoir.

Let us now consider the initial-boundary value problem of the solid domain which is described by,

$${\rho }^s\frac{{\partial }^2\mathbf{u}}{\partial t^2}-\nabla \cdot \mathit{s}-{\rho }^s\mathbf{b}=0\text{in}{\mathrm{\Omega }}^s\forall t\in \left(0,T\right),$$E31

$$\mathbf{u}\left(\mathbf{x},t\right)=\mathbf{g}\left(\mathbf{x},t\right)\mathrm{on}{\mathrm{\Gamma }}^g\forall t\in \times \left(0,T\right),$$E32

$$\mathit{\sigma }\cdot {\mathbf{n}}^s=\mathbf{h}\mathrm{on}{\mathrm{\Gamma }}_{\infty }^{+}\cup {\mathrm{\Gamma }}_{\mathit{fd}}^s\cup {\mathrm{\Gamma }}_{\mathit{fs}}^s\forall t\in \left(0,T\right),$$E33

$$\mathbf{u}\left(\mathbf{x},0\right)={\mathbf{u}}^0\left(\mathbf{x}\right),\frac{\partial \mathbf{u}}{\partial t}\left(\mathbf{x},0\right)={\mathbf{v}}^0\left(\mathbf{x}\right)\text{in}{\mathrm{\Omega }}^s\mathrm{at}t=0.$$E34

Furthermore, following time dependent boundary conditions will be considered in Eq. (33):

$$\mathit{\sigma }\cdot {\mathbf{n}}^s=-{\mathbf{c}}^v\left(\mathbf{v}-{\mathbf{v}}^{\infty }\right)+{\mathit{\sigma }}^{\infty }\cdot {\mathbf{n}}^s\mathrm{on}{\mathrm{\Gamma }}_{\infty }^L\cup {\mathrm{\Gamma }}_{\infty }^R$$E35

$$\mathit{\sigma }\cdot {\mathbf{n}}^s=-{\mathbf{c}}^h\cdot \mathbf{v}+2{\mathbf{c}}^h\cdot {\mathbf{v}}^{\mathit{\text{in}}}\mathrm{on}{\mathrm{\Gamma }}_{\infty }^B$$E36

$$\mathit{\sigma }\cdot {\mathbf{n}}^s=-\left\{p_h\left(\mathbf{x}\right)+p(\mathbf{x},t)\right\}{\mathbf{n}}^s\mathrm{on}{\mathrm{\Gamma }}_{\mathit{fd}}^s\cup {\mathrm{\Gamma }}_{\mathit{fs}}^s$$E37

In addition, solid domain is considered to be an isotropic, homogeneous, linear elastic material with

$${\sigma }_{\mathit{ij}}={\mathit{\lambda \epsilon }}_{\mathit{kk}}{\delta }_{\mathit{ij}}+2{\mathit{\mu \epsilon }}_{\mathit{ij}},$$E38

$${\epsilon }_{\mathit{ij}}=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right).$$E39

In Eqs. (31)–(34), ${\rho }^s$, $\mathbf{u}$, $\sigma$, $\mathbf{b}$, $\mathbf{g}$, $\mathbf{h}$, ${\mathbf{u}}^0$ and ${\mathbf{v}}^0$, denote mass density, displacement, Cauchy’s stress tensor (positive in tension), externally applied body force density, prescribed displacement, external surface traction, initial value of displacement and velocity, respectively. Eq. (37) represents time varying boundary condition due to the hydrostatic, $p_h\left(\mathbf{x}\right)$, and hydrodynamic, $p\left(\mathbf{x},t\right)$, pressure of impounded water acting on the dam-fluid and fluid-soil interface. Furthermore, In Eqs. (35)–(37), ${\mathbf{v}}^{\infty }$ and ${\sigma }^{\infty }$ are the velocity and the stress due to the free-field response of unbounded soil domain. In addition, ${\mathrm{\Gamma }}_{\infty }^L$, ${\mathrm{\Gamma }}_{\infty }^R$, and ${\mathrm{\Gamma }}_{\infty }^B$ represent left, right and bottom truncated boundaries of soil domain, respectively. In Eqs. (35) and (36), first term corresponds to the Lysmer and Kuhlemeyer viscous boundary condition. Physically, it represents a series of dashpots placed at the truncated boundaries of soil domain in parallel and normal directions (see Figure 6b). Further, these terms facilitate the absorption of outgoing scattered wave motion and attempt to model the radiation damping due to the semi-infinite soil domain. Where ${\mathbf{c}}^v$ and ${\mathbf{c}}^h$ are the damping coefficients matrices for the dashpots placed at vertical and horizontal truncated boundaries:

$${\mathbf{c}}^v=\left[\begin{array}{cc}{\rho }^s{c}_L& 0\\ 0& {\rho }^s{c}_T\end{array}\right],{\mathbf{c}}^h=\left[\begin{array}{cc}{\rho }^s{c}_T& 0\\ 0& {\rho }^s{c}_L\end{array}\right],$$E40

in which, $c_L$ and $c_T$ are the speed of longitudinal wave (P-wave) and transverse wave (S-wave) in the unbounded soil domain, respectively. Lastly, in Eqs. (38) and (39), $\lambda$ and $\mu$ are the Lame parameters, and ${\delta }_{\mathit{ij}}$ is the Kronecker delta function.

5. Space-time finite element method

Recently, ST/FEM is employed to solve dam-reservoir-soil interaction problem [10] in which the resultant matrix-vector form is given by

$$\left[{\mathbf{K}}_{11}^f\right]\left\{{\mathbf{Q}}_1\right\}+\left[{\mathbf{K}}_{12}^f\right]\left\{{\mathbf{Q}}_2\right\}-\frac{1}{2}\left[{\mathbf{H}}^f\right]\left\{{\mathbf{V}}_1\right\}+\frac{1}{2}\left[{\mathbf{H}}^f\right]\left\{{\mathbf{V}}_2\right\}=\left\{{\mathbf{J}}_1^f\right\}$$E41

$$\left[{\mathbf{K}}_{21}^f\right]\left\{{\mathbf{Q}}_1\right\}+\left[{\mathbf{K}}_{22}^f\right]\left\{{\mathbf{Q}}_2\right\}-\frac{1}{2}\left[{\mathbf{H}}^f\right]\left\{{\mathbf{V}}_1\right\}+\frac{1}{2}\left[{\mathbf{H}}^f\right]\left\{{\mathbf{V}}_2\right\}=\left\{{\mathbf{J}}_2^f\right\}$$E42

$$\left[{\mathbf{K}}_{11}^s\right]\left\{{\mathbf{V}}_1\right\}+\left[{\mathbf{K}}_{12}^s\right]\left\{{\mathbf{V}}_2\right\}+\frac{3\mathrm{\Delta }{t}_n^2}{24}\left[{\mathbf{H}}^s\right]\left\{{\mathbf{Q}}_1\right\}+\frac{\mathrm{\Delta }{t}_n^2}{24}\left[{\mathbf{H}}^s\right]\left\{{\mathbf{Q}}_2\right\}=\left\{{\mathbf{J}}_1^s\right\}$$E43

$$\left[{\mathbf{K}}_{21}^s\right]\left\{{\mathbf{V}}_1\right\}+\left[{\mathbf{K}}_{22}^s\right]\left\{{\mathbf{V}}_2\right\}+\frac{5\mathrm{\Delta }{t}_n^2}{24}\left[{\mathbf{H}}^s\right]\left\{{\mathbf{Q}}_1\right\}+\frac{3\mathrm{\Delta }{t}_n^2}{24}\left[{\mathbf{H}}^s\right]\left\{{\mathbf{Q}}_2\right\}=\left\{{\mathbf{J}}_2^s\right\}$$E44

where $\left\{{\mathbf{Q}}_1\right\}$ and $\left\{{\mathbf{Q}}_2\right\}$ represent the spatial nodal values of auxiliary variable $q=\partial p/\partial t$ at time $t_n^{+}$ and $t_{n+1}^{-}$, respectively. Similarly, $\left\{{\mathbf{V}}_1\right\}$ and $\left\{{\mathbf{V}}_2\right\}$ are the spatial nodal values of velocity field at time $t_n^{+}$ and $t_{n+1}^{-}$, respectively (for more details see [10]).

Further, In Eqs. (41) and (42),

$$\left[{\mathbf{K}}_{11}^f\right]=\frac{1}{2}\left[{\mathbf{M}}^f\right]+\frac{3\mathrm{\Delta }{t}_n^2}{24}\left[{\mathbf{K}}^f\right]+\frac{\mathrm{\Delta }{t}_n}{3}\left[{\mathbf{C}}_{\infty }^f\right],$$E45

$$\left[{\mathbf{K}}_{12}^f\right]=\frac{1}{2}\left[{\mathbf{M}}^f\right]+\frac{\mathrm{\Delta }{t}_n^2}{24}\left[{\mathbf{K}}^f\right]+\frac{\mathrm{\Delta }{t}_n}{6}\left[{\mathbf{C}}_{\infty }^f\right],$$E46

$$\left[{\mathbf{K}}_{21}^f\right]=-\frac{1}{2}\left[{\mathbf{M}}^f\right]+\frac{5\mathrm{\Delta }{t}_n^2}{24}\left[{\mathbf{K}}^f\right]+\frac{\mathrm{\Delta }{t}_n}{6}\left[{\mathbf{C}}_{\infty }^f\right],$$E47

$$\left[{\mathbf{K}}_{22}^f\right]=\frac{1}{2}\left[{\mathbf{M}}^f\right]+\frac{3\mathrm{\Delta }{t}_n^2}{24}\left[{\mathbf{K}}^f\right]+\frac{\mathrm{\Delta }{t}_n}{3}\left[{\mathbf{C}}_{\infty }^f\right],$$E48

where

$$\left[{\mathbf{M}}^f\right]={\int }_{{\mathrm{\Omega }}^f}{\mathbf{N}}^T\mathbf{N}\frac{1}{c^2}d\mathrm{\Omega }$$E49

$$\left[{\mathbf{K}}^f\right]={\int }_{{\mathrm{\Omega }}^f}\nabla \mathbf{N}\cdot \nabla {\mathbf{N}}^{T}d\mathrm{\Omega }$$E50

$$\left[{\mathbf{C}}_{\infty }^f\right]={\int }_{{\mathrm{\Gamma }}_{\infty }^f}{\mathbf{N}}^T\mathbf{N}\frac{1}{c}\mathit{ds}$$E51

denote mass matrix, diffusion matrix, and viscous boundary at the upstream truncated boundary, respectively. The fluid-solid coupling matrix $\left[{\mathbf{H}}^f\right]$ is given by,

$$\left[{\mathbf{H}}^f\right]={\int }_{{\mathrm{\Gamma }}_{\mathit{fs}}^f\cup {\mathrm{\Gamma }}_{\mathit{fd}}^f}{\mathbf{N}}^T\mathbf{N}{\rho }^f{\mathbf{n}}^f\mathit{ds}$$E52

and right hand side spatial nodal vectors in Eqs. (41) and (42) becomes,

$$\left\{{\mathbf{J}}_1^f\right\}=\left[{\mathbf{M}}^f\right]\left\{{\mathbf{Q}}_0\right\}-\frac{\mathrm{\Delta }{t}_n}{2}\left[{\mathbf{K}}^f\right]\left\{{\mathbf{P}}_0\right\}+\frac{\mathrm{\Delta }{t}_n}{3}\left[{\mathbf{C}}_{\infty }^f\right]\left\{{\mathbf{Q}}_{\mathbf{1}}^{\infty }\right\}+\frac{\mathrm{\Delta }{t}_n}{6}\left[{\mathbf{C}}_{\infty }^f\right]\left\{{\mathbf{Q}}_{\mathbf{2}}^{\infty }\right\}$$E53

$$\left\{{\mathbf{J}}_2^f\right\}=-\frac{\mathrm{\Delta }{t}_n}{2}\left[{\mathbf{K}}^f\right]\left\{{\mathbf{P}}_0\right\}+\frac{\mathrm{\Delta }{t}_n}{6}\left[{\mathbf{C}}_{\infty }^f\right]\left\{{\mathbf{Q}}_{\mathbf{1}}^{\infty }\right\}+\frac{\mathrm{\Delta }{t}_n}{3}\left[{\mathbf{C}}_{\infty }^f\right]\left\{{\mathbf{Q}}_{\mathbf{2}}^{\infty }\right\}$$E54

where $\left\{{\mathbf{Q}}_0\right\}$ and $\left\{{\mathbf{P}}_0\right\}$ correspond to the nodal values of $q$ and $p$ at time $t_n$, and $\left\{{\mathbf{Q}}_1^{\infty }\right\}$ and $\left\{{\mathbf{Q}}_2^{\infty }\right\}$ are the free field hydrodynamic response of reservoir at time $t_n$ and $t_{n+1}$, respectively.

In Eqs. (43) and (44),

$$\left[{\mathbf{K}}_{11}^s\right]=\frac{1}{2}\left[{\mathbf{M}}^s\right]+\frac{3\mathrm{\Delta }{t}_n^2}{24}\left[{\mathbf{K}}^s\right]+\frac{\alpha \mathrm{\Delta }{t}_n}{3}\left[{\mathbf{M}}^s\right]+\frac{\beta \mathrm{\Delta }{t}_n}{3}\left[{\mathbf{K}}^s\right]+\frac{\mathrm{\Delta }{t}_n}{3}\left[{\mathbf{C}}_{\infty }^s\right],$$E55

$$\left[{\mathbf{K}}_{12}^s\right]=\frac{1}{2}\left[{\mathbf{M}}^s\right]+\frac{\mathrm{\Delta }{t}_n^2}{24}\left[{\mathbf{K}}^s\right]+\frac{\alpha \mathrm{\Delta }{t}_n}{6}\left[{\mathbf{M}}^s\right]+\frac{\beta \mathrm{\Delta }{t}_n}{6}\left[{\mathbf{K}}^s\right]+\frac{\mathrm{\Delta }{t}_n}{6}\left[{\mathbf{C}}_{\infty }^s\right],$$E56

$$\left[{\mathbf{K}}_{21}^s\right]=-\frac{1}{2}\left[{\mathbf{M}}^s\right]+\frac{5\mathrm{\Delta }{t}_n^2}{24}\left[{\mathbf{K}}^s\right]+\frac{\alpha \mathrm{\Delta }{t}_n}{6}\left[{\mathbf{M}}^s\right]+\frac{\beta \mathrm{\Delta }{t}_n}{6}\left[{\mathbf{K}}^s\right]+\frac{\mathrm{\Delta }{t}_n}{6}\left[{\mathbf{C}}_{\infty }^s\right],$$E57

$$\left[{\mathbf{K}}_{22}^s\right]=\frac{1}{2}\left[{\mathbf{M}}^s\right]+\frac{3\mathrm{\Delta }{t}_n^2}{24}\left[{\mathbf{K}}^s\right]+\frac{\alpha \mathrm{\Delta }{t}_n}{3}\left[{\mathbf{M}}^s\right]+\frac{\beta \mathrm{\Delta }{t}_n}{3}\left[{\mathbf{K}}^s\right]+\frac{\mathrm{\Delta }{t}_n}{3}\left[{\mathbf{C}}_{\infty }^s\right],$$E58

in which $\left[{\mathbf{M}}^s\right]$, and ${\left[\mathbf{K}\right]}^s$ are the mass and stiffness matrix for the solid domain [12], $\alpha$ and $\beta$ are the coefficients of Rayleigh damping, and matrix $\left[{\mathbf{C}}_{\infty }^s\right]$ is due to the dashpots placed at truncated boundaries of soil domain which has the form,

$$\left[{\mathbf{C}}_{\infty }^s\right]=\left[{\mathbf{C}}_v^s\right]+\left[{\mathbf{C}}_h^s\right]=\left[\begin{array}{cc}{\mathbf{c}}_{11}^v& \mathbf{0}\\ \mathbf{0}& {\mathbf{c}}_{22}^v\end{array}\right]+\left[\begin{array}{cc}{\mathbf{c}}_{11}^h& \mathbf{0}\\ 0& {\mathbf{c}}_{22}^h\end{array}\right]$$E59

$$\left[{\mathbf{c}}_{\mathit{ii}}^v\right]={\int }_{{\mathrm{\Gamma }}_{\infty }^L\cup {\mathrm{\Gamma }}_{\infty }^R}{c}_{\mathit{ii}}^v{\mathbf{N}}^T\mathbf{N}\mathit{ds}i=1,2\left(\mathrm{no}\mathrm{sum}\right),$$E60

$$\left[{\mathbf{c}}_{\mathit{ii}}^h\right]={\int }_{{\mathrm{\Gamma }}_{\infty }^B}{c}_{\mathit{ii}}^h{\mathbf{N}}^T\mathbf{N}\mathit{ds}i=1,2\left(\mathrm{no}\mathrm{sum}\right),$$E61

the solid-fluid coupling matrix,

$$\left[{\mathbf{H}}^s\right]={\int }_{{\mathrm{\Gamma }}_{\mathit{fs}}^s\cup {\mathrm{\Gamma }}_{\mathit{fd}}^s}{\mathbf{N}}^T\mathbf{N}{\mathbf{n}}^s\mathit{ds},$$E62

and right hand side spatial nodal vectors is given by,

$$\begin{array}{l}\left\{{\mathbf{J}}_1^s\right\}=\left[{\mathbf{M}}^s\right]\left\{{\mathbf{V}}_0\right\}-\frac{\mathrm{\Delta }{t}_n}{2}\left[{\mathbf{K}}^s\right]\left\{{\mathbf{U}}_0\right\}+\frac{2\mathrm{\Delta }{t}_n}{3}\left[{\mathbf{C}}_h^s\right]\left\{{\mathbf{V}}_1^{\mathit{\text{in}}}\right\}+\frac{2\mathrm{\Delta }{t}_n}{6}\left[{\mathbf{C}}_h^s\right]\left\{{\mathbf{V}}_2^{\mathit{\text{in}}}\right\}\\ +\frac{\mathrm{\Delta }{t}_n}{3}\left[{\mathbf{C}}_v^s\right]\left\{{\mathbf{V}}_1^{\infty }\right\}+\frac{\mathrm{\Delta }{t}_n}{6}\left[{\mathbf{C}}_v^s\right]\left\{{\mathbf{V}}_2^{\infty }\right\}+\frac{\mathrm{\Delta }{t}_n}{2}\left\{{\mathbf{F}}_1^{\mathit{ext}}\right\}+\frac{\mathrm{\Delta }{t}_n}{2}\left\{{\mathbf{F}}_1^{\infty }\right\}\\ -\frac{\mathrm{\Delta }{t}_n}{2}\left[{\mathbf{H}}^s\right]\left(\left\{{\mathbf{P}}_{\mathbf{h}}\right\}+\left\{{\mathbf{P}}_{\mathbf{0}}\right\}\right),\end{array}$$E63

$$\begin{array}{c}\left\{{\mathbf{J}}_2^s\right\}=-\frac{\mathrm{\Delta }{t}_n}{2}\left[{\mathbf{K}}^s\right]\left\{{\mathbf{U}}_0\right\}+\frac{2\mathrm{\Delta }{t}_n}{6}\left[{\mathbf{C}}_h^s\right]\left\{{\mathbf{V}}_1^{\mathit{\text{in}}}\right\}+\frac{2\mathrm{\Delta }{t}_n}{3}\left[{\mathbf{C}}_h^s\right]\left\{{\mathbf{V}}_2^{\mathit{\text{in}}}\right\}+\frac{\mathrm{\Delta }{t}_n}{6}\left[{\mathbf{C}}_v^s\right]\left\{{\mathbf{V}}_1^{\infty }\right\}\\ +\frac{\mathrm{\Delta }{t}_n}{3}\left[{\mathbf{C}}_v^s\right]\left\{{\mathbf{V}}_2^{\infty }\right\}+\frac{\mathrm{\Delta }{t}_n}{2}\left\{{\mathbf{F}}_2^{\mathit{ext}}\right\}+\frac{\mathrm{\Delta }{t}_n}{2}\left\{{\mathbf{F}}_2^{\infty }\right\}-\frac{\mathrm{\Delta }{t}_n}{2}\left[{\mathbf{H}}^s\right]\left(\left\{{\mathbf{P}}_{\mathbf{h}}\right\}+\left\{{\mathbf{P}}_{\mathbf{0}}\right\}\right),\end{array}$$E64

where $\left\{{\mathbf{U}}_0\right\}$, $\left\{{\mathbf{V}}_0\right\}$, and $\left\{{\mathbf{P}}_0\right\}$ are spatial nodal values of $\mathbf{u}$, $\mathbf{v}$, and $p$ at time $t_n$, $\left\{{\mathbf{P}}_h\right\}$ is nodal values of hydrostatic pressure, and $\left\{{\mathbf{F}}_1^{\mathit{ext}}\right\}$ and $\left\{{\mathbf{F}}_2^{\mathit{ext}}\right\}$ are nodal force vector due to external body forces at time $t_n$ and $t_{n+1}$, respectively,

$$\left\{{\mathbf{F}}_a^{\mathit{ext}}\right\}={\int }_{-1}^{+1}{\int }_{{\mathrm{\Omega }}_h^s}{T}_a\mathbf{N}{\rho }^s\mathbf{b}d\mathrm{\Omega }\mathit{d\theta }.$$E65

Further, the force vector $\left\{{\mathbf{F}}_{a=1,2}^{\infty }\right\}$ is due to the free-field stress at the vertical viscous boundaries of soil domain [13], and described by,

$$\left\{{\mathbf{F}}_{a=1,2}^{\infty }\right\}={\int }_{-1}^{+1}{\int }_{{\mathrm{\Gamma }}_{\infty }^L\cup {\mathrm{\Gamma }}_{\infty }^R}{T}_a\mathbf{N}{\mathit{\sigma }}^{\infty }\cdot {\mathbf{n}}^{\mathbf{s}}\mathit{ds}\mathit{d\theta }.$$E66

Lastly, nodal values of displacement and pressure field at time $t_{n+1}$ are computed by following expression in a post-processing step.

$$\left\{{\mathbf{U}}_2\right\}=\left\{{\mathbf{U}}_0\right\}+\frac{\mathrm{\Delta }{t}_n}{2}\left(\left\{{\mathbf{V}}_1\right\}+\left\{{\mathbf{V}}_2\right\}\right)$$E67

$$\left\{{\mathbf{P}}_2\right\}=\left\{{\mathbf{P}}_0\right\}+\frac{\mathrm{\Delta }{t}_n}{2}\left(\left\{{\mathbf{Q}}_1\right\}+\left\{{\mathbf{Q}}_2\right\}\right)$$E68

6. Numerical examples

In this section, ST/FEM with block iterative algorithm has been employed to study the response of the concrete gravity dam to the horizontal earthquake motion (see [10]). In the numerical modeling two cases are considered; (i) dam-reservoir (DR) system, in which foundation is considered to be rigid, and (ii) dam-reservoir-soil (DRS) system, in which the foundation is an elastic deformable body.

Figure 7 depicts the physical dimensions of the dam-reservoir system. Length of the reservoir in upstream direction is $200$ m, and length of the soil domain in horizontal and vertical direction is $440$ m and $150$ m, respectively. For the dam, elastic modulus, $E$, mass-density, ${\rho }^s$, and Poisson’s ratio, $\nu$, are $28.0$ GPa, $2347.0$ kg/m³, and $0.20$ respectively, and for the foundation, $E=40.0$ GPa, $\rho =2551.0$ kg/m³, and $\nu =0.20$. Material damping in the solid domain is modeled by Rayleigh damping with $\xi =5\%$ viscous damping specified for the foundation and dam separately.

Figure 7.

Physical dimensions of dam-reservoir system (after [10]).

Figure 8 represents the accelerogram (horizontal component) recorded at a control point on the free surface; the maximum and minimum values of acceleration are $396.7$ Gal (at $t=15.19$ s) and $-449.6$ Gal (at $t=14.82$ s), respectively. Further, numerical simulations are performed for a total time duration of $45$ s with a uniform time step size $\mathrm{\Delta }t=0.01$. Time history graphs of acceleration at the crest of the dam in DR and DRS systems are plotted in Figure 9 where it can be seen that the deformation characteristics of underlying foundation significantly decreases the responses for dam. To this end, absolute maximum value of horizontal and vertical component of acceleration obtained for DRS are $1489.89$ Gal and $597.47$ Gal, respectively, and for the DR system these values are equal to $3897.10$ Gal and $1274.65$ Gal. In addition, Fourier spectrum of time-series of acceleration at the crest indicates that in the case of DRS there is a significant decay in the amplitudes and an elongation of time period as compare to the DR system.

Figure 8.

Time history of horizontal component of ground motion recorded at free-surface.

Figure 9.

Acceleration response at the crest of dam in DR and DRS system; time history of (a) horizontal component and (b) vertical component of acceleration, and Fourier spectrum of (c) horizontal component and (d) vertical component of acceleration (after [10]).

Interestingly, in both cases, it is observed that the critical location for pressure is at the base of the dam. Figure 10 presents the evolution of $p$, maximum principle tensile and compressive stresses with time at this location. It is clearly visible that dynamic interactions between the dam-reservoir and the deformable underlying ground significantly lower the hydrodynamic pressure and maximum stresses in the dam. Lastly, the hydrodynamic pressure field and deformed configuration of dam (magnified 500 times) in the DRS system at time $t$ = 18.02 s and $t$ = 18.08 s are presented in Figure 11.

Figure 10.

Temporal response of (a) normalized hydrodynamic pressure, (b) principal tensile stress and (c) principal compressive stress at the base of dam in DR and DRS system (after [10]).

Figure 11.

(a) and (c): Hydrodynamic pressure field in the reservoir, and (b) and (d): magnified deformed configuration of dam at times t = 18:02 s and t = 18:08 s (after [10]).

7. Conclusions

In this chapter, novel concepts of time-discontinuous Galerkin (tDGM) method is presented. A method called v-tDGM is derived to solve second order ODEs in time. In this method velocity is the primary unknown and it remain discontinuous at discrete times. Thereby, the time-continuity of velocity is satisfied in a weak sense. However, displacement is obtained by time-integration of the velocity in a post-processing step by virtue of which it is continuous in time. It is demonstrated that the present method is unconditionally stable and third-order accurate in time for linear interpolation of velocity in time. Therefore, it can be stated that the numerical characteristics of the v-ST/FEM scheme, therefore, make it highly suitable for computing the response of bodies subjected to dynamic loading conditions, such as fast-moving loads, impulsive loading, and long-duration seismic loading, among others.

Subsequently, ST/FEM is used to compute the response of a dam-reservoir-soil (DRS) system to the earthquake loading while considering all types of dynamic interactions. An auxiliary variable $q$ representing the first order time derivative of the pressure is treated as the primary unknown for the reservoir domain. Similarly, velocity $\mathbf{v}$ is the primary unknown for the solid domain. Both $\mathbf{v}$ and $q$ are interpolated such that they remain discontinuous at the discrete times. Hydrodynamic pressure and displacements are the secondary unknowns in the present formulation which are computed by time integration of $q$ and $\mathbf{v}$, respectively. It is concluded that the dynamic interactions between the dam-reservoir system and the underlying deformable foundation significantly dampen the seismic response of the dam and the reservoir and elongate the time period of the acceleration response of the dam.