Abstract

In this paper, a finite element (FE)-based model that efficiently evaluates the dynamic behavior of dam-reservoir-foundation interaction (DRFI) problem was proposed including the radiation of waves to the unbounded rock and reservoir domains. Lagrangian fluid elements were used to discretize the near-field reservoir domain, while the presented infinite fluid elements were used to discretize the far-field reservoir domain. The fully coupled equation of motion for DRFI problem was solved by direct method. A two-dimensional (2D) plane-strain FE formulation of the problem is written in FORTRAN 90 programming language. Investigations were conducted on the effect of near-field domain size (length and depth) on the dynamic behavior of DRFI, dam-foundation interaction (DFI), and dam-reservoir interaction (DRI) problems. The results of this study demonstrate that the proposed model outperforms many other models that have been evolved in the literature in terms of accuracy and speed. The reflected hydrodynamic pressures at the far-field reservoir domain were efficiently absorbed by the suggested infinite fluid elements. The near-field domains size has a noticeable impact on the dynamic behavior of the dam. Making an exact choice about the size is more challenging. However, it was observed that the size of 1.5H is the physically appropriate response.

Keywords
Dam-reservoir-foundation interaction; numerical simulation; near-field domain size; new infinite fluid elements; dynamic response

Graphical Abstract

1 INTRODUCTION

Dams are huge structures that play an important role in human life, providing for many of their essential needs, such as power generation, irrigation, flood prevention, regulating river flow for shipping, etc. In the past studies, the dynamic analysis of DRI, DFI, and DRFI problems was studied by many researchers using several approaches such as the FE method (Chopra and Chakrabarti 1971Chopra, A. and P. Chakrabarti (1971). "The Koyna Earthquake of December 11, 1967 and the Performance of Koyna Dam, Report No. EERC 71-1." Earthquake Engineering Research Center, University of California, Berkeley., Chopra and Chakrabarti 1973Chopra, A. K. and P. Chakrabarti (1973). "The Koyna earthquake and the damage to Koyna dam." Bulletin of the Seismological Society of America 63(2): 381-397., Gutierrez and Chopra 1978Gutierrez, J. A. and A. K. Chopra (1978). "A substructure method for earthquake analysis of structures including structure‐soil interaction." Earthquake Engineering and Structural Dynamics 6(1): 51-69., Akkas, Akay et al. 1979Akkas, N., H. Akay and C. Yilmaz (1979). "Applicability of general-purpose finite element programs in solid-fluid interaction problems." Computers and Structures 10(5): 773-783., Hacıefendioğlu 2009Hacıefendioğlu, K. (2009). "Stochastic response of concrete faced rockfill dams including partially ice-covered reservoir–foundation interaction under spatially varying seismic waves." Cold regions science and technology 58(1-2): 57-67., Akköse and Şimşek 2010Akköse, M. and E. Şimşek (2010). "Non-linear seismic response of concrete gravity dams to near-fault ground motions including dam-water-sediment-foundation interaction." Applied Mathematical Modelling 34(11): 3685-3700., Bayraktar, Sevim et al. 2011Bayraktar, A., B. Sevim and A. C. Altunişik (2011). "Finite element model updating effects on nonlinear seismic response of arch dam–reservoir–foundation systems." Finite elements in analysis and design 47(2): 85-97., Gorai and Maity 2021Gorai, S. and D. Maity (2021). "Numerical investigation on seismic behaviour of aged concrete gravity dams to near source and far source ground motions." Natural Hazards 105(1): 943-966.), Boundary Element (BE) method (Medina and Domínguez 1989Medina, F. and J. Domínguez (1989). "Boundary elements for the analysis of the seismic response of dams including dam-water-foundation interaction effects. I." Engineering Analysis with Boundary Elements 6(3): 152-157., Dominguez and Meise 1991Dominguez, J. and T. Meise (1991). "On the use of the BEM for wave propagation in infinite domains." Engineering Analysis with Boundary Elements 8(3): 132-138., Tsai, Lee et al. 1992Tsai, C., G. Lee and C. Yeh (1992). "Time-domain analyses of three-dimensional dam-reservoir interactions by BEM and semi-analytical method." Engineering Analysis with Boundary Elements 10(2): 107-118., Maeso and Domínguez 1993Maeso, O. and J. Domínguez (1993). "Earthquake analysis of arch dams. I: Dam-foundation interaction." Journal of engineering mechanics 119(3): 496-512., Küçükarslan 2004aKüçükarslan, S. (2004a). "Dynamic analysis of dam–reservoir-foundation interaction in time domain." Computational Mechanics 33(4): 274-281., Küçükarslan 2004bKüçükarslan, S. (2004b). "Time‐domain dynamic analysis of dam–reservoir–foundation interaction including the reservoir bottom absorption." International journal for numerical and analytical methods in geomechanics 28(9): 963-980.), and coupled FE-BE method (Touhei and Ohmachi 1994Touhei, T. and T. Ohmachi (1994). "Modal analysis of a dam–foundation system based on an FE–BE method in the time domain." Earthquake engineering and structural dynamics 23(1): 1-15., Abouseeda and Dakoulas 1998Abouseeda, H. and P. Dakoulas (1998). "Non‐linear dynamic earth dam–foundation interaction using a BE–FE method." Earthquake engineering and structural dynamics 27(9): 917-936., Yaseri and Konrad 2021Yaseri, A. and J. M. Konrad (2021). "Computation of amplification functions of earth dam‐flexible canyon systems by the hybrid FEM‐SBFEM technique." Earthquake Engineering and Structural Dynamics 50(11): 2883-2896.). The modeling of reservoir domain which is a challenging topic has been implemented by using the approaches such as Lagrangian (Akkas, Akay et al. 1979Akkas, N., H. Akay and C. Yilmaz (1979). "Applicability of general-purpose finite element programs in solid-fluid interaction problems." Computers and Structures 10(5): 773-783., Wilson and Khalvati 1983Wilson, E. L. and M. Khalvati (1983). "Finite elements for the dynamic analysis of fluid‐solid systems." International Journal for Numerical Methods in Engineering 19(11): 1657-1668., Hacıefendioğlu 2009Hacıefendioğlu, K. (2009). "Stochastic response of concrete faced rockfill dams including partially ice-covered reservoir–foundation interaction under spatially varying seismic waves." Cold regions science and technology 58(1-2): 57-67., Akköse and Şimşek 2010Akköse, M. and E. Şimşek (2010). "Non-linear seismic response of concrete gravity dams to near-fault ground motions including dam-water-sediment-foundation interaction." Applied Mathematical Modelling 34(11): 3685-3700., Bayraktar, Sevim et al. 2011Bayraktar, A., B. Sevim and A. C. Altunişik (2011). "Finite element model updating effects on nonlinear seismic response of arch dam–reservoir–foundation systems." Finite elements in analysis and design 47(2): 85-97., Hacıefendioğlu, Soyluk et al. 2012Hacıefendioğlu, K., K. Soyluk and F. Birinci (2012). "Numerical investigation of stochastic response of an elevated water tank to random underground blast loading." Stochastic environmental research and risk assessment 26(4): 599-607.), Eulerian (Chopra and Chakrabarti 1971Chopra, A. and P. Chakrabarti (1971). "The Koyna Earthquake of December 11, 1967 and the Performance of Koyna Dam, Report No. EERC 71-1." Earthquake Engineering Research Center, University of California, Berkeley., Chopra and Chakrabarti 1973Chopra, A. K. and P. Chakrabarti (1973). "The Koyna earthquake and the damage to Koyna dam." Bulletin of the Seismological Society of America 63(2): 381-397., Chopra 2020Chopra, A. K. (2020). Earthquake Engineering for Concrete Dams: Analysis, Design, and Evaluation, John Wiley and Sons., Gorai and Maity 2021Gorai, S. and D. Maity (2021). "Numerical investigation on seismic behaviour of aged concrete gravity dams to near source and far source ground motions." Natural Hazards 105(1): 943-966.), and Added Mass (Westergaard 1933Westergaard, H. M. (1933). "Water pressures on dams during earthquakes." Transactions of the American society of Civil Engineers 98(2): 418-433., Kuo 1982Kuo, J. S.-H. (1982). Fluid-structure interactions: Added mass computations for incompressible fluid, University of California, College of Engineering, Earthquake Engineering …., Zhang, Wang et al. 2013Zhang, S., G. Wang, B. Pang and C. Du (2013). "The effects of strong motion duration on the dynamic response and accumulated damage of concrete gravity dams." Soil Dynamics and Earthquake Engineering 45: 112-124.) with appropriate boundary conditions (Sommerfeld 1949Sommerfeld, A. (1949). Partial differential equations in physics, Academic press., Sharan 1985Sharan, S. K. (1985). "Finite element analysis of unbounded and incompressible fluid domains." International Journal for Numerical Methods in Engineering 21(9): 1659-1669., Higdon 1991Higdon, R. L. (1991). "Absorbing boundary conditions for elastic waves." Geophysics 56(2): 231-241., Kellezi 2000Kellezi, L. (2000). "Local transmitting boundaries for transient elastic analysis." Soil Dynamics and Earthquake Engineering 19(7): 533-547., Pelecanos, Kontoe et al. 2013Pelecanos, L., S. Kontoe and L. Zdravković (2013). "Numerical modelling of hydrodynamic pressures on dams." Computers and Geotechnics 53: 68-82., Khazaee and Lotfi 2014Khazaee, A. and V. Lotfi (2014). "Application of perfectly matched layers in the transient analysis of dam–reservoir systems." Soil Dynamics and Earthquake Engineering 60: 51-68.). The dynamic behavior of dam structure is generally evaluated either in the time domain (Chopra 1967Chopra, A. K. (1967). "Reservoir-dam interaction during earthquakes." Bulletin of the Seismological Society of America 57(4): 675-687., Chopra 1968Chopra, A. K. (1968). "Earthquake behavior of reservoir-dam systems." Journal of the Engineering Mechanics Division 94(6): 1475-1500., Nath 1971Nath, B. (1971). "Coupled hydrodynamic response of a gravity dam " Proceedings of the Institution of Civil Engineers 48(2): 245-257., Finn and Varoḡlu 1973Finn, W. L. and E. Varoḡlu (1973). "Dynamics of gravity dam-reservoir systems." Computers and Structures 3(4): 913-924., Hall 1986Hall, J. F. (1986). "Study of the earthquake response of Pine Flat dam." Earthquake engineering and structural dynamics 14(2): 281-295., Guan, Moore et al. 1994Guan, F., I. D. Moore and G. Lin (1994). "Transient response of reservoir–dam–soil systems to earthquakes." International journal for numerical and analytical methods in geomechanics 18(12): 863-880., Batta and Pekau 1996Batta, V. and O. Pekau (1996). "Application of boundary element analysis for multiple seismic cracking in concrete gravity dams." Earthquake engineering and structural dynamics 25(1): 15-30., Samii and Lotfi 2007Samii, A. and V. Lotfi (2007). "Comparison of coupled and decoupled modal approaches in seismic analysis of concrete gravity dams in time domain." Finite elements in analysis and design 43(13): 1003-1012., Burman, Nayak et al. 2012Burman, A., P. Nayak, P. Agrawal and D. Maity (2012). "Coupled gravity dam–foundation analysis using a simplified direct method of soil–structure interaction." Soil Dynamics and Earthquake engineering 34(1): 62-68.) or frequency domain (Mei, Foda et al. 1979Mei, C. C., M. A. Foda and P. Tong (1979). "Exact and hybrid-element solutions for the vibration of a thin elastic structure seated on the sea floor." Applied ocean research 1(2): 79-88., Fenves and Chopra 1983Fenves, G. and A. K. Chopra (1983). "Effects of reservoir bottom absorption on earthquake response of concrete gravity dams." Earthquake engineering and structural dynamics 11(6): 809-829., Fenves and Chopra 1987Fenves, G. and A. K. Chopra (1987). "Simplified earthquake analysis of concrete gravity dams." Journal of engineering mechanics 113(8): 1688-1708., Lotfi, Roesset et al. 1987Lotfi, V., J. M. Roesset and J. L. Tassoulas (1987). "A technique for the analysis of the response of dams to earthquakes." Earthquake Engineering and Structural Dynamics 15(4): 463-489., Humar and Jablonski 1988Humar, J. and A. Jablonski (1988). "Boundary element reservoir model for seismic analysis of gravity dams." Earthquake engineering and structural dynamics 16(8): 1129-1156., Domínguez and Medina 1989Domínguez, J. and F. Medina (1989). "Boundary elements for the analysis of the seismic response of dams including dam-water-foundation interaction effects. II." Engineering Analysis with Boundary Elements 6(3): 158-163., Medina and Domínguez 1989Medina, F. and J. Domínguez (1989). "Boundary elements for the analysis of the seismic response of dams including dam-water-foundation interaction effects. I." Engineering Analysis with Boundary Elements 6(3): 152-157., Bouaanani, Paultre et al. 2002Bouaanani, N., P. Paultre and J. Proulx (2002). "Two‐dimensional modelling of ice cover effects for the dynamic analysis of concrete gravity dams." Earthquake engineering and structural dynamics 31(12): 2083-2102.) and it is essentially affected by important factors like the interaction of the dam and water with bedrock, radiation of waves to the infinity, the size of near-field domains in the water and rock foundation, the mass of rock foundation, the compressibility of water and hydrodynamic pressures. Thus, the inclusion of these parameters in the modeling of the problem is important.

In the modeling of the DRI, DRFI problems using Lagrangian aproach, the coupled equation of motion is generally solved by the methods such as Mode superposition, Newmark, Wilson-θ and HHT-α method. The far-field reservoir domain effect as well as foundation effect is not included in the FE modeling of the problem as radiation damping. Furthermore, the material damping mechanism of the system is generally selected as viscous damping (Bayraktar, Akkose et al. 2005aBayraktar, A., M. Akkose and S. Adanur (2005a). "The effects of sediment on earthquake response of concrete gravity dams to asynchronous ground motion." Journal of Hydraulic Research 43(6): 710-719., Calayir and Karaton 2005aCalayir, Y. and M. Karaton (2005a). "A continuum damage concrete model for earthquake analysis of concrete gravity dam–reservoir systems." Soil Dynamics and Earthquake Engineering 25(11): 857-869., Calayir and Karaton 2005bCalayir, Y. and M. Karaton (2005b). "Seismic fracture analysis of concrete gravity dams including dam–reservoir interaction." Computers & Structures 83(19): 1595-1606., Akköse, Adanur et al. 2007Akköse, M., S. Adanur, A. Bayraktar and A. A. Dumanoğlu (2007). "Stochastic seismic response of Keban dam by the finite element method." Applied Mathematics and Computation 184(2): 704-714., Hacıefendioğlu, Başağa et al. 2007Hacıefendioğlu, K., H. B. Başağa, A. Bayraktar and Ş. Ateş (2007). "Nonlinear analysis of rock-fill dams to non-stationary excitation by the stochastic Wilson-θ method." Applied Mathematics and Computation 194(2): 333-345., Akköse, Adanur et al. 2008Akköse, M., S. Adanur, A. Bayraktar and A. A. Dumanoğlu (2008). "Elasto-plastic earthquake response of arch dams including fluid–structure interaction by the Lagrangian approach." Applied Mathematical Modelling 32(11): 2396-2412., Akkose, Bayraktar et al. 2008Akkose, M., A. Bayraktar and A. A. Dumanoglu (2008). "Reservoir water level effects on nonlinear dynamic response of arch dams." Journal of Fluids and Structures 24(3): 418-435., Bayraktar, Altunişik et al. 2009Bayraktar, A., A. C. Altunişik, B. Sevim, M. E. Kartal, T. Türker and Y. Bilici (2009). "Comparison of near- and far-fault ground motion effect on the nonlinear response of dam–reservoir–foundation systems." Nonlinear Dynamics 58(4): 655., Wang, Zhang et al. 2014Wang, G., S. Zhang, C. Wang and M. Yu (2014). "Seismic performance evaluation of dam-reservoir-foundation systems to near-fault ground motions." Natural hazards 72(2): 651-674.). In these solution methods, the dynamic response of the system is generally obtained by contribution of modes. The number of modes in FE modeling of DRFI and DRI problems is seen to be large in amount. Among this large amount of modes, the contribution of important modes on the dynamic responses of the problems is difficult to be estimated. Many tests and experimental analysis is needed to estimate the number of prominent modes. However, if the contribution of all modes is taken into account, it will take a long time to analyze the problem, which is inefficient and time-consuming. Another drawback with these solution technics is their stability, which needs further work. All of these drawbacks could be removed by solving the coupled equation of motion in the frequency domain using the direct method. In this method, the coupled equation of motion is transformed to the frequency domain. There is no need for further tests to find the number of important modes anymore. The radiation of reflected waves could easily be included by utilizing infinite elements in the model which is the most suitable approach in the FE method. As material damping, either hysteretic or Rayleigh damping models can be included in the coupled equation of motion. The whole system can be solved at the same time. This method is also faster than many other methods. But, the transformation of frequency domain responses to the time domain is the only further work in this method.

The dynamic analysis of roller compacted concrete gravity dam has been investigated by (Ghaedi, Hejazi et al. 2018Ghaedi, K., F. Hejazi, Z. Ibrahim and P. Khanzaei (2018). "Flexible Foundation Effect on Seismic Analysis of Roller Compacted Concrete (RCC) Dams Using Finite Element Method." KSCE Journal of Civil Engineering 22(4): 1275-1287.) considering DRFI effects and flexibility of foundation. It is reported that the flexibility of foundation domain has a great effect on the responses. The effect of base-rock characteristics on the dynamic response of DFI problem was investigated by (Bayraktar, Hançer et al. 2005bBayraktar, A., E. Hançer and M. Akköse (2005b). "Influence of base-rock characteristics on the stochastic dynamic response of dam–reservoir–foundation systems." Engineering Structures 27(10): 1498-1508.) considering rigid base foundation input, massless foundation input and deconvolved bedrock foundation models. In this paper, the responses obtained are compared with each other. By comparing the mean responses, it is reported that the massless-foundation input is inadequate model, but it is a practical technique. Similarly, the dynamic analysis of DRFI problem is investigated by (Ghaemian, Noorzad et al. 2019Ghaemian, M., A. Noorzad and H. Mohammadnezhad (2019). "Assessment of Foundation Mass and Earthquake Input Mechanism Effect on Dam–Reservoir–Foundation System Response." International Journal of Civil Engineering 17(4): 473-480., Salamon, Wood et al. 2021Salamon, J. W., C. Wood, M. A. Hariri-Ardebili, R. Malm and G. Faggiani (2021). Seismic Analysis of Pine Flat Concrete Dam: Formulation and Synthesis of Results, Cham, Springer International Publishing.) incorporating rigid, massless, and massed foundation models and it is reported that the massless foundation model overestimates the dynamic response of dam. The seismic responses of DRFI problem are evaluated by (Gorai and Maity 2019Gorai, S. and D. Maity (2019). "Seismic response of concrete gravity dams under near field and far field ground motions." Engineering Structures 196: 109292.) under near-fault and far-fault ground motion data. In the FE modeling of the problem, the width and depth of the foundation domain is assumed to be 3b and 1.5b, respectively, where b is the width of dam structure. The seismic performance of DRFI problem is investigated by (Chen, Yang et al. 2019Chen, D.-H., Z.-H. Yang, M. Wang and J.-H. Xie (2019). "Seismic performance and failure modes of the Jin'anqiao concrete gravity dam based on incremental dynamic analysis." Engineering Failure Analysis 100: 227-244.). In this FE model, an artificial viscous boundary condition is used at the truncated foundation boundaries with 3 times dam height in each direction. The seismic response evaluation of gravity dams including the effects of foundation mass is implemented by (Asghari, Taghipour et al. 2018Asghari, E., R. Taghipour, M. Bozorgnasab and M. Moosavi (2018). "Seismic Analysis of Concrete Gravity Dams Considering Foundation Mass Effect." KSCE Journal of Civil Engineering 22(12): 4988-4996.). In this FE analysis, viscous dampers are utilized at the truncated boundaries of the foundation. The dimension of foundation in the model is varied from 2-10H, where H is the height of dam. In this paper the length of reservoir domain is assumed to be constant which is 2H in all conditions. And it is reported that the distance 5H is relatively appropriate dimension for truncation of boundaries for foundation domain. The effect of reservoir length on the seismic performance of concrete gravity dams under near and far-fault ground motion is investigated by (Bayraktar, Türker et al. 2010Bayraktar, A., T. Türker, M. Akköse and Ş. Ateş (2010). "The effect of reservoir length on seismic performance of gravity dams to near- and far-fault ground motions." Natural Hazards 52(2): 257-275.). The depth of foundation (H) in the model is assumed to be constant. It is concluded that the length of 3 and 4H is the sufficient length for seismic performance evaluation of dams, where H is the height of dam structure. (Kartal, Cavusli et al. 2017Kartal, M. E., M. Cavusli and A. B. Sunbul (2017). "Assessing seismic response of a 2D roller-compacted concrete dam under variable reservoir lengths." Arabian Journal of Geosciences 10(22): 488.) demonstrated that the 3H is an appropriate length for seismic response evaluation of concrete gravity dams considering concrete nonlinearity. The same result is also reported by (Hariri-Ardebili, Seyed-Kolbadi et al. 2016Hariri-Ardebili, M. A., S. M. Seyed-Kolbadi and M. R. Kianoush (2016). "FEM-based parametric analysis of a typical gravity dam considering input excitation mechanism." Soil Dynamics and Earthquake Engineering 84: 22-43., Wang, Zhang et al. 2021Wang, C., H. Zhang, Y. Zhang, L. Guo, Y. Wang and T. T. Thira Htun (2021). "Influences on the Seismic Response of a Gravity Dam with Different Foundation and Reservoir Modeling Assumptions." Water 13(21).), while the length of 5H is reported by (Pelecanos, Kontoe et al. 2013Pelecanos, L., S. Kontoe and L. Zdravković (2013). "Numerical modelling of hydrodynamic pressures on dams." Computers and Geotechnics 53: 68-82.). While considering all these research papers, the researchers have investigated the most prominent topics regarding the dynamic behavior of gravity dams using FE method. Many different numerical modelings for dynamic analysis of this problem are proposed either in time or frequency domain. So far, very little attention has been paid to the impact of near-field domains size on the dynamic responses. There is no clear and consistent result reported on the precise dimension of near-field domains either for reservoir or foundation domain and both.

This study presents a novel and effective FE model for the dynamic analysis of DRI, DFI, and DRFI. The main goal of this study is determination of a consistent size regarding the near-field domains.

2 DESCRIPTION OF THE FE MODEL

The finite-infinite element model of DRI, DFI, and DRFI problems developed in this study is explained here. In this numerical model, the dam structure and near-field rock foundation domains were discretized by 8-noded iso-parametric quadrilateral finite elements. The reservoir domain was divided into near-field and far-field domains. The near-field reservoir domain includes the sloshing behavior of water with its compressibility. The near-field reservoir domain was discretized by 9-noded Lagrangian fluid elements. The free-surface region of water simulating the sloshing behavior was discretized by 1D 3-noded linear elements. These finite element types are well-known in the literature, and therefore, the basic mathematical formulas that lead to the interpolation functions of the elements are not given. The interpolation functions of the finite elements were taken from (Klaus-Jurgen 1982Klaus-Jurgen, B. (1982). Finite element procedures in engineering analysis, Prentice-Hall.).

The far-field reservoir domain which absorbs the reflected hydrodynamic waves was discretized by 7-noded fluid infinite elements. This infinite element model was taken from (Yerli, Temel et al. 1998Yerli, H., B. Temel and E. Kiral (1998). "Transient infinite elements for 2D soil-structure interaction analysis." Journal of geotechnical and geoenvironmental engineering 124(10): 976-988.) and then, it was modified to use in the unbounded reservoir domain. The foundation of concrete gravity dam and reservoir domain was assumed as bedrock. The far-field foundation domain which radiates the reflected waves was discretized by 7-noded solid infinite elements developed by (Yerli, Temel et al. 1998Yerli, H., B. Temel and E. Kiral (1998). "Transient infinite elements for 2D soil-structure interaction analysis." Journal of geotechnical and geoenvironmental engineering 124(10): 976-988.). The equations of motion for each domain (dam, reservoir and foundation) were obtained by the energy method. Then, the coupled equation of motion for the system was derived using frictionless contact elements at the interfaces between mediums (Akkas, Akay et al. 1979Akkas, N., H. Akay and C. Yilmaz (1979). "Applicability of general-purpose finite element programs in solid-fluid interaction problems." Computers and Structures 10(5): 773-783., Hacıefendioğlu 2009Hacıefendioğlu, K. (2009). "Stochastic response of concrete faced rockfill dams including partially ice-covered reservoir–foundation interaction under spatially varying seismic waves." Cold regions science and technology 58(1-2): 57-67., Hacıefendioğlu, Soyluk et al. 2012Hacıefendioğlu, K., K. Soyluk and F. Birinci (2012). "Numerical investigation of stochastic response of an elevated water tank to random underground blast loading." Stochastic environmental research and risk assessment 26(4): 599-607.). The consistent mass matrices of the elements in the dam, near-field reservoir domain, and near-field foundation domain was calculated by 3x3 Gauss numerical integration method. The stiffness matrices of the elements in the dam and foundation domain were calculated by 3x3 Gauss numerical integration method, while the stiffness matrix of the elements in the near-field reservoir domain and the free surface of the reservoir was calculated by 2x2 reduced Gauss numerical integration method as recommended by (Wilson and Khalvati 1983Wilson, E. L. and M. Khalvati (1983). "Finite elements for the dynamic analysis of fluid‐solid systems." International Journal for Numerical Methods in Engineering 19(11): 1657-1668.). The mass and stiffness matrices of the infinite elements in the far-end boundaries of the reservoir and foundation domain were calculated using the Newton-Cotes numerical integration method (Bettess and Zienkiewicz 1977Bettess, P. and O. Zienkiewicz (1977). "Diffraction and refraction of surface waves using finite and infinite elements." International Journal for Numerical Methods in Engineering 11(8): 1271-1290.).

The coupled equation of motion was transformed into the Laplace domain and then solved by the direct method. In this method, the damping mechanism of the system can be assumed as either frequency-independent (hysteretic damping) or frequency dependent (Rayleigh damping) parameter. In this paper, the damping parameter was assumed as hysteretic. The coupled system was solved at the same time. The complex responses computed in the Laplace domain were finally transformed into time-domain responses by using Durbin's method (Durbin 1974Durbin, F. (1974). "Numerical inversion of Laplace transforms: an efficient improvement to Dubner and Abate's method." The Computer Journal 17(4): 371-376.).

2.1 Finite Element Formulation of the Problem Using Lagrangian Approach

The coupled equation of motion for the DRFI problem was derived using the Lagrangian approach. In this approach, the fluid elements were assumed as linearly elastic, nonrotating, and non-viscous. The basic equations for the fluid domain in the Lagrangian approach were given here according to (Akköse and Şimşek 2010Akköse, M. and E. Şimşek (2010). "Non-linear seismic response of concrete gravity dams to near-fault ground motions including dam-water-sediment-foundation interaction." Applied Mathematical Modelling 34(11): 3685-3700.). In the dynamic motion of the dam, foundation and the reservoir domain, the displacements are defined as unknown parameters. Therefore, the compatibility condition at the dam-reservoir and reservoir foundation interfaces was automatically satisfied. The stress-strain relationships of a two-dimensional fluid element in matrix form can be given as follows:

where P, ${\text{C}}_{\text{11}}$ and ${\text{ε}}_{\text{v}}$ represents the pressure, bulk modulus and volumetric strain (${\text{ε}}_{\text{v}}\text{=}\frac{\partial \text{u}}{\partial \text{x}}\text{+}\frac{\partial \text{v}}{\partial \text{y}}$) of the fluid element, respectively. Since the fluid element was assumed to be irrotational, a constraint parameter can be included in the stress-strain Eq. (1) using penalty methods (Klaus-Jurgen 1982Klaus-Jurgen, B. (1982). Finite element procedures in engineering analysis, Prentice-Hall., Zienkiewicz, Taylor et al. 2000Zienkiewicz, O. C., R. L. Taylor, R. L. Taylor and R. L. Taylor (2000). The finite element method: solid mechanics, Butterworth-heinemann.). In this equation,${\text{P}}_{\text{z}}$, ${\text{C}}_{\text{22}}$and ${\text{W}}_{\text{z}}$ denotes the rotational stress, constraint parameter and the rotation of the fluid element around z axis (${\text{W}}_{\text{z}}\text{=}\frac{\text{1}}{\text{2}}\left(\frac{\partial \text{u}}{\partial \text{y}}\text{-}\frac{\partial \text{v}}{\partial \text{x}}\right)$) normal to the (x-y) plane, respectively. The equation of motion for the fluid domain was derived using the energy principle; the finite element formulation of the total strain energy for the fluid domain can be written as follows:

where ${\underset{¯}{\text{U}}}_{\text{f}}$and $K_{\text{f}}$are nodal displacement vector and stiffness matrix of fluid domain, respectively. An important behavior of water is the capability of displacements without altering its volume. This behavior can be observed at the surface of reservoirs and storage tanks in the form of sloshing waves that move up and down. The effects of sloshing waves at the surface, which are an important aspect, can be considered in Lagrangian approach. It is feasible to consider the effect of free surface waves in the terms of fluid potential energy. The potential energy of fluid system can be written as:

where ${\underset{¯}{\text{U}}}_{\text{sf}}$ and $S_{\text{f}}$ are the up-down nodal displacement vector and stiffness matrix of surface fluid elements, respectively. Moreover, the kinetic energy of the fluid domain can be written as:

where ${\underset{¯}{\stackrel{.}{\text{U}}}}_{\text{f}}$ and $M_{\text{f}}$ are the nodal velocity vector and mass matrix of fluid domain, respectively. The following set of equations that describe the motion of fluid domain can be obtained substituting Eqs. (2-4), into the Lagrange’s (Clough and Penzien 1975Clough, R. and J. Penzien (1975). Dynamics of Structures, McGraw-Hill.) equation as:

where ${\underset{¯}{\stackrel{\mathrm{..}}{\text{U}}}}_{\text{f}}$and ${\widehat{K}}_{\text{f}}$ are the nodal acceleration vector and the stiffness matrix of fluid domain including the free surface stiffness matrix, respectively. The stiffness matrix of the entire system can be calculated by summing the stiffness matrices of all 2D elements ($K_{\text{f}}={\displaystyle \sum K_{\text{f}}^{\text{e}}}$) with the stiffness matrices of all 1D surface elements ($S_{\text{f}}={\displaystyle \sum S_{\text{f}}^{\text{e}}}$). Furthermore, the mass matrix of the whole system calculates by the summation of the mass matrices of all 2D elements ($M_{\text{f}}={\displaystyle \sum M_{\text{f}}^{\text{e}}}$) exist in fluid domain. The stiffness and mass matrices of one element can be calculated as:

where, $B_{\text{f}}^{\text{e}}$ and $C_{\text{f}}$ in Eq. (6) denotes the strain-displacement matrix and elasticity matrix (Eq. (1)) of a fluid element, respectively. Furthermore, ${\text{ρ}}_{\text{f}}\text{, g}$, $h_{\text{f}}^{\text{e}}$ and $H_{\text{f}}^{\text{e}}$ expresses the mass density of water, gravitational constant and a vector containing interpolation functions of 1D surface element and 2D finite fluid element, respectively. The same procedure can be adopted deriving the FE equation of motion for dam structure. The coupled equation of motion for dam-reservoir interaction problem may be obtained by imposing the boundary conditions at the interface between reservoir and dam structure. For this purpose, the nodal displacements perpendicular to the interface between two mediums must be equal:

This boundary condition can be implemented using truss elements (Akkas, Akay et al. 1979Akkas, N., H. Akay and C. Yilmaz (1979). "Applicability of general-purpose finite element programs in solid-fluid interaction problems." Computers and Structures 10(5): 773-783., Hacıefendioğlu 2009Hacıefendioğlu, K. (2009). "Stochastic response of concrete faced rockfill dams including partially ice-covered reservoir–foundation interaction under spatially varying seismic waves." Cold regions science and technology 58(1-2): 57-67., Hacıefendioğlu, Soyluk et al. 2012Hacıefendioğlu, K., K. Soyluk and F. Birinci (2012). "Numerical investigation of stochastic response of an elevated water tank to random underground blast loading." Stochastic environmental research and risk assessment 26(4): 599-607.) between interface nodes or penalty methods (Klaus-Jurgen 1982Klaus-Jurgen, B. (1982). Finite element procedures in engineering analysis, Prentice-Hall.). In this paper, short and axially almost rigid truss elements were used as interface link elements. Therefore, the coupled equation of motion for dam-structure interaction problem subjected to ground acceleration can be written as follows:

where, M, K, and ${\underset{¯}{a}}_g\text{(t)}$ represent the mass matrix, stiffness matrix of the coupled system and time varying ground acceleration vector, respectively. $\underset{¯}{\stackrel{\mathrm{..}}{\text{U}}}$and $\underset{¯}{\text{U}}$express the nodal acceleration vector and nodal displacement vector of the coupled system, respectively. By including the hysteretic damping parameter into the coupled Laplace transformed Eq. (8), a set of linear algebraic equations can be obtained as:

where, s,$\underset{¯}{\stackrel{⏜}{\text{U}}}$, and ${\underset{¯}{\stackrel{⏜}{a}}}_g\text{(s)}$represent complex Laplace variable, Laplace transformed displacement vector and Laplace transformed ground acceleration vector, respectively. Also, $\text{η}$ and i denote the hysteretic damping factor and imaginary number ($\text{i=}\sqrt{\text{-1}}$). The hysteretic damping factor in this equation can be written as $\text{η=2ζ}$ (Gutierrez and Chopra 1978Gutierrez, J. A. and A. K. Chopra (1978). "A substructure method for earthquake analysis of structures including structure‐soil interaction." Earthquake Engineering and Structural Dynamics 6(1): 51-69.). The solution of Eq. (9) at each time step can easily be implemented using Cholesky method.

2.2 Mathematical Formulations of the Infinite Element Model

Only the necessary interpolation functions for infinite elements are given here. The mapping functions of 7-noded infinite element are given as:

The horizontal and vertical displacement fields of 7-noded infinite element in the Laplace domain can be calculated as:

where ${\widehat{\text{N}}}_{\text{i}}$ denotes the displacement shape functions of the infinite element, and can be written as:

In Eq. (12), ${\widehat{\text{P}}}_{\text{k}}\left(\text{ξ}\text{, s}\right)$ expresses the wave-propagation function in the infinite element. This type of function is originally developed by (Chuhan and Chongbin 1987Chuhan, Z. and Z. Chongbin (1987). "Coupling method of finite and infinite elements for strip foundation wave problems." Earthquake engineering and structural dynamics 15(7): 839-851.) in the Laplace domain such as,

The propagation of compression (P), shear (S) and Rayleigh (R) waves in the mediums are included in the infinite element model developed by (Yerli, Temel et al. 1998Yerli, H., B. Temel and E. Kiral (1998). "Transient infinite elements for 2D soil-structure interaction analysis." Journal of geotechnical and geoenvironmental engineering 124(10): 976-988.). Thus, the propagation function is written as,

where α and c_k are the displacement decay parameter and corresponding wave velocity, respectively. The $\text{ξ}$ and L parameters in the function show the infinite direction and characteristic length of the element. Furthermore, the displacement behavior in many unbounded domains could be the form of 1/x function. By selecting an adequate decay coefficient α, it is possible to develop an exponential propagation function as (Bettess and Zienkiewicz 1977Bettess, P. and O. Zienkiewicz (1977). "Diffraction and refraction of surface waves using finite and infinite elements." International Journal for Numerical Methods in Engineering 11(8): 1271-1290.),

where X₁ and X₄ are the x coordinate of the element at the point 1 and 4 (Figure 1). The constants a, b, and c given in the function can be calculated by using the Eq. (11). For this purpose, the following relationships were obtained by considering the nodes 1, 4, and 6 (or 3, 5, and 7) in the infinite element:

Thus, the constants a, b and c in this equation can easily be obtained by assuming ${\widehat{\text{P}}}_{\text{1}}=\mathrm{1,}{\widehat{\text{P}}}_{\text{2}}=\mathrm{0,}{\widehat{\text{P}}}_{\text{3}}=0$,${\widehat{\text{P}}}_{\text{1}}=\mathrm{0,}{\widehat{\text{P}}}_{\text{2}}=\mathrm{1,}{\widehat{\text{P}}}_{\text{3}}=0$, ${\widehat{\text{P}}}_{\text{1}}=\mathrm{0,}{\widehat{\text{P}}}_{\text{2}}=\mathrm{0,}{\widehat{\text{P}}}_{\text{3}}=1$ and their satisfaction with Eq. (11) in three different cases (Yerli, Temel et al. 1998Yerli, H., B. Temel and E. Kiral (1998). "Transient infinite elements for 2D soil-structure interaction analysis." Journal of geotechnical and geoenvironmental engineering 124(10): 976-988.).

The infinite element stiffness and mass matrices were calculated numerically. The typical Gauss-Legendre integration technique was applied for the finite direction. However, Newton-Cotes numerical integration scheme was used to calculate infinite integrals. The typical expression of Newton-Cotes numerical integration scheme is shown in Eq. (17):

where $\text{F(ζ)}$, ${\text{β}}_{\text{i}}$, W_k and ξ_k represent the displacement interpolation function, the type of propagated waves at the infinite direction weight values and abscissa values, respectively (for details: Bettess and Zienkiewicz, 1977Bettess, P. and O. Zienkiewicz (1977). "Diffraction and refraction of surface waves using finite and infinite elements." International Journal for Numerical Methods in Engineering 11(8): 1271-1290.).

3 VALIDATIONS OF THE DEVELOPED PROGRAMS

Many examples from previous studies were selected, and comparative studies were conducted to verify the results of the developed programs in FORTRAN 90. Generally, three types of analyzes were carried out in these investigations. Modal analysis, static, and dynamic analysis were implemented. In the modal analysis section, the natural frequencies of rigid fluid tank (Akköse and Şimşek 2010Akköse, M. and E. Şimşek (2010). "Non-linear seismic response of concrete gravity dams to near-fault ground motions including dam-water-sediment-foundation interaction." Applied Mathematical Modelling 34(11): 3685-3700.), Koyna dam (Chopra and Chakrabarti 1971Chopra, A. and P. Chakrabarti (1971). "The Koyna Earthquake of December 11, 1967 and the Performance of Koyna Dam, Report No. EERC 71-1." Earthquake Engineering Research Center, University of California, Berkeley., Nayak and Maity 2013Nayak, P. and D. Maity (2013). "Seismic Damage Analysis of Aged Concrete Gravity Dams." International Journal for Computational Methods in Engineering Science and Mechanics 14(5): 424-439.), Sariyar dam-reservoir interaction (Köseoğlu 2007Köseoğlu, E. (2007). Sıvı-yapı etkileşim sistemlerinin Lagrange yaklaşımıyla stokastik dinamik analizi, Karadeniz Teknik Üniversitesi/Fen Bilimleri Enstitüsü.), Sariyar dam-reservoir-foundation interaction (Köseoğlu 2007Köseoğlu, E. (2007). Sıvı-yapı etkileşim sistemlerinin Lagrange yaklaşımıyla stokastik dinamik analizi, Karadeniz Teknik Üniversitesi/Fen Bilimleri Enstitüsü.) problems were analyzed. The static analysis of Sariyar dam-foundation interaction (Bayraktar 1991Bayraktar, A. (1991). Beton ağırlık barajlarda baraj-su-zemin etkileşiminin statik ve dinamik analizde değerlendirilmesi, Karadeniz Teknik Üniversitesi/Fen Bilimleri Enstitüsü.) problem was carried out in the static analysis part. Finally, the dynamic analysis of Koyna dam, Koyna dam-reservoir interaction (Chopra and Chakrabarti 1971Chopra, A. and P. Chakrabarti (1971). "The Koyna Earthquake of December 11, 1967 and the Performance of Koyna Dam, Report No. EERC 71-1." Earthquake Engineering Research Center, University of California, Berkeley., Chopra and Chakrabarti 1973Chopra, A. K. and P. Chakrabarti (1973). "The Koyna earthquake and the damage to Koyna dam." Bulletin of the Seismological Society of America 63(2): 381-397.) and soil-structure interaction (Düzgün 2007Düzgün, O. (2007). Effects of topography on the dynamic response of soil structure systems (in Turkish), PhD thesis, Atatürk University, Erzurum, Turkey.) problems were carried out in the dynamic analysis part. In this paper, only the comparative results for natural frequencies of a rigid fluid tank (Akköse and Şimşek 2010Akköse, M. and E. Şimşek (2010). "Non-linear seismic response of concrete gravity dams to near-fault ground motions including dam-water-sediment-foundation interaction." Applied Mathematical Modelling 34(11): 3685-3700.), the dynamic analysis of Koyna dam (Chopra and Chakrabarti 1971Chopra, A. and P. Chakrabarti (1971). "The Koyna Earthquake of December 11, 1967 and the Performance of Koyna Dam, Report No. EERC 71-1." Earthquake Engineering Research Center, University of California, Berkeley., Chopra and Chakrabarti 1973Chopra, A. K. and P. Chakrabarti (1973). "The Koyna earthquake and the damage to Koyna dam." Bulletin of the Seismological Society of America 63(2): 381-397.) and soil-sturucture interaction (Düzgün 2007Düzgün, O. (2007). Effects of topography on the dynamic response of soil structure systems (in Turkish), PhD thesis, Atatürk University, Erzurum, Turkey.) problems were given. The coupled equation of motion for Koyna dam problem was solved by mode superposition method, while the equation of motion for soil-structure interaction problem was solved by direct method in the Laplace domain.

The finite element model of rigid fluid tank is shown in Figure 2. Several studies (Lamb 1924Lamb, H. (1924). Hydrodynamics, University Press., Olson and Bathe 1983Olson, L. G. and K.-J. Bathe (1983). "A study of displacement-based fluid finite elements for calculating frequencies of fluid and fluid-structure systems." Nuclear Engineering and Design 76(2): 137-151., Wilson and Khalvati 1983Wilson, E. L. and M. Khalvati (1983). "Finite elements for the dynamic analysis of fluid‐solid systems." International Journal for Numerical Methods in Engineering 19(11): 1657-1668., Akköse and Şimşek 2010Akköse, M. and E. Şimşek (2010). "Non-linear seismic response of concrete gravity dams to near-fault ground motions including dam-water-sediment-foundation interaction." Applied Mathematical Modelling 34(11): 3685-3700.) have examined the modal analysis of this tank. The first 10 natural frequencies of the tank were calculated and compared in Table 1. According to Table 1, the first 5 modes reflect the sloshing modal frequencies of the fluid tank, which are in the range of 0.3567-0.7653 Hz. The first sloshing frequency of the tank given in Table 1 was analytically calculated by (Lamb 1924Lamb, H. (1924). Hydrodynamics, University Press., Wilson and Khalvati 1983Wilson, E. L. and M. Khalvati (1983). "Finite elements for the dynamic analysis of fluid‐solid systems." International Journal for Numerical Methods in Engineering 19(11): 1657-1668.). Also, the sixth modal frequency given in Table 1 was the first compressive modal frequency of the rigid tank analytically calculated by (Olson and Bathe 1983Olson, L. G. and K.-J. Bathe (1983). "A study of displacement-based fluid finite elements for calculating frequencies of fluid and fluid-structure systems." Nuclear Engineering and Design 76(2): 137-151.). It is clearly seen that the results obtained from the developed program have an excellent agreement with analytical and numerıcal results given in (Lamb 1924Lamb, H. (1924). Hydrodynamics, University Press., Olson and Bathe 1983Olson, L. G. and K.-J. Bathe (1983). "A study of displacement-based fluid finite elements for calculating frequencies of fluid and fluid-structure systems." Nuclear Engineering and Design 76(2): 137-151., Wilson and Khalvati 1983Wilson, E. L. and M. Khalvati (1983). "Finite elements for the dynamic analysis of fluid‐solid systems." International Journal for Numerical Methods in Engineering 19(11): 1657-1668., Akköse and Şimşek 2010Akköse, M. and E. Şimşek (2010). "Non-linear seismic response of concrete gravity dams to near-fault ground motions including dam-water-sediment-foundation interaction." Applied Mathematical Modelling 34(11): 3685-3700.).

Furthermore, the dynamic analysis of the Koyna gravity dam implemented in (Chopra and Chakrabarti 1971Chopra, A. and P. Chakrabarti (1971). "The Koyna Earthquake of December 11, 1967 and the Performance of Koyna Dam, Report No. EERC 71-1." Earthquake Engineering Research Center, University of California, Berkeley., Chopra and Chakrabarti 1973Chopra, A. K. and P. Chakrabarti (1973). "The Koyna earthquake and the damage to Koyna dam." Bulletin of the Seismological Society of America 63(2): 381-397.) is repeated. The problem is solved under the horizontal component of the 1967 Koyna ground motion (Figure 3). Detailed information about this problem can be obtained from (Chopra and Chakrabarti 1971Chopra, A. and P. Chakrabarti (1971). "The Koyna Earthquake of December 11, 1967 and the Performance of Koyna Dam, Report No. EERC 71-1." Earthquake Engineering Research Center, University of California, Berkeley., Chopra and Chakrabarti 1973Chopra, A. K. and P. Chakrabarti (1973). "The Koyna earthquake and the damage to Koyna dam." Bulletin of the Seismological Society of America 63(2): 381-397.). As a result of the analysis, the horizontal displacement time-histories of the crest point at the dam are comparatively given in Figure 4. The discrepancies between the maximum displacement responses given in (Chopra and Chakrabarti 1971Chopra, A. and P. Chakrabarti (1971). "The Koyna Earthquake of December 11, 1967 and the Performance of Koyna Dam, Report No. EERC 71-1." Earthquake Engineering Research Center, University of California, Berkeley., Chopra and Chakrabarti 1973Chopra, A. K. and P. Chakrabarti (1973). "The Koyna earthquake and the damage to Koyna dam." Bulletin of the Seismological Society of America 63(2): 381-397.) and calculated results are seen to be 1.32% for dam and 0.07% for dam-reservoir problem. The dynamic results calculated by the developed program appear to be in good agreement with the results given in (Chopra and Chakrabarti 1971Chopra, A. and P. Chakrabarti (1971). "The Koyna Earthquake of December 11, 1967 and the Performance of Koyna Dam, Report No. EERC 71-1." Earthquake Engineering Research Center, University of California, Berkeley., Chopra and Chakrabarti 1973Chopra, A. K. and P. Chakrabarti (1973). "The Koyna earthquake and the damage to Koyna dam." Bulletin of the Seismological Society of America 63(2): 381-397.).

Finally, the dynamic analysis of soil-structure interaction problem implemented in (Düzgün 2007Düzgün, O. (2007). Effects of topography on the dynamic response of soil structure systems (in Turkish), PhD thesis, Atatürk University, Erzurum, Turkey.) is reanalyzed here in order to control the results of direct method and the correct work of the 7-noded infinite element model in unbounded soil medium. The detailed information about this problem can be obtained from (Düzgün 2007Düzgün, O. (2007). Effects of topography on the dynamic response of soil structure systems (in Turkish), PhD thesis, Atatürk University, Erzurum, Turkey.). The comparative time-history displacement response of the point selected at the top of the structure is given in Figure 5. The discrepancy between the maximum displacement responses is seen to be about 5%.

By considering the comparative results given in Table 1, Figure 4 and Figure 5 it can be concluded that the developed programs in FORTRAN 90 give quite accurate results. It has also been seen that the coupling of finite and infinite elements works accurately to absorb the reflected dynamic waves in the unbounded mediums.

4 PARAMETRIC STUDIES

The dynamic responses of three different problems were examined here. Furthermore, the effect of near-field domain size on the dynamic response of DRI, DFI, and DRFI problems was evaluated by including and not including the unbounded domains in the model. The dynamic responses of these problems were analyzed under two different conditions. The boundaries of the near-field domain in rock foundation were constrained in x and y directions, only the x direction was constrained in the reservoir domain with the name WithOut Infinite Domain (WOID); in the second condition, the infinite elements were used in the far-field boundaries of the model. This condition was named With Infinite Domain (WID). These conditions are further explained in Figures 6-7.

The length and depth of the near-field foundation domain as well as the length of the near-field reservoir domain were assumed to be from 1-15 times dam height. The dimensions of domains in the problems were given in Figures 6-7, and Table 2. The suitable finite element mesh of the coupled problems was obtained by calculating the natural frequencies of the system preparing many different mesh schemes. The material properties of dam, water and rock foundation are given in Table 2. The length and stiffness of interaction elements used at the interfaces between water-dam and water-bedrock were assumed to be 0.001 m and 2x10¹⁴, respectively. This assumption correctly satisfies the dynamic interaction between these domains (Akkas, Akay et al. 1979Akkas, N., H. Akay and C. Yilmaz (1979). "Applicability of general-purpose finite element programs in solid-fluid interaction problems." Computers and Structures 10(5): 773-783.).

The hysteretic damping factor for the coupled system was assumed as $\text{η=10\%}$. This damping factor is equivalent to 5% viscous damping, which was experimentally measured in several dams (Chopra 2020Chopra, A. K. (2020). Earthquake Engineering for Concrete Dams: Analysis, Design, and Evaluation, John Wiley and Sons.). The displacement decay parameter used in the infinite elements was selected as 4.08 (Bettess and Bettess 1984Bettess, P. and J. A. Bettess (1984). "Infinite elements for static problems." Engineering Computations.). The horizontal acceleration time-histories of the Koyna earthquake were imposed to the models as dynamic loads (Figure 3). In the Laplace domain analysis, the number of Laplace parameters (m) is assumed to be 10. The α parameter in the Laplace transformation is assumed to be αT = 6 (Durbin 1974Durbin, F. (1974). "Numerical inversion of Laplace transforms: an efficient improvement to Dubner and Abate's method." The Computer Journal 17(4): 371-376.) where T is the time duration and equal 10.24 sec. The time step in the analysis was selected as 0.01 sec. In each dynamic analysis, there were a total of 1,024 steps. The constraint parameter in the water was assumed to be 100 times its bulk modulus (Wilson and Khalvati 1983Wilson, E. L. and M. Khalvati (1983). "Finite elements for the dynamic analysis of fluid‐solid systems." International Journal for Numerical Methods in Engineering 19(11): 1657-1668.). The velocity of compressive (P), shear (S) and Rayleigh (R) waves in the rock foundation was assumed to be 2,300, 950 and 450 m/sec, respectively. The velocity of compressive waves in the water was assumed to be 2,100 m/sec (Schumann, Stipp et al. 2014Schumann, K., M. Stipp, J. H. Behrmann, D. Klaeschen and D. Schulte‐Kortnack (2014). "P and S wave velocity measurements of water‐rich sediments from the Nankai Trough, Japan." Journal of Geophysical Research: Solid Earth 119(2): 787-805.). The flow charts for expressing the research methodology and demonstration of the steps in the developed programs are given in Figures 8-9, respectively.

5 RESULTS AND DISCUSSIONS

In this section, the applicability of the numerical model is assessed with dynamic responses of DRI, DFI, and DRFI problems under horizontal component of 1967 Koyna earthquake. The dynamic responses of DRI, DFI, and DRFI problems are evaluated in individual subsections. In the last part, the results of these problems are compared and discussed.

5.1 Dynamic Responses of DRI Problem

In the dynamic analysis of the DRI problem, the length of the near-field reservoir domain is varied from 1-7H, where H is the height of the dam. As mentioned before, the boundary condition at the far-field domain was modeled in two conditions; (i) constrained in the x-direction (WOID) and (ii) infinite fluid elements (WID). The foundation was assumed to be rigid (Figures 6a, 7a). The dynamic response obtained from condition (i) is expected to be inaccurate in contrast with the response of condition (ii) due to the reflection of hydrodynamic waves in the bounded reservoir domain. Condition (ii) is more realistic because the transmission of reflected hydrodynamic waves to infinity is accurately provided by using the proposed infinite fluid elements.

The distribution of maximum hydrodynamic pressures at the upstream face of the dam model for different reservoir lengths is compared in Figure 10. According to Figure 10, it can be seen that the hydrodynamic pressure at the upstream face of the dam is varied by variation of the reservoir length. In general, it is seen to be augmented by increasing the reservoir length. The distribution of hydrodynamic pressure at the dam face given in Figure 10a is significantly higher than that response given in Figure 10b. The behavior of the response is seen to be not reasonable due to reflection of hydrodynamic waves in the bounded reservoir domain. In Figure 10b, the response is seen to be more reasonable and the difference between each length response is smaller than those responses given in Figure 10a. The results also reveal that the absorption capacity of the far-field reservoir domain significantly reduces the hydrodynamic pressures at the heel of the dam. Furthermore, the corresponding hydrodynamic pressure approaches infinity at the resonant frequency when the far-field reservoir boundary is constrained, absorbing energy in the far-field reservoir boundary always results in the smallest hydrodynamic response of the dam. For instance, the peak value of hydrodynamic pressure at the heel of the dam for (WOID, α = 1) and (WID, α = 1) conditions are seen to be 0.675 and 0.246 MPa, respectively. The decrease in the responses was founded at about 63%. It can be concluded that the hydrodynamic pressure will be overestimated by constraining the far-field reservoir boundaries (without infinite elements) and the hydrodynamic pressure will be underestimated by including the absorbing far-field reservoir boundary (with infinite fluid elements) in the model.

The horizontal crest displacement time histories of dam structure for reservoir lengths 1H, 4H, and 7H are given in Figure 11. According to Figure 11, it is seen that the displacement time histories are different in each reservoir length. The displacement response for (WOID, α = 1) condition is seen to be the maximum, while in (WID, α = 1) condition it is the minimum response over 1H, 4H, and 7H lengths. Similarly, the absolute maximum displacement responses at the upstream face of the dam are given in Figure 12. According to this figure, it is seen that by extending the length of the near-field reservoir, the crest displacements are also increased. For example, an increase in the near-field reservoir length for (WID, α = 1) and (WID, α = 3) conditions leads to an increase in the absolute maximum displacements from 0.0162 to 0.0178 m with a difference about 9%. With increasing the length of the near-field reservoir domain, the hydrodynamic pressure increases, and thus, the crest displacements are amplified. The maximum compressive and tensile stress responses at the base of the dam model are compared and given in Figure 13. According to this figure, the same behavior is also observed while investigating the stress responses. The stress responses at the base of the dam structure for (WID) condition are seen to be nearly the same for all reservoir lengths. Moreover, the stress responses are smaller than those responses given in (WOID) condition. Thus, it should be noted that the infinite fluid elements used in the far-field reservoir domain effectively absorb the reflected waves and they can represent reality.

5.2 Dynamic Responses of DFI Problem

The dynamic responses of the DFI problem was investigated in this section. The length and depth of rock foundation domain were varied from 1-15H in both upstream and downstream sides of the dam structure (Figures 6b, 7b). As mentioned in the previous section, two conditions (WOID) and (WID) were also considered in this problem. The displacement time histories of the crest point of the dam structure for both conditions with α = 1, 4, 7 are given in Figure 14. According to this figure, it is seen that the displacement responses are different for different α values. Moreover, by increasing the size of near-field domains in the problem, the crest displacement is also increased. While including the radiation of waves in the foundation domain, the displacement responses are significantly decreased (Figure 14b). By imposing constrained boundary conditions to the far-field foundation domain, the responses are increased due to the reflection of waves in the bounded rock foundation domain. The absolute maximum displacement responses at the upstream face of the dam structure are given in Figure 15. Only the responses for α =1-7 are given. According to Figure 15, it is seen that by increasing the near-field domain size, the displacement responses significantly increase in both conditions. But, the responses of (WID) condition for all α values are the smallest due to absorption of the reflected wave in the unbounded foundation domain. The maximum compressive and tensile stress responses at the base of the dam structure is given in Figure 16; in this figure, only the responses for α =1-7 in both conditions are given. It is also observed that by increasing the near-field domain size, the stress responses at the base of the dam significantly increase. This trend in (WID) condition is much slower and smaller than (WOID) condition.

5.3 Dynamic Responses of DRFI Problem

The dynamic responses of the DRFI problem are investigated in this section. The length of the near-field reservoir as well as the length and depth of the near-field rock foundation at upstream face and downstream face sides are varied from 1-15H (Figures 6c, 7c). Only the responses for α = 1-7 are given here. The hydrodynamic pressure at the dam face in (WOID) and (WID) conditions are given in Figure 17. According to this figure, it is seen that the hydrodynamic pressure at the dam face increases by increasing the dimension of near-field domains. The behavior of hydrodynamic response in (WOID) and (WID) conditions can be seen to be different due to different boundary conditions imposed on the far-field boundaries. The response in (WID) condition is seen to be significantly bigger than (WOID) condition due to the reflection of waves in the bounded domains. The hydrodynamic pressure response in (WOID) condition is far from reality because the transmission of reflected waves to unbounded domains (water and rock foundation) is not provided. It can be seen that in the DRI and DRFI problems, the hydrodynamic pressure increases due to including the effect of rock foundation in the analysis. For example, the peak of hydrodynamic pressure at the heel of the dam for DRI and DRFI problems in (WID, α = 3) condition is 0.311 and 0.972 MPa, respectively.

In other words, the peak of hydrodynamic pressure is increased by about 68% for the dam with the rock foundation domain effect. The displacement time-history response for the crest point of the dam structure is given in Figure 18. Only, the responses for α = 1, 4, and 7 are compared. It is observed that the responses are generally increased by increasing the size of near-field domains. Furthermore, the responses for (WID) condition are seen to be smaller than the responses for (WOID) condition. The same behavior is also observed in the absolute maximum displacement responses at the upstream face of the dam structure which is given in Figure 19.

The distribution of maximum compressive and tensile stresses at the base of the dam in (WID) and (WOID) conditions are given in Figure 20. According to this figure, the stress responses vary by variation of near-field domain size in the problem. Furthermore, the responses in (WOID) condition are significantly bigger than those responses obtained in (WID) condition.

The significant stress responses in the dam body generally occur at the heel, toe, and some regions at the below of the neck point. The stress concentration generally occurs in these regions of the dam body. Therefore, the variation of compressive and tensile stress responses at the heel of dam structure from α = 1-15 in (WOID) and (WID) conditions are comparatively investigated here. The stress responses of DRI, DFI, and DRFI problems are individually evaluated in Figures 21a, 21b, 21c, and finally the responses of these problems are compared with each other in Figure 21d.

According to Figure 21a, it is seen that the stress responses vary with the variation of near-field domain size in (WOID) and (WID) conditions. Moreover, the responses after α = 4.5 for both conditions are seen to be nearly the same. The maximum tensile stress for condition (WOID) is calculated as 3.159 MPa and for condition (WID), it is calculated as 3.138 MPa. The difference between responses is seen to be 0.66%. The same evaluation is also conducted on DFI (Figure 21b) and DRFI (Figure 21c) problems. And it is seen that the stress responses after α = 12 are nearly the same and the variations are negligible. The maximum compressive stress responses at the heel of the dam in the DFI problem for (WOID, α = 12) and (WID, α = 12) conditions are seen to be -5.7 and -5.409 MPa, respectively. The difference between responses is calculated as 5.1%. The same evaluation for the DRFI problem is conducted and the difference is found to be 0.24%.

The stress responses at the heel of dam structure for DRI, DFI, and DRFI problems are comparatively evaluated in Figure 21d. According to this figure, the stress response of the DRI problem is significantly smaller than the DRFI problem. Thus, the assumption of a rigid foundation in the DRI problem gives inaccurate response. To further clarify this condition, the compressive stress responses at the heel of the dam in the DRI and DRFI problems for (WID, α = 4.5) condition are seen to be -2.953 and -38.563 MPa, respectively. The difference between responses of the DRFI and DRI problems is calculated as 92.34%. The effect of reservoir on the stress responses at the heel of the dam is also important. To demonstrate the importance, the stress response in the DRFI and DFI problems is evaluated considering (WID, α = 4.5) condition. The tensile stress responses at the heel of the dam in the DRFI problem can be seen to be 44.452 MPa and it is 21.966 MPa in the DFI problem. The rise of stress at that point is seen to be 50.5%. According to Figure 21d, it can be seen that the radiation effect of waves in the bounded domains for the DFI and DRFI problems can effectively be provided if the near-field domains size would be assumed as 12H. By taking the size of near-field domains as 12 times dam height, there is no need for including the far-field domains in the model. Because the near-field domain size is sufficiently big to absorb the reflected dynamic waves. And thus, the truncated boundaries can be constrained in the x and y direction for foundation and only the x-direction for water. At the truncated boundaries, there is no need for infinite elements or transmitting boundaries anymore. The stress responses in (WOID, α = 12) and (WID, α = 12) conditions are seen to be approximately the same. The truncation of boundaries at this size is not economic because the analysis process takes a long time and it is time-consuming. Furthermore, if the stress level in (WOID, α = 12) or (WID, α = 12) condition is the correct response, the same stress level can be obtained by assuming the near-field domain size of 1.5H. This condition is schematically explained in Figure 22 for tensile stress variation at the heel of the dam structure. The same result can also be obtained while observing the compressive stress response.

The near-filed domain size in the DRFI problem is concluded to be three times (Chen, Yang et al. 2019Chen, D.-H., Z.-H. Yang, M. Wang and J.-H. Xie (2019). "Seismic performance and failure modes of the Jin'anqiao concrete gravity dam based on incremental dynamic analysis." Engineering Failure Analysis 100: 227-244., Gorai and Maity 2019Gorai, S. and D. Maity (2019). "Seismic response of concrete gravity dams under near field and far field ground motions." Engineering Structures 196: 109292.) and five times (Asghari, Taghipour et al. 2018Asghari, E., R. Taghipour, M. Bozorgnasab and M. Moosavi (2018). "Seismic Analysis of Concrete Gravity Dams Considering Foundation Mass Effect." KSCE Journal of Civil Engineering 22(12): 4988-4996.) dam height. The size which gives the maximum averaging dam responses is concluded as an appropriate near-field domain size in these papers. This evaluation for determining the size of near-filed domains is inaccurate and has no physical interpretation. The size in which the responses in both (WOID) and (WID) conditions are nearly the same or constant can be the appropriate near-field domain size. This evaluation has a physical meaning and it is naturally correct. Moreover, significant errors can occur in the responses when the near-field domain size is assumed to be 3 or 5 times dam height. For example, the tensile stress responses at the heel of the dam in the DRFI problem for (WID, α = 1.5) and (WID, α = 3) conditions are seen to be 16.116 and 43.307 MPa, respectively. The error is calculated as 62.78%. The tensile stress at the heel of the dam for (WID, α = 1.5) and (WID, α = 5) conditions are seen to be 16.116 and 25.995 MPa, respectively. The error is calculated as 38%.

In this paper, the dynamic analysis of a DRI problem selected from (Küçükarslan 2004bKüçükarslan, S. (2004b). "Time‐domain dynamic analysis of dam–reservoir–foundation interaction including the reservoir bottom absorption." International journal for numerical and analytical methods in geomechanics 28(9): 963-980.) took 150 seconds using a Latitude E6430 x64-based computer, while the same problem took 367 seconds to solve in (Küçükarslan 2004bKüçükarslan, S. (2004b). "Time‐domain dynamic analysis of dam–reservoir–foundation interaction including the reservoir bottom absorption." International journal for numerical and analytical methods in geomechanics 28(9): 963-980.). However, the time duration for dynamic analysis of this problem depends on the number of nodes, elements, time steps in the analysis, and method of analysis. By comparing the time durations regarding the solution of the DRI problem, it was seen that the developed numerical model is faster than many other approaches.

6 CONCLUSIONS

In the Lagrangian modeling of DRI, DFI, and DRFI problems an efficient numerical model was developed in FORTRAN 90 programming language, which is more accurate and faster than many approaches. The coupled equation of motion which gives the complex displacement responses of these problems is solved by the direct method in the Laplace domain. Subsequently, the complex displacement responses were transformed into time domain using Durbin’s method. Other dynamic responses were derived utilizing the reel displacement responses. Furthermore, the effect of the near-field domains sizes on the dynamic responses of these problems were investigated including and not including the far-field domains, and an appropriate size for the near-field domains was recommended. The further results can be listed as follow:

• The proposed numerical model is seen to be more accurate and faster than many other methods developed in the literature.

• The size of the near-field domains affects the stress responses. The stress responses for (WOID) and (WID) circumstances, however, were seen to be the same or constant at a certain size. This size was determined to be 4.5 times dam height for the DRI problem and 12 times dam height for the DFI and DRFI problems. This indicates that if the size of the near-field domains is large enough, the reflected waves in the bounded domains can be totally absorbed.

• If the length of the reservoir domain is assumed to be 4 times the height of the dam, there is no need for transmitting boundaries at the far-field reservoir domain in the DRI problem. Similarly, if the size of the near-field domains is assumed to be 12 times dam height, there is no need for transmitting boundaries at the far-field rock domain in the DFI and DRFI problems.

• The stress level obtained at the heel of dam structure for the DFI and DRFI problems in (WOID, α = 12) or (WID, α = 12) conditions is seen to be equivalent to the stress level for (WID, α = 1.5) condition. Furthermore, the responses obtained in these conditions are the appropriate response.

• Significant errors have been detected in the stress responses of the heel point in the DFI and DRFI problems when the near-field domain size is assumed to be 3 or 5 times dam height. The errors in these problems for sizes 3 and 5 were seen to be 62.78% and 38%, respectively.

• The dynamic responses of the handled problems in this study are amplified by increasing the size of the near-field domains. This increasing trend is bigger in (WOID) condition than in (WID) condition.

• By restricting the far-field reservoir boundaries in the DRI problem, the hydrodynamic pressure will be overestimated, while including the proposed infinite fluid elements at far-field reservoir boundaries the hydrodynamic pressure will be underestimated. The difference between responses at the heel of the dam for these conditions was 63%.

• In the DRI and DRFI problems, the hydrodynamic pressure increases due to inclusion of rock foundation in the model. The assumption of rigid base foundation gives inaccurate hydrodynamic pressure response. The difference between hydrodynamic responses at the heel of the dam for these problems was 68%.

In this paper, the dynamic loads generated by horizontal component of Koyna earthquake data are applied to all FE nodes of the model. The static and the hydrostatic loads were not included in the analysis. Only the severe loads are applied to model. As a future work, the findings of this paper can be combined by considering all types of loads in the problem, as well as spatially varying ground motion data.

Edited by

Publication Dates

History