New developments on BE/BE multi-zone algorithms based on krylov solvers: applications to 3D frequency-dependent problems

TECHNICAL PAPERS

New developments on BE/BE multi-zone algorithms based on krylov solvers - applications to 3D frequency-dependent problems

F. C. de Araújo^I; W. J. Mansur^II; C. Dors^II; C. J. Martins^II

^IDepartamento de Eng. Civil/UFOP 35400-000 Ouro Preto, MG. Brazil [e-mail]

^IIPrograma de Eng. Civil/COPPE/UFRJ C. P. 68506 21945-970 Rio de Janeiro, RJ. Brazil [e-mail]

ABSTRACT

In this paper, new developments concerning the use of BE/BE coupling algorithms for solving 3D time-harmonic problems are reported. The algorithms are derived by considering different iterative solvers. Their chief idea is to work with the global sparse matrix of the coupled system, however without considering the many zero blocks associated with the non-coupled nodes of different subregions. The use of iterative solvers makes it possible to store and manipulate only the block matrices with non-zero coefficients. Preconditioned iterative solvers based on the Lanczos, bi-conjugate gradient, and GMRES (generalized minimal residual) methods are used to derive the BE/BE coupling algorithms. A description of the formulation of these solvers, which are completely general and can be applied to any non-singular, non-hermitian systems of equations, is provided. The performance of the algorithms is verified by solving some foundation-soil interaction problems. Important parameters for estimating the efficiency of the algorithms as required CPU times, matrix sparsity, and accuracy of the obtained responses are presented in the results of the paper.

Keywords: Boundary elements-BE, BE/BE coupling algorithms, frequency-dependent analysis, 3D problems, soil-foundations interaction, iterative solvers for complex-valued systems

Introduction

The Boundary Element Method (BEM) is an important numerical tool to perform frequency-domain dynamic analyses of engineering problems, specially those including unbounded domains. Some situations in which frequency-domain or stationary solutions in boundary element applications are necessary are those related to the estimation of resonance frequencies or when time-dependent fundamental solutions are not known, but the corresponding frequency-dependent ones are. For consideration of complex-valued material parameters as e.g. complex stiffness, used for modeling hysteretic damping (Clough and Penzien, 1993), frequency-domain analyses are also suitable. A general review of important papers up to year 1996 on the topic boundary element methods in dynamic analysis, including contributions on frequency- and time-domain analyses, can be found in Beskos (1987, 1997).

The present paper is concerned with the development of a general and efficient BE/BE coupling algorithm for analyzing time-harmonic problems. Only harmonic loads are considered here. Nevertheless, transient excitations can be also dealt with as a trivial extension of the harmonic case. As a boundary element formulation is involved, the procedure is very suitable for modeling exterior domain problems (infinite domains). On the other hand, the substructuring option (BE/BE coupling) is necessary for taking into account non-homogeneous domains, as e.g. soil layers.

In fact, substructuring strategies extend considerably the range of applications of boundary-integral-based methods. General considerations on subregions techniques can be found in known textbooks on boundary element methods (Banerjee, 1994; Kane, 1994), and in a series of papers published in the last two decades (Crotty, 1982; Kane, Kumar and Saigal, 1990; Rigby and Alliabadi, 1995; Bialecki et al., 1996; Ganguly et al., 1999). These works adopt either non-condensed or condensed strategies and are based on the use of direct solvers.

In previous papers by Araújo and Martins 2000, 2001, and Araújo, Martins and Mansur, 2001, an optimized way for developing generic BE/BE coupling algorithms based on the use of iterative solvers has been reported. In this paper, new advances concerning this BE/BE coupling strategy are presented. The algorithms are derived based on the consideration of different iterative solvers, and their main contribution is that the global matrix of the coupled system is not explicitly assembled; instead, algebraic subsystems corresponding to all its substructures are built and independently manipulated during the solution phase. At the end, the response of the coupled system is assembled. As iterative procedures are used for solving the global system of equations, the coefficients of each subsystem remain unchanged and no additional computer memory, beyond that already allocated for each subsystem, is needed. Non-preconditioned and Jacobi-preconditioned Krylov's solvers, namely the bi-conjugate gradient (Bi-CG), Lanczos, and GMRES (generalized minimal residual) schemes, are applied to derive the generic BE/BE coupling algorithms. A brief description of the formulation of these iterative solvers is provided. In Araújo and Martins (2001), Araújo, Martins and Mansur (2001), and Mansur, Araújo and Malaghini (1992), the complete derivation of the mentioned iterative procedures for real- and complex-valued systems of equations can be found.

The performance of the BE/BE coupling algorithms derived is observed by solving some 3D time-harmonic soil-foundation coupled problems. Important parameters for estimating the efficiency of the computer code, as required CPU times, sparsity of the global matrix, and response accuracy, are presented in the results of the paper.

The Direct Boundary Element Formulation for Frequency-Dependent Problems

The starting expression for frequency-dependent analysis by means of direct boundary element methods is the integral equation

where is the integral-equation jump term, equal to jump terms for elastostatics, and are the frequency-dependent fundamental solutions given by:

and

When the known complex boundary conditions

are considered, the boundary integral equation (1) gives the solution of the frequency-dependent problem (in terms of its complex amplitudes) at whole boundary.

Adopting usual discretization procedures, the boundary integral equation (1) is converted to the following frequency-dependent system of algebraic equations:

After introducing the boundary conditions shown in equations (4) and (5), equation (6) can be written as

where the complex vectors and contain the unknown and known boundary values respectively.

For determining the coefficients of the matrices H and G, a special integration process based on a triangular polar co-ordinate transformation, with optional use of integration subelements (Araújo, 1994), is employed. This integration scheme is very suitable for evaluating the weakly singular integrals involved. In case of Cauchy principal values, connected with the diagonal block matrices coefficients (DBMC) of the H matrix, the rigid body displacement criterion is considered. As known, to apply this criterion to BE models for semi-infinite layers, a special mesh of enclosing elements should be used (Araújo, Martins and Mansur, 2001).

Derivation of Krylov's Solvers for Non-Hermitian Systems

General Aspects

In the computational analysis of engineering systems, continuous models are converted into discrete systems. In this analysis phase, systems of algebraic equations, sometimes sparse, containing a number of equations varying from a few hundreds to a few millions may be originated. Iterative procedures may be very efficient to solve such systems, because, besides preserving the matrix sparsity, they can substantially reduce the CPU time in case of large systems, usual in realistic problems. In Hageman and Young (1981) and Hackbusch (1991), a comprehensive study on general aspects of iterative procedures can be found.

Among the iterative schemes, those classified as Krylov's solvers or polynomial acceleration procedures (Hageman and Young, 1981; Araújo and Martins, 2001) have specially attracted the attention of engineers. The reasons for that is the exceptionally good convergence properties of their preconditioned versions. Important results for non-hermitian BE systems are reported by Araújo and Martins (2001), Mansur, Araújo and Malaghini (1992), Barra et al. (1993), Prasad et al. (1994), and Valente and Pina (2001). In Saad and van der Vorst (2000), a complete review on iterative solvers is provided; more than 200 references are cited.

In this paper, the Lanczos, the biconjugate gradient, and the GMRES algorithms are applied to solve complex-valued, non-hermitian systems of algebraic equations arising in frequency-domain formulations.

The Lancos and Bi-CG Scheme

The starting point for deriving the Lanczos and the bi-conjugate gradient algorithm (Bi-CG) is the Lanczos tridiagonalization algorithm (Araújo, 1989; Wilkinson, 1965). With aid of this algorithm the Lanczos solver is directly derived, and starting from the latter, the Bi-CG scheme can be obtained. Below, a short description of the derivation of their formulation is discussed.

The Lanczos tridiagonalization algorithm can be applied, in the same way as for real matrices, to derive from a matrix and its transpose, , two sequences of complex-valued vectors, respectively, given by

where , are known initial guess vectors. These sequences of vectors are mutually orthogonal to each other, i.e., , and thus, linearly independent for , N being the dimension of the complex space in question. Thus, the following property is verified:

The property expressed by Eq. (10) is a very important one, since it is concerned with the finite termination of the iterative schemes (Hageman and Young, 1981; Araújo, 1989). Moreover, if an usual inner product is used at the orthogonalization of vectors , that is, if for instance a non-hermitian one is considered, expressions for parameters in Eqs. (8) and (9) similar to those for real matrices are obtained. It results (Araújo, 1989; Wilkinson, 1965):

By adopting a three-term iterative formula,

n denoting the iteration order, and forcing that the associated residual,

be a Lanczos vector given by Eqs. (8), the Lanczos solver can be derived. Of course, an auxiliary system of equations associated with must be taken into account; its residual,

is forced to be the Lanzos vector given by (9). The parameters in expressions (14)-(16) are then determined by comparing Eqs. (15) and (16) with Eqs. (8) and (9) respectively. It results:

with .

By adopting a two-term recursive formula, the Bi-CG scheme is derived. One obtains:

The iterative Lanczos and Bi-CG schemes shown above possess the property of being convergent at a maximum of n iterations, where . Nevertheless, as a consequence of truncating errors introduced during the data processing, convergence may not be achieved for . In order to avoid non-convergence cases and also to accelerate the respective convergence rates, preconditioning strategies have been taken into account (Araújo, 1989; Mansur, Araújo and Malaghini, 1992; Barra et al., 1993; Prasad et al., 1994). In this work, the only preconditioning matrix Q considered is the Jacobi one, defined by the (real- or complex-valued) diagonal of the respective system matrices.

The GMRES Scheme

Let be the k-dimensional Krylov's space associated with and , the complex-valued initial residual vector and the system matrix respectively (Hageman and Young, 1981; Araújo, 1989). By taking an orthonormalized basis for this subspace, with given by

a Krylov's iterative scheme can be generated by imposing the condition:

with (Hageman and Young, 1981; Araújo, 1989; Araújo and Martins, 2001), where the parameters are determined according to a certain criterion to be established, and in Eq. (27) is a Krylov's vector (see Fig. 1). Concerning the GMRES solver, is obtained so that be an optimum approximation for , the exact solution of the algebraic system (see Fig. 1). In other words, the following minimization problem is established: find so that

A way to solve this problem is by imposing the condition that be orthogonal to the Krylov's subspace , i.e., by considering

By using the Gram-Schmidt orthogonalization procedure to generate the orthonormalized basis , the following solution for the optimization problem above is obtained:

where is the Hessenberg matrix

With

GMRES solvers also possess a finite termination property, and the same comments previously done regarding the importance of preconditioning remain valid here. The Jacobi preconditioner is also the only one used in this paper for accelerating the GMRES solver. In order to reduce the number of Krylov's vectors used for expressing the system response, restarting option is considered; namely blocks with a fixed size of Krylov's vectors, N being the system order, are adopted.

Formulation of the Multi-Zone BE Algorithm

Related to a certain subregion k, the system of equations (7) can be written as

where

are the volume forces, and the indexes above are defined as follows: ns is the number of subregions, , the outer boundary part of the k-th subregion with unknown displacements, , that one with unknown forces, and , its interface (with unknown displacements and forces). Notice that .

As corresponding to the interface nodes of a given subregion k, both displacements and tractions are unknown, an additional number of equations will be necessary for calculating all its unknown boundary values (including, of course, its interface values). In fact the compatibility and equilibrium equations between common subregion interfaces provide the extra equations which permit assembling together the subsystems corresponding to all substructures of the problem. The number of equations and unknowns, after introducing the boundary conditions, become the same and a unique solution to the problem is obtained.

In fact, an interface is defined by a set of coupled nodes, which in turn are defined according to the following criterion: coupled nodes pertain to two different domains, have the same co-ordinates, and the corresponding unit normal outward vectors at them are opposite to each other. For the generic coupled subdomains shown in Fig. 2, the interfaces between the various subdomains are indicated by a global numbering. For instance, interface is that separating subdomains and , interface , that between subdomains and , and so on. Explicitly, the global coupled system in its most general form is given by

where s is the total number of existing interfaces, and , the coupling submatrix (containing terms of and ) associated with the m-th (global) interface. Of course, if the m-th interface does not pertain to the k-th subregion. Regarding the BE model in Fig. 2, and for instance are nonzero, as only the interfaces and pertain to subdomain .

As continuous boundary elements are considered in the code, special attention must be paid to the simulation of edges and corners. In such cases, one must deal correctly with common interface nodes, defined as nodes common to more than one interface (see Fig. 2). Thus, in order to obtain the same number of equations and unknowns, traction continuity conditions, which make it possible to reduce the number of unknown interface tractions, must be taken into account. In case of smooth interface parts, the respective mesh nodes pertain to an only interface; non-common interface nodes (see Fig. 2) are then present, and equal number of equations and unknowns are obtained only with the consideration of coupling conditions. Notice that the interfaces in Fig. 2 indicated with dot lines (along the first soil layer) are used to simulate corners and edges of the foundations. They are necessary because no discontinuous boundary elements are used.

The continuity and traction conditions mentioned above may be defined as follows: coupling conditions are associated with nodes of the same geometrical position and opposite outward normal vectors; traction continuity conditions, with nodes of same geometrical position and same outward normal vectors. In Fig. 3, these conditions are shown for a generic vector quantity f at nodes i and j. In the computer code, they are automatically introduced into the coupled model.

First, a coupling examination is carried out so as to identify the coupled nodes and to generate the variables indicating which nodes are coupled with which nodes, and, second, the coupled nodes at which additional traction continuity conditions must be considered are determined and the corresponding variables, generated.

If only non-common interface nodes are present at the model, no traction continuity condition is necessary. Otherwise, that is, in case of common interface nodes, the number of necessary traction continuities is calculated by the formula:

Where

In expressions (37)-(39), is the number of variables per common node and , the corresponding number of equations. nnode and nsnode denote, respectively, the number of nodes and subregions involved at the common interface node. In fact, a common interface node is a nest of coupled multiple nodes. For instance, for the common interface node zoomed in at Figure 2, one obtains: , and . As from the coupling conditions, , then to calculate , it is necessary to introduce one continuity condition, for example, .

Parallel to the ideas above, directly concerned with the description of the coupling process, another strategy is used to completely avoid operating with the many zero blocks, necessarily present at the system (36). In fact, the basic idea to do this is to use iterative solvers, which fundamentally work with no matrix transformation. Thus, the system (36) must not be explicitly assembled. Instead, a loop over the subregions is carried out and the independent contributions of each subregion are taken into account by the operation under consideration (matrix-vector or vector-vector multiplications). At end, the final result of the respective operation is gathered. Notice that no zero block is stored and manipulated during the solution of the system.

In this paper, eight different solvers are actually used in conjunction with the coupling technique; namely, the Lanczos, J-Lanczos, Bi-CG, J-Bi-CG, GMRES, J-GMRES, J-Lanczos(-eq), and J-Bi-CG(-eq). The general description of these methods (Lanczos, Bi-CG, GMRES) is given in the section above on Krylov's solvers for non-hermitian systems. The letter "J" preceding the name of the solvers stands for Jacobi-preconditioning, and the particle "eq", designates that the real equivalent version of the respective solver is used (Araújo, Martins and Mansur, 2001). If the particle "eq" is omitted, it means that a complex solver is applied.

Applications

To observe the performance of the coupling algorithms, some problems concerning the dynamic interaction between soil and footings are analyzed (see Fig. 4). As a matter of fact, the following problems are considered: cases of one and two footings, embedded or non-embedded, rigid or flexible, coupled with the soil. In the discretization models adopted, each footing and the soil are considered as one subregion.

Case Study 1: Single Non-Embedded Footing Coupled With the Soil

To validate the computational code developed, a series of analyses in the frequency-domain of a system consisting of a single footing coupled with the soil is carried out. The influence of the footing rigidity on the dynamic response of the system is observed. The relative stiffness of the footing, defined by

where and are the elasticity modulus of the footing and of the soil, respectively, and E is used to measure the relative footing stiffness. The dimensionless thickness is taken to be 0.10, and the footing is considered to be rigid for . The BE model adopted to solve this problem is shown in Fig. 5.

The footing is discretized with 160 eight-noded boundary elements (BE), which correspond to a mesh of 546 nodes. The BE mesh for the soil contains 320 eight-noded boundary elements and 1057 nodes, and a mesh of 48 enclosing elements (EE), with 161 nodes, is additionally used. The concentration of elements along the sides of the footing is necessary to simulate the behavior of the interface tractions as the relative rigidity increases.

This problem was also analyzed by Wittaker and Christiano, 1982, and Qian et. al., 1996. Their responses are plotted together with those obtained with the procedure proposed in this paper in the graphics of Fig. 6. In these graphics, the dimensionless vertical displacement amplitude, D , defined through

where q is the uniformly distributed load and the vertical displacement amplitude D , is plotted versus the dimensionless frequency,

with c_s being the shear wave velocity of the soil and a half of the foundation width.

Case Study 2: One Non-Embedded Footing (NEF) Coupled With the Soil

In this case, which aim to show the computational efficiency of the algorithm, a harmonic distributed load of amplitude and frequency , acting in the vertical direction at the top surface of the footing is considered. The physical constants for the footing and the soil are assumed to be:

Soil parameters:

hysteretic damping ;

Footing: ,

.

Quadrangular eight-noded boundary elements are used in all problems. The thickness of the footings is h = 0.19 m and their half side length, a = 0.76 m , so that a₀ = 0.3746.

The model shown in Fig. 7 is considered. The footing (Fig. 9(b)) is discretized with 128 boundary elements (450 nodes), and the soil (Fig. 9(a)), with 144 boundary elements (529 nodes). The enclosing element mesh (EE mesh) contains 48 enclosing elements (161 nodes).

Case Study 3: One Embedded Footing (EF) Coupled with the Soil

A BE model similar to that shown in Fig. 7 is adopted. The same loading, geometric and physical characteristics used in case study 2 are also considered here. The BE mesh for the footing contains 160 boundary elements (546 nodes) and for the soil, 176 boundary elements (625 nodes). The EE mesh has 48 enclosing elements and 161 nodes.

Case Study 4: Two Non-Embedded Footing (2NEF) Coupled with the Soil

The BE model shown in Fig. 8 is used. Again the same characteristics of case study 2 are adopted. Each footing (Fig. 9(b)) is discretized with 128 boundary elements (450 nodes), and the soil (Fig. 9(c)), with 288 boundary elements (1041 nodes). The respective EE mesh has 104 enclosing elements (333 nodes), and the distributed harmonic load acts on both footings.

The response of the soil-footing(s) interaction problems for the case studies 2-4 described above is shown by the graphics of Figures 10-13 in terms of the vertical displacement amplitudes according to the relations E_f/E_s. The sparsity of the coupled systems (defined here as the ratio between the number of zeroes and the total number of coefficients) in the cases 2, 3 and 4 is 31.6%, 26.5%, and 46.5% respectively. The number of iterations (scaled by the system order) and the CPU times (scaled by the corresponding CPU time of the standard Gauss solver with columns and rows pivoting option) attained in the analyses are plotted against the relationship E_f/E_s in the graphics 14, 15, 17 and 18.

 

Finally, the bar graphics in Figure 16 and ¹⁹ show the CPU time distribution for the different stages of each analysis. This is suitable to emphasize the importance of the separated phases of the whole analysis.

The iterative scheme was stopped whenever being the residual vector at the current iteration. The computer used for carrying out the analyses had an INTEL processor with 1.0GHz and 768Mbytes random access memory, and the code was developed in a COMPAQ VISUAL FORTRAN environment, version 6.5.0.

Conclusions

A general procedure is presented to calculate the dynamic response of 3D foundations, rigid or flexible, of any geometrical shape and under any spatial arrangement, coupled with the soil.

The results presented in Figures 10-13 could not be compared with ones furnished by other authors in all cases analyzed, as it was not possible to find all respective results in the technical literature consulted. Nevertheless, the good agreement between the responses obtained with the method proposed here and those furnished by Wittaker and Christiano, 1982, and Qian et. al., 1996 (Fig. 6), ensure the quality of the other responses first published in this paper. This case study was used to validate the procedure implemented.

Concerning the performance of the substructuring algorithm, one sees that convergence was reached at low rates of CPU times, if compared with those attained by using a direct standard Gauss solver (Figs. 15 and ¹⁸). Additionally, Figs. 14, 15, ¹⁷ and ¹⁸ given an interesting insight on how the performance of the solver varies according to the relationship E_f/E_s. On the other hand, the coupling procedure proposed is very suitable for developing parallelized BE computer codes. Thus, the procedure might constitute a very promising technique for treating realistic 3D engineering problems via Boundary Element Methods, usually represented by means of algebraic systems of equations with hundreds of thousands of equations. In these cases the factorization or inversion of the resulting system matrix is a very hard computational task, so that the proposed strategy might be very attractive.

Acknowledgements

The first author gratefully acknowledges the financial support from the Brazilian Research Council - CNPq and from the Research Foundation of Minas Gerais- FAPEMIG

References

Araújo, F. C., 1994, "Zeitbereichslösung Linearer Dreidimensionaler Probleme der Elastodynamik mit einer gekoppelten BE/FE-Methode" (in german), Ph.D Thesis, Technische Universität Braunschweig, Germany, 182p.

Araújo, F. C., 1989, "Iterative Techniques for Solving Linear Systems of Equations Originated from the Boundary Element Method" (in Portuguese), M.Sc. Thesis, COPPE - Federal University of Rio de Janeiro, Brazil, 116p.

Araújo, F. C., Martins, C. J. and W. J. Mansur, 2001, "An Efficient BE Iterative-Solver-Based Substructuring Algorithm for 3D Time-Harmonic Problems in Elastodynamics", Engineering Analysis with Boundary Elements, Vol. 25, pp. 795-803.

Araújo, F. C. and Martins, C. J., 2001, "A Study of Efficiency of Multizone BE/BE Coupling Algorithms Based on Iterative Solvers - Applications to 3D Time-Harmonic Problems", In M. Denda, M., H. Aliabadi and A. Charafi (eds.), Proc. Joint Meeting of Boundary Element Techniques and International Association for Boundary Elements, New Brunswick-NJ-USA, 16-18 July, Advances in Boundary Element Techniques II, Geneva: Hoggar, Vol.1. pp.21-29.

Araújo, F. C. and Martins, C. J., 2000, "Efficient Generalized Boundary Element Substructure Algorithm for Analysis of 3D Machinery Foundations", In J. Tassoulas (ed.), Proc. Fourteenth Engineering Mechanics Conference, ASCE, University of Texas at Austin - USA, 21 - 24 May.

Banerjee, P. K., 1994, "The Boundary Element Method in Engineering", McGraw-Hill Book Company, New York-USA, 512p.

Barra, L. P. S., Coutinho, A. L. G. A., Telles, J. C. F. and Mansur, W. J., 1993, "Multi-Level Hierarchical Preconditioners for Boundary Element Systems", Engng. Anal. Boundary Elements, Vol. 12, pp. 103-109.

Beskos, D. E., 1987, "Boundary Element Methods in Dynamic Analysis", Appl. Mechanics Reviews, Vol. 40, pp. 1-23.

Beskos, D. E., 1997, "Boundary Element Methods in Dynamic Analysis: Part II (1986-1996)", Appl. Mechanics Reviews, Vol. 50, pp. 149-197.

Bialecki, R. A., Merkel, M., Mews, H. and Kuhn, G., 1996, "In- and Out-of-Core BEM Equation Solver with Parallel and Non-Linear Options", Int. J. Num. Methods in Engineering, Vol. 39, 4215-4242.

Clough, R. W. and Penzien J., 1993, "Dynamics of Structures", McGraw-Hill, Inc., second edition, New York-USA, 739p.

Crotty, J. M., 1982, "A Block Equation Solver for Large Unsymmetric Matrices Arising in the Boundary Element Method", Int. J. Num. Methods in Engineering, Vol. 18, pp. 997-1017.

Ganguly, S., Layton, J. B., Balakrishna, C. and Kane, J. H., 1999, "A Fully Symmetric Multi-Zone Galerkin Boundary Element Method", Int. J. Num. Methods in Engineering, Vol. 44, 991-1009.

Hackbusch, W. 1991, "Iterative Lösung Grosser Schwachbesetzter Gleichungssysteme", B. G. Teubner Stuttgart, Germany, 382p.

Hageman, L. A. and Young, D. M., 1981, "Applied Iterative Methods", Academic Press, Inc., 386pp.

Kane, J. H., Kumar, B. L. and Saigal, S., 1990, "An Arbitrary Condensing, Noncondensing Solution Strategy for Large Scale, Multi-Zone Boundary Element Analysis", Comput. Meth. Appl. Mech. Engng, Vol. 79, 219-244.

Kane, J. H., 1994, "Boundary Element Analysis in Engineering Continuum Mechanics", Prentice-Hall, Englewood Cliffs, NJ-USA.

Mansur, W. J., Araújo, F. C. & Malaghini J. E. B., 1992, "Solution of BEM Systems of Equations via Iterative Techniques", Int. J. Num. Methods in Engineering, Vol. 33, pp. 1823-1841.

Prasad, K. G., Kane, J. H., Keyes, D. E. and Balakrishna, C., 1994, "Preconditioned Krylov Solvers for BEA", Int. J. Num. Methods in Engineering, Vol. 37, pp. 1651-1672.

Qian, J., Tham, L. G. and Cheung, Y. K., 1996, "Dynamic Cross-Interaction between Flexible Surface Footings by Combined BEM and FEM", Earthquake Engineering and Structural Dynamics, Vol. 25, pp. 509-526.

Rigby, R. H. and Alliabadi, M. H., 1995, "Out-of-Core Solver for Large, Multi-Zone Boundary Element Matrices", Int. J. Num. Methods in Engineering, Vol. 38, 1507-1533.

Saad, Y and van der Vorst, H. A., 2000, "Iterative Solution of Linear Systems in the 20th Century", Journal of Computational and Applied Mathematics, Vol. 123, pp. 1-33.

Valente, F. P. and Pina, H. L., 2001, "Iterative Techniques for 3-D Boundary Element Method Systems of Equations", Engineering Analysis with Boundary Elements, Vol. 25, pp. 423-429.

Wilkinson, J. H., 1965, "The Algebraic Eigenvalue Problem", Claredon Press, Oxford.

Whittaker, W. L. and Christiano, P., 1982, "Dynamic Response of Plate on Elastic Half-Space", J. Eng. Mech. Div. ASCE, Vol. 108, pp. 133-154.

Paper accepted January, 2004. Technical Editor: José Roberto de F. Arruda

Publication Dates

• Publication in this collection
  12 Aug 2004
• Date of issue
  June 2004

History

• Accepted
  Jan 2004