A Novel Method for Modal Analysis of Dam-Reservoir Systems

Civil Engineering Infrastructures Journal 2024, 57(1): 33-59 DOI: 10.22059/CEIJ.2023.349505.1873 RESEARCH PAPER A Novel Method for Modal Analysis of Dam-Reservoir Systems Rezaiee-Pajand, M.1* , Kazemiyan, M. S.2 and Mirjalili, Z.3 1 Professor, Department of Civil Engineering, Ferdowsi University of Mashhad, Mashhad, Iran. 2 Assistant Professor, Department of Civil Engineering, Eqbal Lahoori Institute of Higher Education, Mashhad, Iran. 3 Ph.D. Candidate, Department of Civil Engineering, Ferdowsi University of Mashhad, Mashhad, Iran. © University of Tehran 2023 Received: 05 Oct. 2022; Revised: 03 Feb. 2023; Accepted: 29 May 2023 ABSTRACT: For dynamic modal analysis of the gravity dams, it is required to solve the non-symmetric eigenvalue problem which is a time-consuming process. This paper aims to propose a new procedure for this purpose. In this novel method, there is no need to solve the non-symmetric coupled eigenproblem. Instead, two novel eigenvalue problems are formulated and solved. They are simultaneously applied for dynamic modal analysis of concrete gravity dams. They represent the cubic-symmetric forms of the respective non-symmetric Eigenvalue problem, and they are entitled “cubic ideal-coupled eigenproblems”. Moreover, it is proved that the decoupled and ideal-coupled schemes presented in the previous works can be considered as special cases of the current more general procedure. For solving the aforesaid cubic eigenproblems, the classical subspace algorithm is generalized. To assess the accuracy of the suggested technique, it is employed for the dynamic analysis of two well-known benchmark gravity dams in the frequency domain. The dam crest responses are calculated under upstream and vertical excitations for two different wave reflection coefficients. Then, the obtained results are compared with those of the decoupled and ideal-coupled approaches. Findings corroborate the fact that the authors' formulation is more accurate than the other two aforesaid tactics under all practical conditions. Keywords: Concrete Gravity Dam, Coupled Method, Cubic Ideal-Coupled Method, Decoupled Method, Fluid-Structure Interaction, Ideal-Coupled Method. 1. Introduction The dynamic behavior analysis of a concrete gravity dam-reservoir system can be effectively carried out using the Finite Element-(Finite Element-Hyper Element) technique, commonly abbreviated as FE(FE-HE) (Aftabi Sani and Lotfi, 2010). In * Corresponding author E-mail: [email protected] other words, the dam is modeled with the help of solid finite elements, while the reservoir is divided into a near-field region and a far-field one. The former is near the dam and has an irregular shape, while the latter one is including rectangular strips extending to infinity. These two regions are modeled by the fluid finite element and the 34 fluid hyper-elements, respectively (Aftabi Sani and Lotfi, 2010). Generally, the dynamic analysis process can be performed in the time or frequency domain (Chandravanshi and Mukhopadhyay, 2017). For dam-reservoir systems, many nonlinear constitutive models have been developed for time domain analysis. On another note, the analysis can be conducted in the frequency domain either by the direct approach (Lotfi, 2005) or the modal one. In fluid-structure interaction problems, such as the damreservoir, various alternatives exist for conducting modal analysis (Lotfi, 2005). Some of them utilized the true coupled mode shapes of the system. Within these methodologies, a significant portion of the computational time is dedicated to solve the asymmetric Eigenvalue problem that governs the free vibration behavior of the dam-reservoir systems. To remedy this difficulty, various methods have been proposed for symmetrizing this problem (Rezaiee-Pajand et al., 2021). In general, the earlier ones used extra unknowns (in addition to the pressure) for the fluid domain to symmetrize the problem (Olson and Vandini, 1989). Moreover, they are not efficient, and some of them are not able to calculate the hydrostatic pressure (Everstine, 1981). In this condition, the “decoupled” (Samii and Lotfi, 2007) and “ideal-coupled” (Aftabi Sani and Lotfi, 2010) modal strategies have been proposed to defeat the deficiencies of the aforesaid methods. In the "decoupled" technique, the dam and reservoir Eigen-vectors are separately calculated. They are also applied in the solution procedure instead of the coupled Eigen-vectors. Similarly, the idealcoupled method separately uses the Eigenvectors of the dam and reservoir with modifications in comparison to the decoupled tactic. Based on the related literature, among various techniques symmetrizing the Eigen-problem solved in the dynamic analysis of the dam-reservoir systems, only these two methods use the Eigen-vectors of each domain for Rezaiee-Pajand et al. developing a symmetric version of the originally non-symmetric coupled Eigenproblem. They rely on the mode shapes of two symmetric Eigenvalue problems, which are relatively straightforward from a programming perspective. In a study conducted by Hariri-Ardebili and Mirzabozorg (2013), a direct time-domain approach was proposed for the dynamic stability analysis of the coupled damreservoir-foundation system in threedimensional space. This approach takes into account the impact of the duration of ground motion on the system's dynamic structural stability. Gu et al. (2014) investigated the degradation and safety evaluation of a concrete gravity dam by employing a deterministic and a probabilistic method. Chen et al. (2014) investigated the process of damage and rupture in concrete gravity dams subjected to strong ground motions. Afterward, Lokke and Chopra (2015) suggested a response spectrum analysis strategy estimating the peak response directly from the earthquake design spectrum. Mandal and Maity (2016) proposed a two-dimensional method considering both the fluid-structure and soil-structure interaction in finding the transient response of concrete gravity dams. In another research, Ansari and Agarwal (2017) proposed a new damage index for gravity dams. Furthermore, Guo et al. (2019) used the Lagrange multiplier method for including the dead loads of the arch dam in the dynamic analysis procedure. Moreover, Sotoudeh et al. (2019) conducted a seismic analysis of a system comprising a reservoir, gravity dam, and layered foundation, considering the effects of a vertically propagating earthquake. The methodology developed by Casas and Pavanello (2017) obtained optimal dynamic structural shape through parameter changing, in order to maximize the gap between two adjacent Eigenvalues and also avoid the resonance phenomena at a specific natural frequency interval in Civil Engineering Infrastructures Journal 2024, 57(1): 33-59 coupled fluid-structure systems. Nariman et al. (2019) considered dam-reservoirfoundation interaction and used an extended finite element approach for damage detection of gravity dams. In a study by Liang et al. (2019) a probabilistic analysis was carried out to assess the seismic stability performance of a high arch dam. The analysis method incorporated considerations for contraction joints, boundaries of potential sliding rock masses, and the interaction between the dam and its foundation. Recently, Sotoudehnia et al. (2021) developed an iterative method for reducing the order of the coupled Eigenvalue problem related to fluid-structure interaction systems. This paper aims to introduce a novel modal procedure for the mentioned issue, which is referred to as the cubic-ideal coupled approach. It is considered the enhancement of the “decoupled” and “ideal-coupled” modal techniques. This paper text is structured as follows. Section 2 provides a concise overview of the analysis approach. Then, the coupled, decoupled, and ideal-coupled Eigenproblems are reviewed in Section 3. Afterward, a new cubic-ideal coupled scheme is thoroughly introduced. Furthermore, it is proved that the decoupled and ideal-coupled techniques can be envisaged as special cases of authors' procedures. In Section 4, the cubic eigen vectors are employed for dynamic modal analysis. Section 5 deals with developing a new approach for solving the aforesaid Eigenvalue problem by generalizing the well-known subspace iteration algorithm. In Section 6, the dynamic responses of the triangle ideal dam and Pine Flat dam are achieved by using the special program developed by the authors. Finally, the discussion and conclusions are reported in Sections 7 and 8, respectively. 2. Analysis Method The modal analysis technique is employed in this study (Chandravanshi and Mukhopadhyay, 2017). The FE-(FE-HE) 35 approach is employed to discretize both the dam and fluid domains. For simplicity, the formulation is initially explained without considering the far-field region of the reservoir. Then, the impacts of this region are added to the general case. Therefore, the coupled governing equation of the system takes on the following form (RezaieePajand et al., 2022) : 𝐂 𝟎 𝐫̇ 𝐌 𝟎 𝐫̈ ][ ] + ][ ] + [ [ 𝟎 𝐋 𝐩̇ 𝐁 𝐆 𝐩̈ 𝐓 𝐫 −𝐌𝐉𝐚g [𝐊 −𝐁 ] [𝐩] = [ ] −𝐁𝐉𝐚g 𝟎 𝐇 (1) in which 𝐊, 𝐌, and 𝐂: represent the stiffness, mass, and damping matrices of the dam body, respectively. Furthermore, 𝐇, 𝐆, and 𝐋: correspond to the generalized stiffness, mass, and damping of the fluid domain. Moreover, 𝐁: is the interaction matrix; it emerges in the finite element formulation as a result of vibrating the structure in contact with the water (Aftabi Sani and Lotfi, 2010). The provided matrix establishes a correlation between fluid pressure and structural acceleration. Furthermore, vectors 𝐫 and 𝐩 consist of undetermined nodal displacements and pressures, respectively. It should be added, J is a matrix which each of its two rows are a 2 × 2 identity matrix. It is worthwhile to mention that each column of this matrix is related to a unit rigid body motion in the stream and vertical direction. Additionally, 𝐚g is the vector of ground accelerations. By performing the Fourier transform, the matrix Eq. (1) can be transformed into the following form. [ 𝐫 −ω2 𝐌 + 𝐊(1 + 2βd i) −𝐁 T ] [𝐩] 2 2 −ω 𝐁 −ω 𝐆 + iω𝐋 + 𝐇 −𝐌𝐉𝐚g =[ ] −𝐁𝐉𝐚g (2) where i: represents the imaginary unit and 𝜔: denotes the natural frequency of the system. It is important to note that the provided relation utilizes the hysteretic damping matrix, which takes on the following form (Aftabi Sani and Lotfi, 2010). 36 𝐶= Rezaiee-Pajand et al. 2𝛽𝑑 𝐊 𝜔 (3) where 𝛽𝑑 : denotes the constant hysteretic factor associated with the dam body. It is important to highlight that Eq. (2) represents the coupled equation of a dam within a finite reservoir system in the frequency domain. 3. Free Vibration Analysis It is evident that the Eigenvalue problem corresponding to Eq. (2) can be formulated as follows (Casas and Pavanello, 2017; Rezaiee-Pajand et al., 2023; Sotoudehnia et al., 2021). T 𝐫 𝟎 𝐌 𝟎 (ω2 [ ] + [−𝐊 𝐁 ]) [𝐩] = [ ] 𝟎 𝐁 𝐆 𝟎 −𝐇 (4) Obviously, this linear Eigenvalue problem is similar to the equation governing the free vibration of undamped systems. In contrast, it is not symmetric. To calculate the eigen pairs of the dam-reservoir system, solving this unsymmetrical linear Eigenvalue problem is necessary. Although it is preferred to solve the actual coupled equation of the dam-reservoir system, there are several more efficient alternatives, which will be presented in the following Sub-sections. 3.1. Coupled Eigenproblem Through direct solution of the original eigenvalue problem (4), the actual coupled eigen pairs can be obtained. Usage of the achieved eigen vectors in modal analysis leads to more precise responses when contrasted with other available options. However, standard Eigen-solvers cannot be used to solve the mentioned equation due to their unsymmetrical nature. Based on the studies of other researchers, it has been found that methods for solving unsymmetrical Eigenvalue problems tend to be more time-consuming compared to symmetrical ones. From a programming perspective, they are also more intricate (Aftabi Sani and Lotfi, 2010; Felippa, 1985; Lotfi and Samii, 2012). It should be reminded that the symmetric shapes of the aforesaid Eigenvalue problem can be achieved by introducing new variables. Nevertheless, these extra variables cause complexity in computer programming. 3.2. Decoupled Eigenproblem A symmetrical form of the initial Eigenvalue problem (4) can be achieved by omitting the interaction matrix 𝐁. It is referred to as the "decoupled" form and has the succeeding appearance (Lotfi, 2005): (ω2 [ 𝐫 𝟎 𝐊 𝟎 𝐌 𝟎 ]) [𝐩] = [ ] ]−[ 𝟎 𝟎 𝐇 𝟎 𝐆 (5) The symmetry of the decoupled Eigenvalue problem is quite apparent. This property allows for the utilization of standard Eigen-solvers to efficiently solve the problem. Note that; the eigen vectors obtained from these symmetric equations do not correspond to the actual mode shapes of the real system. Nonetheless, these modes can find application in a modal analysis strategy termed the "decoupled modal approach". It is worth mentioning that the decoupled eigen vectors can be regarded as the Ritz vectors. Therefore, it can be demonstrated that utilizing all of these modes leads to precise solutions. It is important to emphasize that the Eigenvalues obtained from the decoupled Eigenproblem represent the natural frequencies of the dam and reservoir individually (Aftabi Sani and Lotfi, 2010). 3.3. Ideal-Coupled Eigenproblem Herein, the Eigenvalue problems associated with two ideal dam-reservoir systems are solved, rather than the actual coupled system. In the first ideal system, the fluid is considered incompressible, and in the second one, the dam is assumed to be massless. The Eigenvalues obtained from these idealized problems exhibit a higher degree of proximity to the natural frequencies of the real coupled damreservoir system, in contrast to the Eigenvalues derived from the decoupled approach. Additionally, the eigen vectors obtained from these computations are more Civil Engineering Infrastructures Journal 2024, 57(1): 33-59 analogous to the actual mode shapes of the system. These vectors can be employed within a modal analysis approach referred to as the "ideal-coupled modal strategy" (Aftabi Sani and Lotfi, 2010). The first ideal Eigenproblem is presented in a simplified form as follows. [ω2 (𝐌 + 𝐌a ) − 𝐊] 𝐫 = 0 (6) 𝐌a = 𝐁T 𝐇−1 𝐁 (7) 𝐩 = ω2 𝐇 −1 𝐁𝐫 (8) [ω2 (𝐆 + 𝐆a ) − 𝐇]𝐩 = 𝟎 (9) where Ma: represents the added mass matrix and can be obtained as below. Thus, by utilizing Eq. (8), it becomes possible to derive the pressure vector. It is obvious that the dimension of this Eigenproblem corresponds to the number of unknown nodal displacements. The formulation of the second ideal Eigenvalue problem is as follows. in which 𝐆a = 𝐁𝐊 −1 𝐁 T (10) The displacement vector can be computed with the help of the next relation. 𝐫 = 𝐊 −1 𝐁 T 𝐩 (11) Clearly, the dimensions of the second ideal Eigenvalue problem match the count of unknown nodal pressures in the fluid domain. The previously mentioned ideal Eigenproblems can be reformulated as Eq. (12). 𝐌 + 𝐌a (ω [ 𝟎 𝟎 =[ ] 𝟎 2 𝟎 𝐫 𝐊 𝟎 ]) [ ] ]−[ 𝐆 + 𝐆a 𝟎 𝐇 𝐏 (12) This Eigenproblem is a linear and symmetric one. As a result, the solution to this problem can be obtained through the application of commonly used standard methods. Obviously, eliminating 𝐌𝐚 and 37 𝐆𝐚 from Eqs. (6-9) results in the decoupled Eigenvalue problem. Hence, the decoupled version of the actual Eigenproblem represents a specific instance of the idealcoupled Eigenproblem. It is worthwhile to remark that the ideal-coupled approach is more accurate compared to the decoupled one (Aftabi Sani and Lotfi, 2010). 3.4. New Cubic Ideal-Coupled Eigenproblem At this stage, a new symmetric form of the Eigenproblem (4) is introduced. It includes two different cubic Eigenvalue problems, which are separately discussed in this section. It is shown that both decoupled and ideal-coupled strategy can be envisaged as special cases of authors' formulation. Moreover, they are less accurate than the current method. Using the lower partition equation of Eq. (4) and solving the pressure vector in terms of the displacement vector results in the subsequent relation. 𝐩 = ω2 (𝐇 − ω2 𝐆)−1 𝐁𝐫 (13) Obviously, (𝐇 − 𝜔2 𝐆) is the subtraction of two matrices, and it is required to be inverted for calculating the pressure vector. Recall that, Eq. (13) is the exact form of Eq. (8) which plays an important role in the ideal-coupled approach. It is worth mentioning that Eq. (8) can be obtained by removing 𝐆 from Eq. (13). By employing the second-order approximation of the Taylor series, this matrix inversion can be computed as follows (Bakhtiari-Nejad et al., 2005). (𝐇 − ω2 𝐆)−1 ≅ 𝐇 −1 + ω2 𝐇−1 𝐆𝐇−1 +ω4 𝐇−1 𝐆𝐇−1 𝐆𝐇−1 (14) Substituting this relation into Eq. (13) leads to the next result. 𝐩 ≅ ω2 (𝐇 −1 + ω2 𝐇 −1 𝐆𝐇 −1 + ω4 𝐇−1 𝐆𝐇−1 𝐆𝐇 −1 ) 𝐁𝐫 (15) By inserting this equality into the upper partition of Eq. (4), the coming cubic Eigenvalue problem is achieved. 38 Rezaiee-Pajand et al. [ω6 𝐐𝐆𝐇 −1 𝐆𝐐T + ω4 𝐐𝐆𝐐T + ω2 (𝐌 + 𝐌a ) − 𝐊]𝐫 = 𝟎 (16) 𝐒 = 𝐁𝐊 −1 in which 𝐐 = 𝐁T 𝐇 −1 (17) 𝐫 = (𝐊 − ω2 𝐌)−1 𝐁 T 𝐩 (18) Certainly, the dimension of this cubic Eigenproblem is equivalent to the count of unknown nodal displacements in the system. Note that; eliminating the first two terms of Eq. (16) leads to Eq. (6). Accordingly, the first form of the idealcoupled method is a special case of the first cubic ideal-coupled approach. In what follows, the second cubic ideal Eigenproblem is established. To achieve this goal, the displacement vector is solved in terms of the pressure vector by utilizing the upper partition equation of Eq. (4). Consequently, the displacement vector can be computed as below. In fact, Eq. (11) is the approximate form of the last relation, in which 𝐌 is neglected. It should be added that Eq. (11) is one of the key formulas in the ideal-coupled technique. Similarly, (𝐊 − 𝜔2 𝐌) can be inverted with the help of the second-order approximation of the Taylor series. In this way, the succeeding relation can be written (Bakhtiari-Nejad et al., 2005). (𝐊 − ω2 𝐌)−1 ≅ 𝐊 −1 + ω2 𝐊 −1 𝐌𝐊 −1 +ω4 𝐊 −1 𝐌𝐊 −1 𝐌𝐊 −1 (19) Substitution of the aforementioned equation into Eq. (18) yields the next equality. −1 2 −1 −1 𝐫 ≅ (𝐊 + ω 𝐊 𝐌𝐊 + ω4 𝐊 −1 𝐌𝐊 −1 𝐌𝐊 −1 )𝐁T 𝐩 (22) It is clear that the dimension of this cubic Eigenproblem is equivalent to the count of unknown nodal displacements in the system. It is worthwhile to highlight that neglecting the first two terms of Eq. (21) leads to Eq. (9). Hence, the second form of the ideal-coupled scheme is a special case of the second cubic ideal-coupled approach. A 𝑛 × 𝑛 cubic Eigenproblem has Eigenvalues. According to the characteristics of the coefficient matrices, the Eigenvalues may be infinite or finite, and the finite values may be real or complex (Tisseur and Meerbergen, 2001). Obviously, real values are the approximate natural frequencies of the dam-reservoir system, and the other ones are fictitious . The aforesaid two cubic ideal-coupled Eigenvalue problems, i.e., Eqs. (16-21), can be expressed totally as the next shape. (ω6 [ 𝐐𝐆𝐇 −1 𝐆𝐐T 𝟎 𝟎 ] 𝐒𝐌𝐊 −1 𝐌𝐒T 𝐐𝐆𝐐T 𝟎 + ω4 [ ] 𝟎 𝐒𝐌𝐒 T 𝐌 + 𝐌a 𝟎 + ω2 [ ] 𝟎 𝐆 + 𝐆a 𝟎 𝐊 𝟎 𝐫 −[ ]) [𝐩] = [ ] 𝟎 𝟎 𝐇 (23) By solving two separate cubic Eigenvalue problems, the solution of this combined symmetric Eigenproblem can be calculated. It is worthwhile to mention the current relationship can be changed into Eq. (12) by ignoring 𝜔4 𝐐 𝐆 𝐐T, 𝜔4 𝐒 𝐌 𝐒 T, 𝜔6 𝐐 𝐆𝐇−𝟏 𝐆 𝐐T and 𝜔6 𝐒 𝐌 𝐊 −𝟏 𝐌 𝐒 Tterms. 4. Cubic Ideal-Coupled Modal Analysis (20) Introducing this relationship into the lower partition of Eq. (4) leads to the next equation. [ω6 𝐒𝐌𝐊 −1 𝐌𝐒T + ω4 𝐒𝐌𝐒T + ω2 (𝐆 + 𝐆a ) − 𝐇]𝐩 = 𝟎 where (21) Herein, it is assumed that Eigenproblem (23) is solved, and the mode shapes are found. Consequently, the nodal displacements and pressures can be written as follows. 𝒓 𝑿 [ 𝒑] = [ 𝑆 𝟎 𝟎 𝒀𝑆 ][ ] 𝑿𝐹 𝒀𝐹 (24) Civil Engineering Infrastructures Journal 2024, 57(1): 33-59 where 𝐗 S and 𝐗 F : are matrices, which include the eigen vectors of Eigenvalue problems (16) and (21), respectively. These vectors are regarded as the Ritz vectors stemming from the original coupled Eq. (2). They can be employed in combination to provide an approximate solution for the exact problem. It is worth noting that the cubic-ideal coupled mode shapes do not exhibit orthogonality concerning the original mass and stiffness matrices. Nevertheless, the subsequent matrices can be defined as below. 𝑲∗ = 𝑿𝑆 𝑇 𝑲𝑿𝑆 ; 𝑴∗ = 𝑿𝑆 𝑇 𝑴𝑿𝑆 𝑯∗ = 𝑿𝐹 𝑇 𝑯𝑿𝐹 ; 𝑮∗ = 𝑿𝐹 𝑇 𝑮𝑿𝐹 ; 𝑳∗ = 𝑿𝐹 𝑇 𝑳𝑿𝐹 (25) (26) Inserting Eq. (24) into Eq. (2) and performing some simple mathematical operations lead to the next result (Aftabi Sani and Lotfi, 2010). −𝜔2 𝑴∗ + 𝑲∗ (1 + 2𝛽𝑑 𝑖) −𝑿𝑆 𝑇 𝑩𝑇 𝑿𝐹 𝒀 ] [ 𝑆] 𝑇 −2 (−𝜔2 ∗ 𝑮 + 𝑖𝜔𝑳∗ + 𝑯)∗ 𝒀𝐹 −𝑿𝐹 𝑩𝑿𝑆 𝜔 −𝑿𝑆 𝑇 𝑴𝑱𝒂𝑔 = [ −2 𝑇 ] −𝜔 𝑿𝐹 𝑩𝑱𝒂𝑔 [ (27) It is obvious that the vector, which includes the modal participation factors, can be computed with the help of this relation. At each frequency, the response vector can be calculated by introducing the modal participation factor into Eq. (24) in the case of cubic ideal-coupled modal analysis for a dam-finite reservoir system. 4.1. Dam-Reservoir System Up until this point, the formulation of the dam within a finite reservoir was presented. However, when considering a reservoir that extends infinitely, it becomes necessary to incorporate hyper-elements in conjunction with the fluid finite elements. By assembling the hyper-elements matrices, Eq. (2) is converted into the next form. [ −ω2 𝐌∗ + 𝐊 ∗ (1 + 2βd i) −𝐗F T 𝐁𝐗S −𝐗S T 𝐁 T 𝐗F 𝐘 ] [ S] ω−2 (−ω2 𝐆∗ + iω𝐋∗ + 𝐇 ∗ + 𝐗F T 𝐇̄h (ω)𝐗F ) 𝐘F −𝐗S T 𝐌𝐉𝐚g ] = [ −2 T ω 𝐗F (−𝐁𝐉𝐚g + 𝐑̄ p (ω)𝐚g ) (28) 𝐇 (ω) 𝟎 𝐇̄h (ω) = [ h ] 𝟎 𝟎 (29) 𝐑 (ω) 𝐑̄ p (ω) = [ p ] 𝟎 39 (30) 𝐇h (𝜔) and 𝐑p (𝜔) are obtained by expanding 𝐇h (𝜔) and 𝐑 p (𝜔), respectively. These matrices include all pressure degrees of freedom. Note that; Eqs. (27-28) are utilized to determine the vector of participation factors in cases where the reservoir is finite and extends to infinity, respectively. 4.2. Linearized Forms for Solving the Cubic Eigenproblems In the process of the numerical solution of the Standard Eigen Problem (SEP) and the generalized one (GEP), the matrices involved are generally reduced to some simpler forms, which reveal the Eigenvalues. For nonlinear Eigenproblems, these forms cannot be developed. Numerical approaches applied for finding the solution of the cubic Eigenproblems are divided into two categories. The first group directly solves the cubic Eigenproblem, and the second one works with the linearized forms (Afolabi, 1987; Tisseur and Meerbergen, 2001). Note that; most of the numerical tactics, which belong to the first category, are the variants of Newton’s methods whose rate of convergence is highly related to the closeness of the starting guess to the actual solution (Higham and Kim, 2001; Long et al., 2008). These algorithms are able to calculate one Eigen-pair at a time. In practice, it is impossible to guarantee that the scheme converges to the desired Eigenvalue even for an appropriate initial guess. In the techniques based on the linearized forms, a 𝑛 × 𝑛 cubic Eigenproblem is transformed into a 3𝑛 × 3𝑛 linear Eigenvalue problem. In this way, common linear Eigen-solvers incorporated in commercial and non-commercial software packages can be employed. It should be highlighted that the Eigenvalues of a cubic Eigenproblem are similar to their linear form. Furthermore, the eigen vectors can be obtained from the corresponding linear 40 Rezaiee-Pajand et al. problem. Based on the characteristics of the coefficient matrices of a given Eigenvalue problem, various linear forms can be presented for a given cubic Eigenvalue problem. The most important drawback of linearization is that the linearized Eigenproblem's dimension is three times the original cubic one. Based on the linear forms presented in Mackey et al. (2006), suitable symmetric linear forms of the aforementioned cubic ideal-coupled Eigenvalue problems, i.e., Eqs. (16-21), are respectively introduced, as follows. 𝐐𝐆𝐐T − 𝐐𝐆𝐇 −1 𝐆𝐐T 𝐌 + 𝐌a − 𝐐𝐆𝐇 −1 𝐆𝐐T −𝐊 ([𝐌 + 𝐌a − 𝐐𝐆𝐇 −1 𝐆𝐐T 𝐌 + 𝐌a − 𝐐𝐆𝐐T − 𝐊 −𝐊] −𝐊 −𝐊 −𝐊 𝐐𝐆𝐇 −1𝐆𝐐T 𝐐𝐆𝐇 −1 𝐆𝐐T 𝐐𝐆𝐇 −1 𝐆𝐐T 𝐫 2[ T −1 T −1 T (𝐌 − ω 𝐐𝐆𝐇 𝐆𝐐 𝐐𝐆𝐐 + 𝐐𝐆𝐇 𝐆𝐐 − + 𝐌a ) 𝐐𝐆𝐐T + 𝐊 ]) [𝐫̄ ] −1 T T 𝐫̄̄ 𝐌 + 𝐌a + 𝐊 𝐐𝐆𝐇 𝐆𝐐 𝐐𝐆𝐐 + 𝐊 𝟎 = [𝟎] 𝟎 (31) 𝐒𝐌𝐒 T − 𝐒𝐌𝐊 −1 𝐌𝐒 T 𝐆 + 𝐆a − 𝐒𝐌𝐊 −1 𝐌𝐒 T −𝐇 ([𝐆 + 𝐆a − 𝐒𝐌𝐊 −1𝐌𝐒 T 𝐆 + 𝐆a − 𝐒𝐌𝐒 T − 𝐇 −𝐇] −𝐇 −𝐇 −𝐇 𝐩 𝐒𝐌𝐊 −1 𝐌𝐒 T 𝐒𝐌𝐊 −1 𝐌𝐒 T 𝐒𝐌𝐊 −1 𝐌𝐒 T − ω2 [𝐒𝐌𝐊 −1 𝐌𝐒 T 𝐒𝐌𝐒 T + 𝐒𝐌𝐊 −1 𝐌𝐒 T − (𝐆 + 𝐆a ) 𝐒𝐌𝐒 T + 𝐇 ]) [𝐩̄ ] 𝐩̄̄ 𝐆 + 𝐆a + 𝐇 𝐒𝐌𝐊 −1 𝐌𝐒 T 𝐒𝐌𝐒 T + 𝐇 𝟎 = [𝟎] 𝟎 (32) where 𝐩, 𝐩, 𝐫 and 𝐫: contain fictitious entries. It is worthwhile to remark that the dimensions of these problems are equal to three times the unknown nodal pressures and displacements, correspondingly. Obviously, the coefficient matrices are symmetric. As a consequence, these linear Eigenproblems can be easily solved by employing common linear symmetric Eigenvalue solution routines. 5. Generalized Subspace Method Various algorithms have been proposed for estimating the mode shapes and natural frequencies of the linear symmetric Eigenproblems. One of the famous schemes extensively applied is entitled subspace iteration technique developed by Bathe (1996). This method is very popular in the finite element analysis of huge structures (Rezaiee-Pajand et al., 2019). With the help of this procedure, any arbitrary number of structural Eigenvalues and eigen vectors can be approximately calculated. Herein, this well-known approach is generalized for solving the cubic ideal-coupled problems. In each iteration of the generalized approach, a set of vectors is achieved. It should be added that the number of these vectors is less than the size of the initial cubic problem, and the original problem is projected into the corresponding vector space. As a result, a smaller cubic Eigenvalue problem is established. Afterward, it is linearized in an analog manner to the previous subsection. Then, the common linear symmetric Eigenvalue solution routines are utilized for finding the eigen pairs of this smaller problem. Recall that; the obtained responses are the approximations of the Eigenvalues and eigen vectors of the initial cubic eigenproblem. Eventually, the Eigenvalues and eigen vectors of the projected Eigenproblem converge to the eigen pairs of the initial cubic one. It is worth emphasizing that the decoupled mode shapes are applied for establishing the starting set of vectors, which forms the basis of the vector space in the first iteration. In Figures 1 and 2, the steps of this algorithm are proposed for Eigenproblems (16) and (21), respectively. In these flowcharts, 𝑀𝑎𝑥𝐼𝑡𝑒𝑟 and 𝜀: are the maximum allowable iteration and error, correspondingly. 6. Numerical Examples In this study, the finite element method was employed as the initial approach to conduct the analysis. To accomplish this task, a computer program was created by implementing the theories elucidated in the preceding sections. As previously mentioned, solid finite elements were employed to model the dam. Furthermore, the near-field and far-field fluid domains were discretized using fluid finite elements and hyper-elements, respectively. The computer program provides various options for dynamic modal analysis of gravity Civil Engineering Infrastructures Journal 2024, 57(1): 33-59 dams, including the true coupled, decoupled, ideal-coupled, and new cubic ideal-coupled techniques. To solve the Eigenvalue problems with these approaches, different Eigen-solvers are employed. In what follows, the Eigenvalue solution routine of each scheme is introduced. The linear symmetric subspace iteration tactic (Bathe, 1996), which is denoted by SS in the coming sections, is deployed for solving the decoupled and ideal-coupled Eigenproblems. Recall that; the true coupled problem is not symmetric. Hence, its Eigenproblem is solved by the pseudosymmetric subspace iteration strategy abbreviated by PS (Arjmandi and Lotfi, 2011). Moreover, two methods are utilized for the cubic Eigenvalue problems. The first one uses linearization and symmetric subspace iteration strategy (LS), and the second one takes advantage of the suggested Generalized Subspace iteration algorithm (GS). In the subsequent sections, to prove the high accuracy of the proposed new method, it is utilized for conducting dynamic analysis of the ideal triangle and Pine Flat gravity dams in the frequency domain. In these examples, the dynamic responses of the dam crests are calculated in response to both upstream and vertical excitations. This analysis takes into account two different values of wave reflection coefficients (𝛼), specifically 1 and 0.5. It should be reminded that 𝛼 = 1 represents the full reflection and 𝛼 = 0.5 allows for the partial reflection of waves, which influences the reservoirfoundation boundaries (Bougacha and Tassoulas, 1991; Jafari and Lotfi, 2018). In each case, the amplitude of the complexvalued accelerations for the dam crest point is plotted versus the dimensionless frequency 𝜔/𝜔1𝑠 . It should be added that 𝜔 and 𝜔1𝑠 denote the excitation frequency and the first frequency of the dam on the rigid foundation with no water in the reservoir, respectively. The results obtained from the analysis are then compared with the exact solutions, which are derived using a direct 41 method that incorporates all the true coupled mode shapes. Additionally, the same comparison is conducted for the decoupled and ideal-coupled approaches. Moreover, the accuracy and consumed time of the above-cited Eigen-solution routines in finding the eigen pairs are compared. For this purpose, the next efficiency and error indices are introduced. 𝑇𝐼𝑖 = 100 × 𝑇𝑚𝑖𝑛 𝑇𝑖 (33) 𝑛𝑚 𝑗 𝑗 − 𝑓𝑖 | 1 |𝑓 𝐸𝐼𝑖 = 100 × ( ∑ 𝑒𝑥𝑎𝑐𝑡 ) 𝑗 𝑛𝑚 𝑓 𝑗=1 𝑒𝑥𝑎𝑐𝑡 (34) where the consumed time of the fastest Eigen-solution routine and the i-th one are demonstrated by 𝑇 𝑚𝑖𝑛 and 𝑇𝑖, 𝑗 correspondingly. Furthermore, 𝑓𝑖 and 𝑗 𝑓𝑒𝑥𝑎𝑐𝑡 : are the jth natural frequency of the ith tactic and the true coupled one, respectively. Besides, nm: denotes the number of computed natural frequencies. 6.1. Ideal Triangle Gravity Dam In this subsection, the mentioned methods are utilized for the dynamic analysis of a famous gravity dam named the ideal triangle gravity dam in the frequency domain. In what follows, the finite element model and basic parameters of this system are introduced, and the obtained results are presented. 6.1.1. Model At this stage, the focus is on the finite element model of the ideal triangle dam situated on a rigid foundation. To represent this dam, a discretization technique is employed, utilizing 20 isoparametric 8node plane-solid finite elements. As it was previously mentioned, the water domain includes near-field and far-field regions. The former one continues up to a specific length (𝐿), which is measured at the dam crest point in the upstream direction. Herein, it is assumed that 𝐿 = 0.2𝐻. It should be added that 𝐻 is the dam height or maximum water depth in the reservoir. 42 Rezaiee-Pajand et al. Fig. 1. Flowchart of generalized subspace method for Eigenproblem (16) Following the near-field region, the farfield portion commences and stretches to infinity in the upstream direction. The nearfield region is simulated using 5 isoparametric 8-node plane-fluid elements, while the far-field segment is modeled with a fluid hyper-element consisting of 5 isoparametric 3-node sub-elements. It is worthwhile to mention that the used mesh pattern has been previously applied by other researchers (Sotoudehnia et al., 2021; Ziaolhagh et al., 2016). Figure 3 depicts the finite element model of the ideal triangle dam and its reservoir. Civil Engineering Infrastructures Journal 2024, 57(1): 33-59 Fig. 2. Flowchart of generalized subspace method for Eigenproblem (21) (a) Dam with the near-field fluid region (b) Dam with the near far-field fluid regions Fig. 3. The finite element model of the ideal triangle gravity dam 43 44 Rezaiee-Pajand et al. 6.1.2. Basic Parameters The concrete dam is presumed to be homogenous with isotropic linearly viscoelastic behavior. Its elastic modulus, unit weight, and Poisson's ratio are 27.5 Gpa, 24.8 kN/m3 and 0.2, respectively. Additionally, the hysteretic damping factor is 0.05. Furthermore, the impounded water is assumed to be irrotational, compressible, and inviscid, with a pressure wave velocity of 1440 m/s and a unit weight of 9.81 kN/m³. 6.1.3. Free Vibration Responses It should be reminded that each aforementioned formulation includes two cases whose Eigenproblems are not similar. Consequently, each method has two sets of modes, except for the true coupled technique. It is worth emphasizing that the set of mode shapes associated with the nodal displacements can be computed by solving the first Eigenvalue problems, and the corresponding set of mode shapes related to the nodal pressures are calculated by solving the second Eigenvalue problems. Accordingly, the frequencies of the first and second cases are listed in Tables 1 and 2, respectively. The provided tables reveal that the natural frequencies of the true coupled problem generally appear to be lower than the two sets of natural frequencies computed in each instance of the decoupled, ideal-coupled, and cubic ideal-coupled approaches. For comparison, Figure 4 illustrates the error indices of the decoupled, ideal-coupled, cubic idealcoupled, and true coupled approaches. With the help of this figure, the accuracy of the aforesaid tactics can be compared. It is obvious that the natural frequencies of the cubic ideal-coupled approach are more accurate compared to the decoupled and ideal-coupled methods. In other words, the most accurate tactic is the authors' technique, and the error-index of the decoupled method is higher than those of others. Table 1. The first five natural frequencies for the ideal triangle dam-reservoir system with 𝐿 = 0.2𝐻 Natural frequencies fi (Hz) Cubic idealDecoupled Ideal-coupled coupled Mode True number dam First ideal case coupled First cubic ideal (Ziaolhagh et al., (Incompressible fluid case 2016) assumption) 1 2.29 1.49 1.28 1.25 2 5.19 4.08 3.44 2.54 3 6.04 5.94 5.90 4.96 4 8.93 7.84 6.68 5.65 5 13.17 11.29 9.96 6.13 Table 2. The second five natural frequencies for the ideal triangle dam-reservoir system with 𝐿 = 0.2𝐻 Natural frequencies fi (Hz) Cubic idealDecoupled Ideal-coupled Mode coupled True Number coupled Reservoir Second ideal case Second cubic (Ziaolhagh et al., (Incompressible fluid ideal case 2016) assumption) 1 1.80 1.35 1.26 1.25 2 5.40 3.61 2.92 2.54 3 9.03 7.09 5.89 4.96 4 12.76 10.45 9.59 5.65 5 16.66 12.68 10.05 6.13 Civil Engineering Infrastructures Journal 2024, 57(1): 33-59 45 Fig. 4. Error index In this example, the efficiency indices of the decoupled approach with SS, the idealcoupled approach with SS, and the cubic ideal-coupled approach with GS and LS are equal to 100, 100, 100, and 4.49, respectively. For the truly coupled algorithm with PS, the efficiency index is 22.70. Clearly, the efficiency indices of the first three ones are the same, and the cubic ideal coupled with LS is the slowest one. Note that; the ratio of the efficiency index to the error-index is a key parameter for comparing the performance of numerical techniques. In other words, the performance of a numerical method is dependent on both its efficiency and accuracy. Accordingly, the cubic ideal-coupled technique with GS performs more successfully in comparison to the other schemes in this numerical example. 6.1.4. Forced Vibration Responses At this stage, the magnitudes of the complex acceleration values for the dam crest point are plotted versus the dimensionless frequency 𝜔/𝜔1𝑠 . To achieve this goal, the aforementioned strategies are used. It is important to note that the number of modes to perform the dynamic analysis in each case is the same. Recall that; other researchers previously presented the responses by using the true coupled tactic (Hojati and Lotfi, 2011; Samii and Lotfi, 2011). For 𝛼 = 1, the outcomes are illustrated in Figures 5 and 6 for the upstream and vertical excitations, respectively. For all two types of excitations considered, it can be observed that the cubic ideal-coupled approach performance is mainly better than the decoupled and idealcoupled strategies. Now, the results for 𝛼 = 0.5 are illustrated in Figures 7 and 8 for the upstream and vertical excitations, respectively. Obviously, for these cases, the cubic ideal-coupled scheme's responses are also closer to the exact response (i.e., the direct method with the true coupled mode shapes) in comparison to the decoupled and idealcoupled techniques' results. 6.2. Pine Flat Gravity Dam Herein, the mentioned strategies are applied for conducting a dynamic analysis of the Pine Flat gravity dam in the frequency domain. Subsequent sections will provide information about the finite element model, the essential parameters of this system, and the obtained results. 6.2.1. Model The finite element model of the Pine Flat dam on a rigid foundation has been examined. The dam is discretized using 40 isoparametric 8-node plane-solid finite elements. It is important to note that the water domain includes near-field and farfield regions. In this example, 𝐿 = 200 𝑚. The far-field portion initiates after the nearfield region and extends infinitely in the upstream direction. For modeling the nearfield region, 90 isoparametric 8-node plane- 46 Rezaiee-Pajand et al. fluid elements are utilized, while the farfield section is represented by a fluid hyperelement consisting of 9 isoparametric 3node sub-elements. It is worth mentioning that the mesh pattern used here has been previously employed by other researchers (Ganji and Lotfi, 2021; Omidi and Lotfi, 2017). Figures 9 and 10 provide a visualization of the finite element model of the Pine Flat dam and its corresponding reservoir. Horizontal Ground Motion Horizontal Acceleration 30 𝛼=1 20 10 0 0 1 2 3 𝛚/𝛚𝐬𝟏 (a) Decoupled method 4 5 4 5 4 5 Horizontal Ground Motion Horizontal Acceleration 30 𝛼 20 10 0 0 1 2 3 𝝎/𝝎𝒔𝟏 (b) Ideal-coupled method Horizontal Acceleration Horizontal Ground Motion 30 𝛼 20 10 0 0 1 2 3 𝝎/𝝎𝒔𝟏 (c) Cubic ideal-coupled method Fig. 5. Frequency response function at the dam crest resulting from horizontal excitation with 𝛼 = 1 Civil Engineering Infrastructures Journal 2024, 57(1): 33-59 47 Vertical Ground Motion Horizontal Acceleration 40 𝛼=1 30 20 10 0 0 1 2 𝛚/𝛚𝐬𝟏 3 4 5 4 5 (a) Decoupled method Vertical Ground Motion Horizontal Acceleration 40 𝛼=1 30 20 10 0 0 1 2 𝛚/𝛚𝐬𝟏 3 (b) Ideal-coupled method Vertical Ground Motion Horizontal Acceleration 40 𝛼=1 30 20 10 0 0 1 2 𝛚/𝛚𝐬𝟏 3 4 5 (c) Cubic ideal-coupled method Fig. 6. Frequency response function at the dam crest resulting from vertical excitation with 𝛼 = 1 48 Rezaiee-Pajand et al. Horizontal Ground Motion Horizontal Acceleration 30 𝛼 = 0.5 20 10 0 0 1 2 𝛚/𝛚𝐬𝟏 3 4 5 (a) Decoupled method Horizontal Ground Motion Horizontal Acceleration 30 𝛼 = 0.5 20 10 0 0 1 2 𝛚/𝛚𝐬𝟏 3 4 5 (b) Ideal-coupled method Horizontal Ground Motion Horizontal Acceleration 30 𝛼 = 0.5 20 10 0 0 1 2 𝛚/𝛚𝐬𝟏 3 4 5 (c) Cubic ideal-coupled method Fig. 7. Frequency response function at the dam crest due to horizontal excitation with 𝛼 = 0.5 Civil Engineering Infrastructures Journal 2024, 57(1): 33-59 49 Vertical Ground Motion Horizontal Acceleration 20 𝛼 = 0.5 10 0 0 1 2 3 𝛚/𝛚𝐬𝟏 (a) Decoupled method 4 5 Vertical Ground Motion Horizontal Acceleration 20 𝛼 = 0.5 10 0 0 1 2 𝛚/𝛚𝐬𝟏 3 4 5 4 5 (b) Ideal-coupled method Vertical Ground Motion Horizontal Acceleration 20 𝛼 = 0.5 10 0 0 1 2 3 𝛚/𝛚𝐬𝟏 (c) Cubic ideal-coupled method Fig. 8. Frequency response function at the dam crest resulting from vertical excitation with 𝛼 = 0.5 50 Rezaiee-Pajand et al. Fig. 9. Dam body with the near-field fluid region Fig. 10. Dam body with the near-field and far-field fluid regions It is important to note that the Pine Flat dam features a sloped upstream face. To enhance accuracy, the hyper-elements need to be connected to the vertical sides of the finite elements. Consequently, within the finite element region, this slope should gradually decrease before establishing connections between the hyper-elements and the finite elements. 6.2.2. Basic Parameters The dam body is constructed from homogeneous concrete with isotropic linearly viscoelastic behavior, possessing an elasticity modulus of 22.75 Gpa, unit weight of 24.8 kN/m3, and a Poisson's ratio of 0.2. Moreover, a hysteretic damping factor of 0.05 is considered for the material. In addition, the impounded water is treated as irrotational, compressible, and inviscid, having a unit weight of 9.81 kN/m3 and a pressure wave velocity of 1440 m/s. 6.2.3. Free Vibration Responses At first, the frequencies of the first and second Eigenproblems of this damreservoir system are proposed in Tables 3 and 4, respectively. Clearly, the natural frequencies of the true coupled problem are generally lower than the corresponding sets of natural frequencies obtained using the decoupled, ideal-coupled, and cubic ideal-coupled approaches in each case. Furthermore, the natural frequencies of the cubic strategy are closer to those of the true coupled one in comparison to the other approaches. At this stage, the first and second pressure mode shapes are demonstrated in Figures 11 and 12, respectively. Recall that; the true coupled mode shapes were previously proposed in other works (Samii and Lotfi, 2007). Table 3. The first five natural frequencies for the Pine Flat dam-reservoir system with 𝐿 = 200 𝑚 Natural frequencies fi (Hz) Decoupled Cubic idealMode Ideal-coupled True coupled (Samii and Lotfi, coupled number (Samii and 2007) Lotfi, 2007) First ideal case (Incompressible First cubic Dam fluid assumption) ideal case 1 3.15 2.67 2.58 2.53 2 6.48 5.77 4.95 3.27 3 8.74 8.66 8.45 4.67 4 11.25 10.35 9.27 6.22 5 16.99 15.98 13.51 7.92 Table 4. The second five natural frequencies for the Pine Flat dam-reservoir system with 𝐿 = 200 𝑚 Natural frequencies fi (Hz) Decoupled Cubic ideal(Samii and Mode Ideal-coupled True coupled coupled Lotfi, 2007) number (Samii and Second ideal case (Incompressible Second Cubic Lotfi, 2007) Reservoir fluid assumption) ideal case 1 3.12 2.94 2.68 2.53 2 4.75 4.24 3.52 3.27 3 7.80 6.05 5.01 4.67 4 9.30 7.92 7.28 6.22 5 9.96 9.46 8.89 7.92 Civil Engineering Infrastructures Journal 2024, 57(1): 33-59 (a) Decoupled (c) Cubic ideal-coupled Fig. 11. First pressure mode shapes (a) Decoupled 51 (b) Ideal-coupled (d) True coupled (b) Ideal-coupled (c) Cubic ideal-coupled (d) True coupled Fig. 12. Second pressure mode shapes 52 Rezaiee-Pajand et al. Obviously, the mode shapes of the cubic ideal-coupled are more similar to true coupled ones. For brevity, the dam mode shapes are not presented. However, the scheme of this paper is more successful in calculating these mode shapes than the other tactics. For the aforesaid dam-reservoir system, Figures 13 and 14 show the error and efficiency indices of the decoupled, idealcoupled, cubic ideal-coupled and true coupled strategies, correspondingly. With the help of these figures, the accuracy and analysis duration of the aforesaid schemes can be compared. As it was previously mentioned, the ratio of the efficiency index to the error-index is a key parameter for comparing the performance of numerical techniques. Accordingly, the ideal-coupled scheme and cubic ideal-coupled technique with GS perform more successfully in comparison to other algorithms. Obviously, the natural frequencies of the cubic ideal-coupled method are more closely aligned with the true coupled frequencies compared to the decoupled and ideal-coupled approaches. Consequently, if the same number of modes is employed, the cubic ideal-coupled approach is expected to offer improved accuracy in dynamic response compared to the decoupled and ideal-coupled techniques. Similarly, the decoupled and ideal-coupled tactics are faster than the other algorithms. The cubic ideal coupled with GS technique is ranked as second. Besides, the cubic ideal coupled with GS is much faster than the true coupled with PS, and the cubic ideal coupled with LS is the slowest one. Fig. 13. Error index Fig. 14. Time index Civil Engineering Infrastructures Journal 2024, 57(1): 33-59 6.2.4. Forced Vibration Responses In a similar manner to the previous example, the chart depicts the changes in the magnitudes of complex-valued accelerations at the crest of the dam in relation to the dimensionless frequency 𝜔/𝜔1𝑠 . When considering 𝛼 = 1, the outcomes are illustrated in Figures 15 and 53 16 for the upstream and vertical excitations, respectively. Analogously, the aforesaid methods utilize the same number of modes to perform the dynamic analysis in each case. In previous research, true coupled responses were proposed (Chopra et al., 1980). Horizontal Acceleration Horizontal Ground Motion 30 𝛼=1 20 10 0 0 1 2 3 4 5 𝛚/𝛚𝐬𝟏 (a) Decoupled method Horizontal Ground Motion Horizontal Acceleration 30 𝛼=1 20 10 0 0 Horizontal Acceleration 30 1 2 3 4 5 3 4 5 𝝎/𝝎𝒔𝟏 (b) Ideal-coupled method Horizontal Ground Motion 𝛼=1 20 10 0 0 1 2 𝛚/𝛚𝐬𝟏 (c) Cubic ideal-coupled method Fig. 15. Frequency response function at the dam crest resulting from horizontal excitation with 𝛼 = 1 54 Rezaiee-Pajand et al. Vertical Ground Motion Horizontal Acceleration 30 𝛼=1 20 10 0 0 1 2 𝛚/𝛚𝐬𝟏 3 4 5 (a) Decoupled method Vertical Ground Motion Horizontal Acceleration 30 𝛼=1 20 10 0 0 1 2 𝛚/𝛚𝐬𝟏 3 4 5 4 5 (b) Ideal-coupled method Vertical Ground Motion Horizontal Acceleration 30 𝛼=1 20 10 0 0 1 2 𝛚/𝛚𝐬𝟏 3 (c) Cubic ideal-coupled method Fig. 16. Frequency response function at the dam crest resulting from vertical excitation with 𝛼 = 1 Civil Engineering Infrastructures Journal 2024, 57(1): 33-59 55 At this stage, the results for 𝛼 = 0.5 are depicted in Figures 17 and 18 for the upstream and vertical excitations, correspondingly. For all two types of excitations considered, the cubic ideal-coupled scheme performs more successfully than the decoupled and ideal-coupled approaches. Horizontal Ground Motion Horizontal Acceleration 20 𝛼 = 0.5 10 0 0 1 2 𝛚/𝛚𝐬𝟏 3 4 5 (a) Decoupled method Horizontal Ground Motion Horizontal Acceleration 20 𝛼 = 0.5 10 0 0 Horizontal Acceleration 20 1 2 3 4 5 3 4 5 𝛚/𝛚𝐬𝟏 (b) Ideal-coupled method Horizontal Ground Motion 𝛼 = 0.5 10 0 0 1 2 𝛚/𝛚𝐬𝟏 (c) Cubic ideal-coupled method Fig. 17. Frequency response function at the dam crest resulting from horizontal excitation with 𝛼 = 0.5 56 Rezaiee-Pajand et al. Vertical Ground Motion Horizontal Acceleration 20 𝛼 = 0.5 10 0 0 1 2 3 𝛚/𝛚𝐬𝟏 (a) Decoupled method 4 5 Vertical Ground Motion Horizontal Acceleration 20 𝛼 = 0.5 10 0 0 1 2 𝛚/𝛚𝐬𝟏 3 4 5 4 5 (b) Ideal-coupled method Vertical Ground Motion Horizontal Acceleration 20 𝛼 = 0.5 10 0 0 1 2 𝛚/𝛚𝐬𝟏 3 (c) Cubic ideal-coupled method Fig. 18. Frequency response function at the dam crest resulting from vertical excitation with 𝛼 = 0.5 Civil Engineering Infrastructures Journal 2024, 57(1): 33-59 Obviously, for these cases, the responses obtained using the cubic ideal-coupled approach demonstrate a higher level of agreement with the exact responses obtained through the direct method when compared to the results obtained from the decoupled and ideal-coupled methods. 7. Discussion One of the challenges which exists in the dynamic analysis of dam-reservoir systems is solving the corresponding non-symmetric Eigenproblem. Two well-known approaches, which are used for symmetrizing this problem, are decoupled and ideal-coupled methods. This paper has presented a novel method that is more accurate than both methods. However, although it is faster than the true coupled approach, it is not as fast as the decoupled and ideal-coupled are. Hence, further research activities are suggested for developing more accurate and faster methods, in comparison to the decoupled and ideal-coupled tactics. 57 considering two different values of wave reflection coefficients. The obtained results were compared with those of the decoupled and ideal-coupled strategies. By thoroughly investigating the findings, it is concluded that: • The novel approach can find more accurately the forced and free vibration responses of the gravity dams in comparison to the decoupled and idealcoupled approaches. This is because its eigen pairs are closer to those of the true coupled ones. In other words, this paper suggested a modal dynamic analysis strategy that is more accurate than the other aforementioned available ones. • Moreover, it is observed that the cubic ideal-coupled scheme with the suggested Eigen-solver algorithm is faster than the true coupled one with the pseudosymmetric method while it requires more time in comparison to the decoupled and ideal-coupled techniques, which are less accurate than authors' tactic. 9. References 8. Summary and Conclusions In this paper, a novel frequency-domain approach for performing modal analysis of concrete gravity dam-reservoir systems was presented. This method was developed based on two cubic Eigenvalue problems. To solve them, the well-known subspace algorithm was generalized. Moreover, their solution can be found with the combination of the classical subspace scheme with linearization. To achieve this goal, the linearized forms of the aforesaid Eigenproblems were proposed. This tactic was utilized for the dynamic analysis of two famous gravity dams, namely the ideal triangle dam and the Pine Flat Dam. It is worth noting that the far-boundary condition of the reservoir was considered by employing the hyper-elements. Furthermore, the dynamic responses of the dam crests were calculated in response to both upstream and vertical excitations, Afolabi, D. (1987). “Linearization of the quadratic Eigenvalue problem”, Computers and Structures, 26(6), 1039-1040, https://doi.org/10.1016/0045-7949(87)90120-9. Aftabi Sani, A., and Lotfi, V. (2010). “Dynamic analysis of concrete arch dams by ideal-coupled modal approach”, Engineering Structures, 32(5), 1377-1383, https://doi.org/10.1016/j.engstruct.2010.01.016. Ansari, M.D.I., and Agarwal, P. (2017). “Damage index evaluation of concrete gravity dam based on hysteresis behavior and stiffness degradation under cyclic loading”, International Journal of Structural Stability and Dynamics, 17(01), 1750009. https://doi.org/10.1142/S0219455417500092. Arjmandi, S.A., and Lotfi, V. (2011). “Computing mode shapes of fluid-structure systems using subspace iteration methods”, Scientia Iranica A, 18(6), 1159-1169, https://doi.org/10.1016/j.scient.2011.09.011. Bakhtiari-Nejad, F., Rahai, A. and Esfandiari, A. (2005). “A structural damage detection method using static noisy data”, Engineering Structures, 27(12), 1784-1793, https://doi.org/10.1016/j.engstruct.2005.04.019. Bathe, K.J. (1996). Finite Element procedure, 2nd 58 Edition, Prentice-Hall, INC. Bougacha, S. and Tassoulas, J.L. (1991). “Seismic response of concrete gravity dams, II: Effects of sediments”, Journal of Engineering Mechanics Division, ASCE, 117(8), 1839-1850, https://doi.org/10.1061/(ASCE)07339399(1991)117:8(1839). Casas, W.J.P. and Pavanello, R. (2017). “Optimization of fluid-structure systems by Eigenvalues gap separation with sensitivity analysis”, Applied Mathematical Modelling, 42, 269-289, https://doi.org/10.1016/j.apm.2016.10.031. Chandravanshi, M.L. and Mukhopadhyay, A.K. (2017). “Modal analysis of a vertically tapered frame”, International Journal of Structural Stability and Dynamics, 17(03), 1771001, https://doi.org/10.1142/S0219455417710018. Chen, H.Q., Li, D.-Y. and Guo, S.-S. (2014). “Damage–rupture process of concrete dams under strong earthquakes”, International Journal of Structural Stability and Dynamics, 14(07), 1450021, https://doi.org/10.1142/S0219455414500217. Chopra, A.K., Chakrabarti, P. and Gupta, A.S. (1980). Earthquake response of concrete gravity dams including hydrodynamic and foundation interaction effects, University of California, California. Everstine, G.C. (1981). “A symmetric potential formulation for fluid structure interaction”, Journal of Sound and Vibration, 79(1), 157-160, https://doi.org/10.1016/0022-460X(81)90335-7. Felippa, C.A. (1985). “Symmetrization of the contained compressible-fluid vibration Eigenproblem”, Communications in Applied Numerical Methods, 1, 241-247, https://doi.org/10.1002/cnm.1630010509. Ganji, H.T., and Lotfi, V. (2021). “Transient analysis of concrete gravity dams by Wavenumber-TD method for general excitation”, Proceedings of the Institution of Civil Engineers-Structures and Buildings, 174(4), 259-275, https://doi.org/10.1680/jstbu.19.00041. Gu, Q., Yu, C., Lin, P., Ling, X., Tang, L. and Huang, S. (2014). “Performance assessment of a concrete gravity dam at Shenwo reservoir of China using deterministic and probabilistic methods”, International Journal of Structural Stability and Dynamics, 14(05), 1440002, https://doi.org/10.1142/S0219455414400021. Guo, S., Liang, H., Li, D., Chen, H. and Liao, J. (2019). “A comparative study of cantilever-and integral-type dead loads on the seismic responses of high arch dams”, International Journal of Structural Stability and Dynamics, 19(03), 1950021, https://doi.org/10.1142/S0219455419500214. Hariri-Ardebili, M.A., and Mirzabozorg, H. (2013). Rezaiee-Pajand et al. “A comparative study of seismic staibility of coupled arch dam-foundation-reservoir systems using infinite elements and viscous boundary models”, International Journal of Structural Stability and Dynamics, 13(6), 1350032, https://doi.org/10.1142/S0219455413500326. Higham, N.J. and Kim, H.M. (2001). “Solving a quadratic matrix equation by Newton’s method with exact line searches”, SIAM Journal on Matrix Analysis and Applications, 23(2), 303316, https://doi.org/10.1137/S0895479899350976. Hojati, M. and Lotfi, V. (2011). “Dynamic analysis of concrete gravity dams utilizing twodimensional modified efficient fluid hyperelement”, Advances in Structural Engineering, 14(6), 1093-1106, https://doi.org/10.1260/13694332.14.6.1093. Jafari, M. and Lotfi, V. (2018). “Dynamic analysis of concrete gravity dam-reservoir systems by a wavenumber approach for the general reservoir base condition”, Scientia Iranica, 25(6), 30543065. https://doi.org/10.24200/SCI.2017.4227. Liang, H., Guo, S., Tu, J. and Li, D. (2019). “Seismic stability sensitivity and uncertainty analysis of a high arch dam-foundation system”, International Journal of Structural Stability and Dynamics, 19(07), 1950066, https://doi.org/10.1142/S0219455419500664. Lokke, A. and Chopra, A.K. (2015). “Response spectrum analysis of concrete gravity dams including dam-water-foundation interaction”, Journal of Structural Engineering, ASCE, 141(8), 1-9, https://doi.org/10.1061/(ASCE)ST.1943541X.0001172. Long, J.H., Hu, X.Y. and Zhang, L. (2008). “Improved Newton’s method with exact line searches to solve quadratic matrix equation”, Journal of Computational and Applied Mathematics, 222, 645-654, https://doi.org/10.1016/j.cam.2007.12.018. Lotfi, V. (2005). “Significance of rigorous fluidfoundation interaction in dynamic analysis of concrete gravity dams”, Structural Engineering and Mechanics, 21(2), 137-150, Lotfi, V. and Samii, A. (2012). “Frequency domain analysis of concrete gravity dam-reservoir systems by wavenumber approach”, In: The 15th World Conference on Earthquake Engineering, Lisbon. Mackey, D.S., Mackey, N., Mehl, C. and Mehrmann, V. (2006). “Vector space of linearizations for matrix polynomials”, SIAM Journal on Matrix Analysis and Applications, 28(4), 971-1004, https://doi.org/10.1137/050628350. Mandal, K.K., and Maity, D. (2016). “Transient response of concrete gravity dam considering dam-reservoir-foundation interaction”, Journal of Earthquake Engineering, 22(2), 211-233, Civil Engineering Infrastructures Journal 2024, 57(1): 33-59 https://doi.org/10.1080/13632469.2016.121780 4. Nariman, N.A., Lahmer, T. and Karampour, P. (2019). “Uncertainty quantification of stability and damage detection parameters of coupled hydrodynamic-ground motion in concrete gravity dams”, Frontiers of Structural and Civil Engineering, 13(2), 303-323, https://doi.org/10.1007/s11709-018-0462-x. Olson, L.G. and Vandini, T. (1989). “Eigenproblems from finite element analysis of fluid-structure interactions”, Computers and Structures, 33(3), https://doi.org/10.1016/0045679-687, 7949(89)90242-3. Omidi, O. and Lotfi, V. (2017). “A symmetric implementation of pressure-based fluid– structure interaction for nonlinear dynamic analysis of arch dams”, Journal of Fluids and Structures, 69, 34-55, https://doi.org/10.1016/j.jfluidstructs.2016.12.0 03. Rezaiee-Pajand, M., Aftabi Sani, A. and Kazemiyan, M.S. (2019). “An efficient Eigen-solver and some of its applications”, International Journal for Computational Methods in Engineering Science and Mechanics, 20(2), 130-152, https://doi.org/10.1080/15502287.2018.153415 0. Rezaiee-Pajand, M., Esfehani, S. and Ehsanmanesh, H. (2022). “An explicit and highly accurate Runge-Kutta family”, Civil Engineering Infrastructures Journal, 56(1), 51-78, https://doi.org/10.22059/CEIJ.2022.330788.179 2. Rezaiee-Pajand, M., Kazemiyan, M.S., and Aftabi Sani, A. (2021). “A literature review on dynamic analysis of concrete gravity and arch dams”, Archives of Computational Methods in Engineering, vol?? Issue?? 1-16, https://doi.org/10.1007/s11831-021-09564-z. Rezaiee-Pajand, M., Mirjalili, Z., and Kazemiyan, M.S. (2023). “Analytical 2D model for the liquid storage rectangular tank”, Engineering Structures, 289, 116215, https://doi.org/10.1016/j.engstruct.2023.116215. Samii, A. and Lotfi, V. (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, 1003-1012, https://doi.org/10.1016/j.finel.2007.06.015. Samii, A. and Lotfi, V. (2011). “High-order adjustable boundary condition for absorbing evanescent modes of waveguides and its application in coupled fluid-structure analysis”, Wave Motion, 49(2), 238-257, https://doi.org/10.1016/j.wavemoti.2011.10.001. Sotoudeh, P., Ghaemian, M. and Mohammadnezhad, H. (2019). “Seismic analysis of reservoir-gravity dam-massed layered 59 foundation system due to vertically propagating earthquake”, Soil Dynamics and Earthquake Engineering, 116, 174-184. https://doi.org/10.1016/j.soildyn.2018.09.041. Sotoudehnia, E., Shahabian, F. and Sani, A.A. (2021). A dynamic order reduction method for fluid-structure systems”, Applied Mathematical Modelling, 89, 136-153, https://doi.org/10.1016/j.apm.2020.06.071. Tisseur, F. and Meerbergen, K. (2001). “The quadratic Eigenvalue problem”, SIAM Review, 43(2), 235-286, https://doi.org/10.1137/S0036144500381988. Ziaolhagh, S.H., Goudarzi, M. and Aftabi Sani, A. (2016). “Free vibration analysis of gravity damreservoir system utilizing 21 node-33 Gauss points triangular elements”, Coupled Systems Mechanics, 5(1), 59-86, https://doi.org/10.12989/csm.2016.5.1.059. This article is an open-access article distributed under the terms and conditions of the Creative Commons Attribution (CC-BY) license.