WYD method for an eigen solution of coupled problems

Designing efficient and stable algorithm for finding the eigenvalues and eigenvectors is very important from the static as well as the dynamic aspect in coupled problems. Modal analysis requires first few significant eigenvectors and eigenvalues while direct integration requires the highest value to ascertain the length of the time step that satisfies the stability condition. The paper first presents the modification of the well known WYD method for a solution of single field problems: an efficient and numerically stable algorithm for computing eigenvalues and the corresponding eigenvectors. The modification is based on the special choice of the starting vector. The starting vector is the static solution of displacements for the applied load, defined as the product of the mass matrix and the unit displacement vector. The starting vector is very close to the theoretical solution, which is important in cases of small subspaces. Additionally, the paper briefly presents the adopted formulation for solving the fluid-structure coupled systems problems which is based on a separate solution for each field. Individual fields (fluid and structure) are solved independently, taking in consideration the interaction information transfer between them at every stage of the iterative solution process. The assessment of eigenvalues and eigenvectors for multiple fields is also presented. This eigen problem is more complicated than the one for the ordinary structural analysis, as the formulation produces non-symmetrical matrices. Finally, a numerical example for the eigen solution coupled fluidstructure problem is presented to show the efficiency and the accuracy of the developed algorithm.


INTRODUCTION
In dynamic analysis of structures such as: dams, underwater structures, shore and off shore structures and similar, it is necessary to simulate the fluid-structure interaction to ascertain the real behaviour of such a complex system.This problem is commonly referred to as a coupled problem.The initial development in dynamic coupled (multi field) problems had taken place in aerospace and nuclear industry, and has lately expanded to all aspects of engineering, including civil engineering [1,2] The most universal, most powerful and the most widely applied engineering tool for this analysis are undoubtedly numerical methods [1,9]: finite element technique, boundary element technique, mesh less techniques or some other.Although there are different approaches to solving the fluid-structure coupled problems, the most natural formulation is the Lagrange-Euler formulation with structure's unknown node displacements and unknown fluid's pressures as degrees of freedom.
The dynamic coupled problems are often nonlinear and thereby very laborious and expensive processes and consequently there has been much effort directed.This problem is still current and recently published papers [4,5,6,7,8] are a testament of intensive research in this field.
Eigen problem solution in dynamic analysis is important in determining the system's dynamic characteristics [1, 15,16,18].Solving dynamic problems through model analysis is based on knowing the system's first few eigenvalues and eigenvectors.Direct integration methods require the highest or the lowest value to ascertain the length of the time step that satisfies the stability condition (explicit methods) or accuracy (implicit methods).
The problem of determining eigenvalues is still a field of intensive research [9,14].There is no universal method that could solve the eigenvalues problems with adequate speed and accuracy.Choosing the method depends on a number of factors: number of the system's degrees of freedom, number of eigenvalues to be calculated, shape and matrices density and where in the eigen spectrum the eigenvalues of interest are located [9,15].There are also cases when some methods for determining eigenvalues fail [9,15].Some are focused on computing all eigenvalues and vectors, some pass over near eigenvalues.Some algorithms are applicable only for symmetrical matrices and so on [15].
The greatest need for efficient and effective algorithms exists in the area of solving large eigen problems, those of the order of several thousand degrees of freedom and greater, such as coupled problems.For such large eigen problems it is computationally prohibitive to calculate the exact system eigen solution.Many of the methods used to solve smaller problems are not feasible.This has resulted in the development of approximate eigen solution techniques that provide only a partial eigen solution.
To date many procedures have been developed for computing partial eigen solutions of such large eigen problems [12], for example: simultaneous iteration method, subspace iteration method, determinant search method, Lanczos method, Ritz vector method etc.
Each of these methods have their advantages and disadvantages.The following is a presentation of the modified WYD method which is essentially a variant of the Ritz vector method [13].

EIGEN SOLVING TECHNIQUE
The WYD method for computing eigenvalues and eigenvectors was first published in [13].It is an efficient algorithm for determining the first "k" eigenvalues/vectors, where "k" is a freely chosen number.The WYD method does not compute the eigenvalues/vectors but transforms the system so that another method can be applied, such as the Jacobi method, vector iteration method and similar.The bases of the numerical procedure is looking for the solution in only one subspace, which is significantly faster than iteration across subspaces [11,12,13].The procedure produces multiple, 2k times, static solutions of the problem thus forming the Ritz base vectors, which is practical for programming on microcomputers.The problem of searching eigenvalues is reduced from a large n × n problem (n being several million in large problems) to a 2k × 2k problem which greatly reduces the number of computational operations and the size of the accumulated calculation error.
The WYD method is stable and reliable, eigen values are not skipped which is something that occurs with a squeeze procedure vector iteration.Generally, for the required "k" eigenvalues/vectors, "2k" Ritz vectors are needed.In the solution the first "k" vectors are exact and the other "k" are approximations [12].
Eigen problem of the structural dynamics is given in: [11,12,13] where K is the well known system stiffness matrix and M is the well known system mass matrix.The procedure for forming the "2k" Ritz space is as follows [13]: 1. Computation of the first Ritz vector x 1 : where x 0 is the initial vector that preconditions the first vector.That is, x 1 extracts as the first vector the one closest to it by being the next perpendicular to Mx 0 .If we start from the beginning, i.e. we look for the lowest eigen vector, then value x 0 is chosen as the vector with unit components.The method is then referred to as the modified WYD method [11].This selection insures that the eigen values are not passed over if the analyzed dynamic system is asymmetrical.If the system is symmetrical, anti-symmetrical eigen vectors and corresponding values can't be recorded.This flaw can be rectified by deliberately transforming the symmetrical system into an asymmetrical one by small modifications.M-norming follows: (3) 2. Computation of the other Ritz vectors x i (i = 1, 2, ..., 2k): determining the constants c j (j = 1, 2, ..., i-1) and determining the new vector orthogonal to the previous (Gramm-Schmidt procedure): and its M-norming: (7) 3. K-orthogonalization of X Ritz vectors and forming the projective subspace: with: where is a full matrix and E is unit matrix.This produces a standard eigen problem: that can be solved applying the e.g.Jacobi method.Eigenvalues of this "compressed" problem are 2k eigenvalues of the initial problem (where the first "k" values are exact and the other "k" approximations).Eigenvectors of the initial problem can be obtained from: where X is the Ritz vector matrix (n × 2k), and Q the eigen vector matrix obtained in the projective subspace.

SHORT DESCRIPTION OF COUPLED PROBLEM
As previously emphasized, in dynamic analysis of structures such as dams, underwater structures, shore and off shore structures and similar, it is necessary to simulate the groundfluid-structure interaction.In these coupled problems, the structure's behaviour is under the fluid's constant influence which leads to the conclusion that these two mediums have to be analyzed interdependently.
Behavior of the fluid-structure interaction problem can be expressed with general, well known second order differential equations in the matrix form [1]: Equation ( 12) can be written as: (13) where: x 1 , x  13) and requires large computing resources.This problem is usually solved through the so called partitioned solution scheme where the combined system is divided into single fields that are then observed independently taking into consideration the interaction forces on contact surfaces (here the system's global matrices are symmetrical thus simplifying the system computation).The procedure is iterative until the convergence criterion is satisfied.This approach is not applicable for solving eigen problems where the entire system has to be analyzed.
By using the fluid pressure formulation and structure's displacement formulation, the behaviour of the fluid-structure system can be described with a system of two second order differential equations: that define the dynamic equilibrium of the system, with: M s i M f are the structure's and fluid's mass matrices, f cs is the fluid on structure forces vector and f cf is the structure on fluid forces vector, Q is the interaction matrix.
Equation system ( 14) can be written in matrix form: (16) making it evident that the global system's mass and stiffness matrices are asymmetrical.
Eigen coupled problem can be defined from equation ( 16), with an analogous expression, as: Equation ( 17), can be written as: where: In accordance with postulations given in chapter 2, we are looking for the lowest eigenvector, so x 0 is chosen as the vector with unit components.
Fluid-structure interaction surface with fluid and structure elements is shown in Figure 1.Interaction matrix Q includes only the surface integration and is defined as: (20) N ui i N pj are the structure's and fluid's base function matrices, and → n is the unit outer norm on the interaction surface shown in Figure 1.

EXAMPLE
An analysis of the Grančarevo dam in Bosnia and Hezegovina (Figure 2) has been performed to show the efficiency of the method.The Grančarevo Arch Dam is a double-curvature concrete dam with perimetral joint.The height of the dam is 123 meters and the crest length is 439 meters.Its bottom thickness is 27 meters and its top thickness 4.60 meters.The dam's foundation dig is 230.000m 3 and the volume of poured concrete is 376.000m 3 .The head of the dam is 100 meters.The dam created the Bile ća reservoir with the maximum water depth of 51 meters and the available storage capacity of 1100 million cubic meters.The Bileća reservoir is the largest storage lake on Balkan.Its dimensions are: total storage volume: 1280 hm 3 and surface of the reservoir on normal top water level: 2764 ha [19,20,21,22].
Figure 3 shows the plan of the dam's body with land topology.Geometrical characteristics of the dam's body were based on the data in lit.[20,21,22].Geometrical data tables show basic geometrical characteristics for individual arches some of which are shown in Figure 4.
Figure 5 shows the finite element mesh of the structure (the dam and the surrounding rock) and the accumulated water.For the discretization of the dam and the terrain 27 "brick" elements with three degrees of freedom in every node (node displacements) were used and for the accumulated water 27 "brick" elements with one degree of freedom in every node (node pressures).Due to the dam's geometry, finite elements on the edges of the dam are significantly deformed.The dam and the surrounding terrain were modeled with 2720 elements and 27329 nodes and the fluid with 1068 elements and 11324 nodes.
Table 1 shows the material characteristics of the concrete and water used in the analysis.The dam was analyzed for the maximum water level, 8 m below the dam's head [20,21].
Table 2 shows the values of eigen frequencies for the first five modes recorded during experimental testing (C1) [21]       Figure 6 shows the first four eigenvectors that are the result of the modified WYD method analysis.It is evident that they match the results in literature [19], with the first mode frequency being nearly exact.The second mode slightly deviates, whilst the third and the fourth are close matches.

CONCLUSION
As presented the modified WYD method can be easily applied for solving large eigen tasks in coupled fluid/structure problems, particularly for a small number of eigenvalues and vectors.The procedure universally and automatically generates the starting vectors and the required size of the subspace.There are vectors present in all the needed directions for the required number of eigenvectors.Solution convergence and high accuracy is assured.Eigenvectors are not skipped.In cases of system symmetry or if identical properties occur in two orthogonal directions, corresponding paired solutions are produced as a single solution.A small difference in the properties results with a separation of a single solution.
The method is very fast, it doesn't requires large computing resources and is easily programmable.

Figure 3
Figure 3 Plan of dam's body with land topology [20].

Figure 6
Figure 6  Eigenvalues and eigenvectors for the dam-water-foundation rock system.
• 1 , ẍ1 are displacement, velocity and acceleration vectors, M 11 , C 11 , K 11 are mass, dampening and stiffness matrices, and f 1 is the outer nodal forces vector of the first field, x 2 , x • 2 , ẍ2 , M 22 , C 22 , K 22 , f 2 are corresponding values of the second field, whilst: M 12 , C 12 , K 12 , M 21 , C 21 , K 21 are the matrices of the field interactions [1, 6, 7].If there are no simplifications, the above mentioned global matrices are nonsymmetrical, which makes it difficult to directly solve the equation (