Eigenvalue solution for arch dams: ADAD-IZIIS Software

3D earthquake analyses of complex engineering structures such as arch dams are based on the determination of a range of eigenvalues and associated eigenvectors. Despite the variety of practical numerical eigensolutions available in the literature, the method of subspace iteration was found to offer time efficient and accurate eigenvalue solutions appropriate for large systems with high degrees of freedom. This paper presents the theoretical principles of the well known subspace iteration method implemented in the existing ADAD-IZIIS software for finite-element analyses of arch dams.


Introduction
Eigenvalues and eigenvectors are inherent characteristics of any engineering structural model. These parameters also referred to as 'mode shapes and natural frequencies' represent major dynamic features of the structure. In many cases, the mode shapes and natural frequencies of structures are the basis of the design criteria. Their evaluation is important: in the planning process during the pre-construction stage for accurate seismic design in the post-construction period for the long term monitoring of structural seismic properties, vital to help engineers improve the safety and maintainability of critical structures.
In addition, an update of the model stiffness parameters, obtained by comparing the calculated eigenvalues with those measured by ambient or force vibration tests, is important to provide insight in the current health status of the structure. Eigenvalues and eigenvectors [1][2][3], have many applications in both pure and applied engineering mathematics. The eigenproblems related to steady-state and dynamic analyses of structures are of particular interest in the civil engineering practice. In general, dynamic analyses of structures can be based on the mode superposition method for linear problems, or on the direct integration of uncoupled equation of dynamic motion for nonlinear problems. The application of either one of these techniques requires solution of the eigenproblem. The complex structural systems usually generate large matrices such that it becomes impractical to evaluate all eigenvectors and eigenvalues that theoretically participate in the response. This is especially important in the computation of the linear elastic dynamic response based on the mode superposition method. Knowing that for rigid structures, such as gravity dams, only the first several eigenvalues and associated eigenvectors are significant to the seismic response [4][5][6][7][8], it is necessary to define an optimal number of modes that should be included for an accurate modal analyses. This number varies depending on the type of the structure, type of dynamic loading, and desired accuracy. It should satisfy the criterion of the cumulative effective mass participation greater than 90 % of the total structural mass. An optimum number of modes also depends on the frequency content of the selected excitation. The three dimensional (3D) dynamic analyses of important complex flexible structures, such as arch dams, usually require the determination of a larger number of eigenvalues and eigenvectors [9][10][11][12][13]. This is especially true when the applied excitation has a higher level of power spectral density and a wide range of dominant frequencies of the Fourier amplitude spectrum. The application of the direct integration method for solving dynamic motion of dams constitutes an imperative for the analysis of nonlinear behavior [11,[14][15][16]. In this case, the eigensolution is required as the direct integration method is based on the application of the Rayleigh damping concept. The Rayleigh damping coefficients, a and b, are determined using known damping ratios in uncoupled mode, x i , along with the corresponding circular frequencies, ω i , chosen to represent the energy dissipation ability of the structural system. For a given coupled dam-reservoir-foundation system, the Rayleigh damping coefficients provide equivalent damping ratios capable to cover correctly the range of modes (frequencies) of interest. The simulated energy dissipation efficiency can be measured by plotting the equivalent modal damping as a function of the eigenvalues for the entire eigenspace [17]. The selection of new values is needed for cases when the representation of the effective damping ratios is not sufficiently accurate. To overcome this problem, Indrajit et al. [17] suggested an iterative solution for the best-fit values of a and b [17]. Therefore, the definition of a wider range of eigenvalues is essential even in cases of nonlinear dynamic analyses. The objective of the present study is to provide a critical stateof-the-art literature review of the available practical numerical eigensolutions and to select, implement, and validate a method that is most appropriate for large systems with high degrees of freedom using the existing ADAD-IZIIS software for the finiteelement analysis of arch-dams [18].

Eigenvalue solution methods
The fundamentals and basic considerations used for the solution of a typical eigenproblems have been extensively presented in literature [1][2][3]. The standard eigenvalue formulation can be solved using the linear equation of the following form: where, x is a nonzero vector that satisfies the above equation, and l represents the eigenvalue of the positively definite matrix A. The eigenvalue l is the factor by which the eigenvector x changes when multiplied by the matrix A. An eigenspace of a matrix is the set of all pairs of eigenvectors and eigenvalues. The objective of the dynamic response analyses of an assemblage of structural elements is the solution of the generalized eigenvalue problem written in the following form: When the finite element discretization is used for domains for which the general form of eigenvalues is sought, the matrices A and B are positively definite and represent the stiffness and mass matrix, respectively. Both matrices are of the same order and, for cases when the mass matrix is not diagonal, both have same bandwidth. Experience has shown that consistent mass formulation is not always necessary and that a good accuracy can also be obtained using a lumped mass matrix [19].
The buckling instabilities of structures also constitute an eigenvalue problem [20] requiring definition of only the first several eigenvalues. In this case, the standard linear equation has the following form: GRAĐEVINAR 70 (2018) 10, 881-890 Eigenvalue solution for arch dams: ADAD-IZIIS Software where K is the deflection stiffness matrix, K G is the geometric stiffness matrix corresponding to the matrix B in equation (2), l is the buckling load, and u is the corresponding buckling mode. In this case, the eigenvalue problem is defined with the general equation (2). However, to obtain the required solution, equation (2) needs to be transformed as follows: where the target eigenvalue solution is the absolute maximum value of χ = 1/l that gives the lowest buckling load l. The consideration of fluid compressibility in hydrodynamics also involves application of the standard eigenvalue formulation. Following the transformation, the second order Poisson differential equation assumes the following form [21][22]: where p = p(x,y,z,t) is the hydrodynamic pressure distribution in space and time. In the frequency domain, it can be written as p = p o sin(ωt), where the second derivative is -p o ω 2 sin(ωt). This formulation leads to the standard eigenvalue definition: It can be concluded from the above discussion that eigensolution is required for numerous engineering problems. Several numerical methods are available in the current practice. The selection of the most appropriate method to solve a particular problem relies first of all on the system size, matrix bandwidth and the number of required eigenvalues and eigenvectors. In addition, the efficiency of the solution involving parameters such as: accuracy, required number of numerical operations, storage requirements, and central processor unit (CPU) time, should also be considered.

Low degree of freedom systems
The existing efficient and accurate methods for eigensolution of systems with low degrees of freedom and small bandwidth can be regrouped into three major groups: vector iteration methods (inverse iteration, forward iteration, Gram-Schmidt iteration, shifting in vector iteration, Rayleigh quotient iteration, Jacobi iteration) transformation methods (Jacobi method, Hausholder-QRinverse iteration method) explicit and implicit polynomial iteration methods.
Among the above methods, the well-known Householder-QR-inverse iteration technique [2] merits attention as the generally accepted method for solving eigenvalues based on standard formulation of the problem without the need for its transformation. The generalized Jacobi iteration is also an interesting and often considered technique. It provides a direct solution for all eigenvalues and eigenvectors of the established generalized eigenvalue problem [22][23][24][25][26]. In addition, the static condensation of the massless degrees of freedom in a lumped mass analysis is not essential for a direct solution and it can be particularly efficient when the off-diagonal elements are small or sparse [22,26].

High degree of freedom systems
Two basic numerical techniques have been developed for eigensolution of large systems with high degrees of freedom: the Determinant Search Method [24][25][26][27] and the Subspace Iteration Method [22,28,29]. Both techniques use some of the methods applied for the low degree of freedom systems, either in the process of transformation or in the process of iteration required for finding solution to an eigenproblem. The algorithm of Subspace Iteration Method, originally developed by Bathe [22], has recently been improved by parallel processing enabled by the use of new generation of computers [29].

Determinant search method
The determinant search is a widely used and a more advanced method. It effectively incorporates the triangular factorization and the Gram-Schmidt inverse vector iteration method. This combination offers an efficient solution in the case of large systems where the number of required eigenvalues and corresponding vectors is much smaller than the order of matrices, and in cases where matrices have a small bandwidth. The method is based on the fact that eigenvalues actually represent the roots of the characteristic polynomial. The applied procedure first shifts the solution in the vicinity of the subsequent unknown root and then uses the inverse iteration to determine the eigenpair. The number of negative pivots in the course of factorization indicates whether the shift is in accordance with the required root. The convergence of the applied iterations can be sometimes slow, especially when iterations yield convergence towards a multiple root or a cluster of roots. In such cases, the eigenvalue separation theorem (Sturm sequence property) is applied to accelerate the iterative procedure [30][31]. The required eigenvalues and vectors are then obtained in succession from the least dominant eigenpair upwards. Two serious setbacks accompany the determinant search method under certain circumstances. The first one is the failure to assess the accuracy of the approximations for the eigenvectors. The second setback is the possibility of losing the lower eigenvalues and eigenvectors and inaccurate estimation of the higher eigenvalues which are mostly congested. The GRAĐEVINAR 70 (2018) 10, 881-890 Violeta Mircevska, Miroslav Nastev, Viktor Hristovski, Ivana Bulajic method's inefficiency is amplified by the need of one triangular factorization for each determinant evaluation, which needs additional computational time.

Subspace iteration method
The subspace iteration method solves the generalized form of the eigenvalue problem to directly define the first several eigenvalues, indicated herein with p, and the associated eigenvectors. In opposition to the determinant search method, there is no need for transformation into the standard eigenvalue form. The subspace iteration method gives accurate and time effective solutions for a wide range of eigenvalues and is convenient for large bandwidth systems with hundreds of degrees of freedom. To decrease the computational time, it uses primarily vector iterations and only a few factorizations. The method's main feature is the reduction of dimensionality of the analysed problem from N-dimensional space to q = 2p dimensional subspace for solving the q dimensional eigenvalue problem, where N is the number of degrees of freedom of the analysed domain and p is the number of least dominant sought eigenvalues with the required accuracy ( where p < q << N ). Generally, the convergence can be obtained in only a few iterations when the starting subspace spanned by the q vectors is a "good" approximation of the least dominant p-dimensional subspace (i.e., when the initial transformation vectors are close to the sought eigenvectors). It has also been shown that a higher number of starting iteration vectors contributes towards spanning a larger initial subspace and a higher ultimate convergence rate of the iteration vectors. Therefore, the matrix X in the generalized eigenvalue equation comprises q eigenvectors to guarantee monotonic convergence of the p least dominant eigenvectors. The reasonable number of iteration vectors leading to monotonic convergence is defined as q = min{2p, p+8}. The technique is based on the simultaneous iteration of a desired number of vectors and of operators A and B, taking advantage of the fact that both tend towards the diagonal form within subspace iterations. Obviously, the main difficulty of the subspace iteration method lies in the selection of the starting subspace. An optimum eigensolution algorithm should also prove numerical stability in generating an orthogonal basis in each subspace within the smallest number of required simultaneous iterations that still provide a satisfactory convergence rate. A simultaneous iteration over the generalized eigenvalue formulation assumes solution of the following equation: where the matrix X k+1 stores iteration vectors after performing the first k iteration steps, and the matrix R k+1 is the upper triangular matrix which ensures that the vectors in X k+1 are B-orthogonal.
Different iteration schemes exist for finding solution to equation (7). Bauer [32] proposed an original bi-iteration method for an arbitrary form of the matrix A. On the other hand, Rutishauser [34] focused on solving a special case of the symmetric matrix A. The simultaneous iteration method introduced by Jennings is based on calculation of the linear prediction matrix for improved approximations by means of the current iteration vectors [34][35][36]. Although these advanced iteration schemes incorporate different techniques, they all represent a convenient combination of basic methods such as: the Ritz analysis that provides the best approximation of the lowest eigenvalues [37], the Rutishauser subspace iteration algorithm [33], or the method of inverse iteration with Gram-Schmidt orthogonalisation.

Algorithm for subspace iteration method
There are several reasons why the generalized Jacobi method was applied in the ADAD-IZIIS program [18] for solving a generalized form of the eigenvalue equation. It is a rotational method for elimination of off-diagonal terms based on the threshold iteration. In the process, the off-diagonal terms are zeroed only if they are a magnitude smaller than the threshold for the current iteration. The coupling factors and are used to represent the threshold within each iteration as a measure of the coupling between generalized degrees of freedom i<j. The coupling factors are assumed to be smaller than 10 -10 within the fifth iteration and the threshold is set to 10 -2k , where k is the iteration number. The relative change in the eigenvalue estimates is set to be smaller than 10 -8 . The main difficulty of the subspace iteration method lies in the selection of initial transformation vectors as a starting subspace. Hence, the choice of the initial transformation vectors X represents an important step in the subspace iteration method and a challenging topic for scientific investigations. Only several schemes for selecting the starting vectors in the p-dimensional subspace have been proposed in the literature, and so this selection still remains highly subjective. The general approach consists in establishing the matrix X columns using elements from the matrices A and B only [22]. The matrix R = B·X, representing the right-hand side of the equation A·X = l·B·X, was found to be the most effective when the first vector in the matrix X is simply a diagonal of the matrix B. The remaining vectors are unit vectors with +1 at the coordinate with the largest a ii /b ii ratio. The rationale behind the selection of the above type of starting vectors is to avoid potential missing of a mode by exciting all the mass degrees of freedom in the first vector of matrix B. At the same time, the remaining vectors must be linearly independent and should excite points of maximum mass and flexibility ratio only. Moreover, to ensure better convergence, the unit entries of the second to the last vector should not be very close together. The algorithm based on the above approach is closely related to the static condensation analysis. However, it may sometimes lead GRAĐEVINAR 70 (2018) 10, 881-890 Eigenvalue solution for arch dams: ADAD-IZIIS Software to a significant number of iterations. Another potential option is to include the conventional Ritz analysis in which the specified load patterns can represent a reasonable first approximation. The algorithm proposed by Uhrig [38] based on the component mode synthesis is also sometimes applied.

Step by step solution
The applied algorithm comprises the following steps over the generalized eigen formulation (2).
Step 1 -Definition of the starting eigensubspace: the eigenspace is first determined with q=min{2p, p+8}. For cases when 2p<p+8, the eigensubspace is spanned by the starting iteration vectors q, where q=2p. These vectors constitute the matrix R 1 of the order of N·2p with the following form: where B is the diagonal mass matrix of the order of N·N; X 1 is the initial transformation matrix (starting iteration vectors) spanning the eigensubspace q. The order of the matrix X 1 is N·2p and it is defined according to the Bathe's scheme [22]; and p is desired number of eigenvalues and vectors.
Step 2 -Solution of the matrix : triangular factorization is applied to solve the matrix using the Gaussian elimination method over the equation (2). In the method, the matrix represents the starting matrix for the next iteration k = 1, 2, 3, .... ∞ (9) where the operator A remains in an unaltered form (original stiffness matrix) of the order of N·N, is the new iteration matrix of the order of N·2p, and R k is a previously defined iteration matrix of the order of N·2p.
Step 3 -Transformation vectors: according to the applied technique, the iteration vectors assume the role of transformation vectors over the operators, matrices A and B. The transformation actually represents a projection of these matrices onto solution space e k+1 , as follows: The above transformation reduces the degrees of freedom from N to 2p, and eliminates the coordinates out of the selected smaller q subspace. The order of the operators A k+1 and B k+1 is 2p·2p.
Step 4 -Solution of the eigensystem: the solution of the eigensystem of the projected operators using the generalized eigenformulation and generalized Jacobi method is as follows: A k+1 · Q k+1 = B k+1 · Q k+1 · l k+1 (11) where, the operators A k+1 and B k+1 are defined in the q subspace with a reduced number of degrees of freedom. Therefore, any of the available methods for eigensolution with high or low degrees of freedom could be applied in this step. The generalized Jacobian rotational method for elimination of offdiagonal terms is applied in this case. The required eigenvectors approximations are calculated in a single step. The matrix Q k+1 comprises the eigenvectors of the subspace q and is of the order of 2p·2p. With the iteration vectors approaching the eigenvectors, the operators A and B shift to a diagonal form.
Step 5 -Solution of the R k+1 matrix: the next step involves determination of improved approximation of the eigenvectors estimates and generation of a new iteration matrix R k+1 used in the consecutive iteration. (12) where the revised solution of the eigenvector matrix, equation 7, is of the order of N·2p. The improved iteration matrix R k+1 of the same order is obtained as follows: Step 6 -Testing accuracy of the solution: the eigenvalue approximations l of the last two subsequent iterations k and k+1 are used to test convergence of the iterative process as follows: < tolerance, for i = 1, 2, ..., p The usual error tolerance is 10 -6 . Although the iterations are performed for q vectors, q>p, the convergence is measured only for approximations obtained for the first p eigenvalues.
Step 7 -Repeat steps 2 to 6: the numerical procedure is repeated until the convergence within the assigned accuracy is achieved. The final eigensolution can be expressed as follows: i , for i = 1, 2, ..., p (15) To ensure that the lowest p eigenvalues and vectors have been determined, the Sturm sequence check is applied at the end of the iterative process based on the eigenvalue separation theorem. The check comprises computation of eigenvalues histograms for symmetric three-diagonal matrices where the number of eigenvalues is less than the assumed cut-off eigenvalue [30,31]. GRAĐEVINAR 70 (2018) 10, 881-890 Violeta Mircevska, Miroslav Nastev, Viktor Hristovski, Ivana Bulajic

ADAD-IZIIS software
The ADAD-IZIIS software [18] is used to conduct an analytical procedure for the three-dimensional dynamic analysis of arch dams including the effects of dam-water interaction (water incompressibility), soil-structure interaction, and the nonlinear effect of contraction joint opening. A powerful feature of the program is the option for computer design of the dam body following the exact topology of the canyon where the dam is built in a process of connection and embedment of arches in the terrain [39][40]. The program automatically generates: the finite element mesh of the dam along with a portion of the foundation mass to account for the dam-foundation interaction, and the boundary element mesh of the fluid domain boundaries to account for the fluid-structure interaction [40]. The magnitude of hydrodynamic pressures at the dam-fluid interface is dependent on the amount of energy transmitted to the fluid by vibration of the dam and the surrounding terrain [40][41]. The generation of the combined finite element and boundary element mesh is sufficiently accurate for replicating the terrain topology [40].
The step-by-step solution for subspace iteration eigensolution described in the previous section was developed for this study in an algorithm written in Fortran language. It was then embedded in the ADAD-IZIIS software. The software has an option to use a multithread processing engaging the MKL Intel and OMP libraries. However, the parallel processing in case of eigensolution is not treated in this paper.

Eigenvalues and eigenvectors of a complex arch dam system
To test its capabilities, the ADAD-IZIIS software with the embedded subspace iteration algorithm was used to evaluate the eigenvalues and eigenvectors for a double-curved arch dam.
The validation was carried out by comparing the simulation results with those obtained with the commercially available SAP 2000 v.14 software [42].  Figure 2. The mathematical model of the dam, shown in Figure  3, is generated automatically, using defined shapes of the arcs at all selected elevations in accordance with the topology of terrain. The dam body is composed of 27 monolithic blocks. The FE model of the arch dam consists of 199 substructures further discretized into 1592 finite elements. Twenty and fifteennode solid elements were used in the automatic discretization process. The built model contains 3,805 external nodes located between the substructures and has 11,415 degrees of freedom. The ADAD-IZIIS software considers banded matrices as a special type of sparse matrices. The matrix band could contain some zero elements. In the presented analyses, the system of mass and stiffness matrices has a bandwidth of 1220 for the particular case of 11,415 degrees of freedom. The ADAD-IZIIS software uses a special algorithm for bandwidth minimization Eigenvalue solution for arch dams: ADAD-IZIIS Software of 0.85). The fundamental period of the free vibration of the dam (empty reservoir) is Te=0.318s, while it amounts to Tf=0.345s for the coupled dam-fluid system. Accordingly, the mass of water in the reservoir interacts with the solid structure during the vibration, increasing the length of the vibration period (eigenvalues) of the system. In both ADAD-IZIIS and SAP 2000 v.14 software, the hydrodynamic pressure (HDP) was calculated according to the added mass principle [43] using the following equation: (16) where H is the water level in the reservoir and z is the coordinate of the point on the dam face measured in the coordinate system in which H is defined. Table 1 presents comparison of the eigenvalues, i.e. the vibration periods of the empty and full reservoirs. The ratios between the assessments of the vibration periods computed by both ADAD-IZIIS and SAP 2000 v.14 software packages were obtained by means of the following simple equations: R ADAD/SAP = 100 · (T ADAD -T SAP )/T SAP , and Ref = 100 · Te/Tf. It can be observed from Table 1 that the R ADAD/SAP ratios vary in the narrow range between 0 and 10% for the empty reservoir and the reservoir filled with water. R ADAD/SAP ratios different from zero value are more frequent for comparisons with the full reservoir, which is most probably due to the fact that the eigenvalues were calculated outside of the ADAD-IZIIS software using the added mass concept. Such results are more than encouraging because they show high agreement, and ratios different from zero are merely an artefact resulting from the numbers rounded to the nearest hundredth of a second rather than from any numerical inaccuracy. On the other hand, to obtain a matrix with a lower bandwidth. Large and complex structures, such as dams, are usually associated with large bandwidth matrices even if the minimization process is applied. If we express the banded matrix in terms of density, the analysed system exhibits a density of 21 The dynamic behaviour of the coupled arch dam-reservoir system was significantly affected by the fluid-dam interaction and, therefore, the vibration periods of the entire system had to be defined. Since the ADAD-IZIIS software provides the finiteelement boundary element (FEM-BEM) orientated solution of the FSI effects, and as the fluid domain is not discretized with fluid finite elements, the added mass method [41] proved to be the only possibility to evaluate the eigenvalues of the coupled system. The added mass concept is effective and sufficiently accurate for calculation of the eigenvalues of the coupled damreservoir systems. The eigenvalues for the coupled system were calculated for the impounded water depth of 110m (water depth to dam height ratio  Eigenvalue solution for arch dams: ADAD-IZIIS Software the Ref ratios are, as expected, within the 80 and 90% range for both codes. Possible discrepancies coincide with those for the R ADAD/SAP ratios. The generated FE models and obtained vibration patterns and mode shapes for the first, second, and tenth mode, are given in Figures 4 to 6, respectively. The presented FE models of the dam are slightly different. They mismatch at the end of the arches, as the visualization software included with the ADAD-IZIIS code draws arches along the whole length without exclusion of the embedded part. This, however, does not affect significantly the accuracy of the eigensolution as the free vibration patterns are associated only with the degrees of freedom, with no specified restraints.
The few characteristic mode shapes which give a preview of the vibration patterns obtained by both software packages are also given in Figures 4 to 6. It can be observed that the mode shapes coincide with each other.

Conclusion
The computational algorithm of the subspace iteration method embedded in the ADAD-IZIIS software for finite-element analyses of arch dams is presented in this paper. The method is especially efficient in defining a wider range of eigenvalues and eigenvectors for complex systems with large stiffness and mass matrices and with hundreds of degrees of freedom. The algorithm is based on solution of the generalized form of the eigenvalue problem directly without the need for transformation into the standard eigenvalue form. The main features of this method are reduction of space dimensionality and simultaneous iteration of a desired number of vectors and corresponding operators benefiting from the matrix diagonalization during the iteration process.
The main advantage of the applied subspace iteration technique is its time effectiveness and accurate approximation of the modes, practically excluding the possibility of missing the dominant mode. The adapted code in the ADAD-IZIIS software was tested on a complex double-curved arch dam to evaluate eigenvalues and eigenvectors. The validation procedure was carried out by comparing simulation results with those obtained with the commercially available SAP 2000 v.14 software. Only negligible discrepancies were observed in the lengths of vibration periods of the considered first 12 modes.