An Adaptive Multigrid Solver for Plate Vibration and Buckling Problems

PERGAMON

Computers and Structures 69 (1998) 625±637

An adaptive multigrid solver for plate vibration and buckling problems Salvatore Lopez a, Stefania Fortino b, Ra€aele Casciaro b a

Istituto di Scienza delle Costruzioni, UniversitaÁ di Bologna, 40136, Bologna, Italy b Dipartimento di Strutture, UniversitaÁ della Calabria, 87030, Rende (CS), Italy Received 9 January 1997; accepted 5 February 1998

Abstract The paper describes a ®nite element multigrid strategy for the solution of plate vibration and buckling problems. The solution procedure combines the residual iteration scheme used in Casciaro R, Aristodemo M. International Conference on Finite Elements Nonlinear Solid and Structural Mechanics, Gelio, Norway, 1977, (main loop) with the adaptive multigrid scheme proposed in Lopez S, Casciaro R. International Journal of Numerical methods in Engineering 1997;40:919±96. A subspace iteration technique allows several principal eigensolutions to be obtained simultaneously. The sequence of meshes used by the multigrid process is generated through an adaptive local re®nement based on the Zienkiewicz±Zhu estimate of the discretization error Zienkiewicz OC, Zhu JZ. International Journal of Numerical Methods in Engineering 1987;24:337±57. A pointer-based data structure and an HC Aristoderno M. Computers and Structures 1985;21(5):987±93, ®nite element discretization make the solution process highly ecient. Some numerical tests show the e€ectiveness of the algorithm. # 1998 Elsevier Science Ltd. All rights reserved. Keywords: Eigenvalue analysis; Multigrid methods; Local re®nement

1. Introduction Vibration and buckling problems are usually reduced, through a ®nite element discretization, to a generalized matrix eigenvalue equation of the type: ‰A ‡ lBŠu ˆ 0; where A and B are [n  n] symmetric matrices, the ®rst being positive de®nite. In vibration problems, the equation is written Ku ÿ lMu = 0, where K and M are the sti€ness and mass matrix, respectively, u is the natural vibration mode vector and l: = o2 is the square of the natural frequency of the system. In buckling problems, the equation is written Keu + lKgu = 0 where Ke and Kg are the so-called elastic and geometrical sti€ness matrices, respectively, u is the buckling mode vector and l is the load factor. In both cases, only a few p principal modes, i.e. pWn, are required.

Finite element analysis of these problems requires ®ne discretization grids to obtain accurate numerical solutions. For buckling problems, in particular, very ®ne ®nite element meshes are needed for an adequate representation of the post-buckling behaviour [1±3]. Furthermore, in the presence of singularities, a localized re®nement may be convenient in order to improve the solution. Standard methods for solving F.E. eigenvalue problems based on matrix manipulations are computationally expensive for large dimension problems with very ®ne discretization grids. Multigrid strategies represent a convenient alternative in this case because of their ability to solve problems with n unknowns in O(n) arithmetical operations [4, 5]. This goal is achieved by analyzing a (theoretically in®nite) sequence of discretization meshes of increasing density and by obtaining the solution on the ®ner grids through iterative corrections performed, essentially, on the coarser ones. The

0045-7949/98/$ - see front matter # 1998 Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 4 9 ( 9 7 ) 0 0 0 9 0 - X

626

S. Lopez et al. / Computers and Structures 69 (1998) 625±637

strategy fully exploits the smoothing properties of the relaxation scheme and the transfer of information between meshes of di€erent size to achieve optimal computational eciency with low storage demand [6± 8]. Even greater eciency can be achieved by combining multigrid solutions with adaptive re®nement techniques, designed to by-pass the domain zones characterized by a suciently accurate coarse solution. The solution strategy presented in this paper combines the adaptive multigrid solver described in [9, 10] with the residual iteration scheme used in [11] for the analysis of eigenvalue problems coming from the bifurcation theory. The iterative nature and the similarity of the multigrid (inner) loop and the residual iteration (main) scheme allows their ecient combination. Further features of the procedure are: . A subspace orthogonalization technique is used to obtain several eigensolutions simultaneously. . The re®nement is based on standard bisection but covers only the progressively smaller subdomains where higher resolution is needed according to a Zienkiewicz±Zhu [12, 13] local error indicator (adaptive local re®nement). . All quantities relevant to a node of the grid are collected in a single structure which uses pointers for obtaining adjacencies and inter-grid links. The use of dynamic lists allows a simple treatment of unstructured and locally re®ned grids with a marginal storage demand. The paper is organized as follows. Section 2 describes the main aspects of the solution strategy. Section 3 discusses some implementation details of the local re®nement process, including the data structure used and the error indicator. Section 4 reports some numerical results related to plate vibration and buckling problems. Some ®nal comments are given in the conclusions.

2. A Multigrid Strategy for Eigenvalue Problems We consider a set of grids {Gh, h = 0


m}, with decreasing mesh-size, approximating the same domain O where the self-adjoint di€erential generalized eigenvalue problem: Au ‡ lBu ˆ 0; u 6ˆ 0

…1†

is assigned. A ®nite element discrete ®eld uh[x] $ Uh being de®ned on each grid, the discrete solution in Gh is given by the condition Kh ‰lh Šuh :ˆ …Ah ‡ lh Bh †uh ˆ 0; uh 6ˆ 0;

…2†

where uh and lh are the Gh-eigenvectors and eigenvalues, respectively.

The eigenvalue Eq. (2) can be solved using the following iteration scheme (main loop): ( …k‡1† …k† …k† ˆ u…k† uh h ÿ Hh Kh ‰lh Šuh ; …3† l…k‡1† ˆ 1=kuh…k‡1† k h where the matrix Hh>0 represents an approximation of the inverse of Kh[0]. The scheme, which has been extensively used for solving nonlinear buckling problems (see [1, 2]), converges to the solution with smallest positive lh under weak hypotheses. When Kh[lh] is linear in lh, as in the present case, it corresponds to a rewriting of the inverse power method of which it retains excellent convergence behaviour, even allowing quite rough approximations of the inverse. The latter feature can be exploited by implicitly de®ning this approximation by an iterative multigrid rough solution of the problem: Ah uh…k† ˆ fh…k† ;

…4†

where uh…k† :ˆ u…k‡1† ÿ u…k† h h ; fh…k† :ˆ ÿK‰lh…k† Šu…k† h ; lh…k† :ˆ 1=kuh…k† k:

…5†

That is, the solution of the inner problem (4) is obtained by means of multigrid cycles (inner loop) where the forcing term fh is periodically updated from the Eq. (5) (main loop). The analysis is extended to the simultaneous search for the ®rst p eigensolutions by introducing a subspace orthogonalization phase in the iterative process.The overall solution process, shown in Fig. 1 is de®ned as follows: 2.1. Main loop 1. Starting from the ®rst grid G0, the process is initialized by performing a standard matrix-based eigenvalue analysis which provides all the p required eigenvectors (due to the small number of unknowns, this phase is computationally inexpensive). The reference values: Fi ˆ kK0 ‰0Šui k; i ˆ 1


p

…6†

are computed and the eigenvectors are transferred by interpolation to the grid Gm, m = 1. The eigenvector index i is initialized to 0. 2. The index i is updated to i + 1, then the current (k) primary residual fm is computed from Eq. (5): …k† …k† fm ˆ ÿKm ‰l…k† m;i Šum;i :

…7†

Note that this operation does not require the actual construction of matrix Km[l(k) m,i] but is performed

S. Lopez et al. / Computers and Structures 69 (1998) 625±637

627

Fig. 1. Pseudo-code list of the solution process.

through direct assembling of the vector fm using its element contributions. 3. The multigrid inner loop is then performed providing the correction u(k) h as described in the sequel. 4. The ith eigensolution is improved according to Eq. (5): …k‡1† …k† um;i ˆ u…k† m;i ‡ u m;i ;

…k‡1† li…k‡1† ˆ 1=kum;i k:

…8†

5. The steps from (1) to (4) are cycled until all the eigenvectors have been treated. Then: (a) If the solution error k f(k) m k/Fi is less than the assigned tolerance Z for all eigenvectors, then a new grid is generated covering domain zones that need a ®ner discretization, if these exist (see Section 3). Then m is updated to m + 1 and the solution is transferred to the new grid by interpolation. (b) Otherwise, a subspace orthogonalization is performed as described in the sequel, i is reset to 0 and the process restarts from step (1). 6. The steps from (1) to (5) are repeated until there are no further zones to be re®ned. Then the process is stopped. The multigrid inner loop and the subspace orthogonalization phase [steps (3) and (5b) of the main process] are detailed in the following: 2.2. Multigrid inner loop (step 3) The multigrid loop is performed as detailed in [10], i.e. c1<<1, c2<1 and c3<1 being appropriate factors, let

:ˆ f…k† …k;j† r…k;j† h h ÿ Ah u h

…9†

be the residual at the jth relaxation sweep on the current Gh grid, and let Zm: = c1 kf(k) h k the relaxation tolerance. Successive relaxation cycles (j = 1, 2


) are performed on Gh until one of the following conditions prevails: (a) the relaxation convergence deteriorates due to the presence of prevailing low frequency components, (k,j ÿ 1) k>c2 ; then the residual is transi.e. kr(k,j) h k/krh k ferred to the coarser grid Gh ÿ 1, Zh ÿ 1: = c3kr(k,j) h is set and the process continues on the coarse grid; (b) h < m and the current residual is suciently small, i.e. kr(k,j) h k < Zh; then the correction obtained is transferred to the next ®nest grid Gh + 1 where the iterative process continues; k is less than (c) h = m and the residual error k r(k,j) h Zh; then the inner process is stopped and the obtained correction u(k) is provided to the main h loop.

2.3. Subspace orthogonalization phase (step 5b) Let uh,i, i = 1


p be the current estimate for the ith eigenvector, the following matrices are computed:  M0 :ˆ fm0ij g; m0ij ˆ uTi Auj ; i; j ˆ 1::p …10† M1 :ˆ fm1ij g; m1ij ˆ uTi Buj ; i; j ˆ 1::p: A standard matrix solution for the reduced p-dimensional eigenvalue problem

628

…M0 ‡ lM1 †x ˆ 0

S. Lopez et al. / Computers and Structures 69 (1998) 625±637

…11†

is then performed, whose solution xi, i = 1


p is transferred in full space by linear combination: ui 3‰u1


up Šxi :

…12†

Note that the link between the main and the inner iterations is only represented by the forcing term f(k) m,i. Apart from the updating at each main step of this term, a standard multigrid strategy for linear problems is used. An adaptive local re®nement strategy can be introduced because of a local evaluation of the discretization error, according to Zienkiewicz±Zhu [12], is available before the generation of the new grid. The constants c1, c2 and c3 which control the inner loop noticeably a€ect the overall eciency of the process. In particular, the value of c1 in¯uences the main loop behaviour: smaller values of c1 produce more inner cycles, but less main cycles, and vice versa. Values in the range 10ÿ2±10ÿ3 may represent a good compromise. The constants c2 and c3 a€ect the intergrid transfer conditions. They can be precisely tuned by reference to the multigrid operators used in order to optimize the solution process (see [4, 6]). The values c2=0.60 and c3=0.15 have been used in the applications. 3. Implementation Details 3.1. Local re®nement The discretization is based on rectangular elements and uses a standard bisection as re®nement rule. This provides uniform ®ne grids and allows the use of

recurrent, programmed once and for all, procedures for both relaxation and inter-grid transfer operators. An appreciable simpli®cation in code and in data structure design is thus achieved. Non-uniform discretization is obtained by con®ning the re®ned grid to progressively smaller subdomains covering only the zones which actually need re®nement according to a local discretization error indicator. This retains the local nature of the algorithm with a considerable saving in storage and algebraic operations. Note that, during the re®nement process, the overall discretization grid Mm is composed of several local grids Ghm with di€erent levels of re®nement collecting the ®nest elements present in the zones of the domain. A further re®nement step involves all grids of Mm, i.e. the ®nest grid Gm and the coarse grids Ghm are all potentially subject to re®nement. In this way, a constant level of accuracy is assured [9, 14]. The new ®ne grid contains all nodes appertaining to the domain parts to be re®ned (nodes denoted as F in Fig. 2) and also includes a boundary band which collects all nodes connected to the previous ones through the relaxation or transfer operators (nodes R and T in Fig. 2). The boundary band dimension obviously depends on the ®nite elements used. The relaxation process is performed on the F nodes, and also involves R nodes. The interpolation and the restriction processes involve F, R and T nodes. The transfer between grids is realised by means of canonical transfer operators and R and T nodes are initialized by direct interpolation from the coarse solution (C nodes) and unmodi®ed during the relaxation process. These features allow both equilibrium and continuity inter-grid conditions to be satis®ed.

Fig. 2. h ÿ 1 and h grids: node classi®cation.

S. Lopez et al. / Computers and Structures 69 (1998) 625±637

It is worth mentioning that a similar approach is followed in [15]. However, in that paper the T nodes are not explicitly introduced and the values of R nodes are computed directly from the coarse solution. In this way, a little storage is saved, but at the cost of either computational time or coding. 3.2. Data structure During the multigrid analysis, the generic grid Gh, while composed of several parts Ghm (m>h), created by the local re®nement process in successive steps, has a recursive structure. In fact, each node is linked with the adjacent ones (North, South, East and West) on the same grid and to the correspondent `®ne' node in the next ®nest grid. The C++ language is very suitable for handling such a situation. In particular, data can be organized through two basic structures: Node and Grid. The ®rst one collects all the quantities de®ned on the single node of the grid, such as the boundary conditions (ct), the nodal coordinates (x and y), the forcing vector (f), the eigenmodes (uk) and its iterative corrections (du). Five Node-pointers are also included to express the nodal adjacencies on the same grid (the horizontal connectors (N, S, E and W) ) and the link with the ®ner grid (the vertical connector D), as shown in Fig. 3. Note that vectors fh and rh used by the process are not required at the same time, so they can all be stored in the same memory location (f). All nodes of the grid are collected in a multiple-segments array which is constructed in successive re®nement steps and is accessed by the pointer to the ®rst element of each segment. The Grid structure collects all quantities relevant to the overall grid, such as the number of nodes (nn) contained in each segment, the mesh-size (MeshSize), the tolerances (Tol relax) and (Tol re®n) used by the relaxation and the re®nement procedures, the current values of the residual and the discretization errors (Err relax) and (Err re®n) and the pointers (nodes) to

629

the `®rst' Node's of the grid. All Grid's used by the analysis are collected in a single array which is accessed by the pointer to its ®rst element. A special Node called Node0 is used for the description of the boundaries: a node is on the boundary if one or more of its horizontal pointers points to Node0. All the pointers in Node0 point to Node0 itself and ct is set to be constrained. In this way the multigrid procedures that work with recognition on the adjacent nodes, remain unchanged even when the process runs to boundary nodes. The node Node0, moreover, is used to identify the subdomains which are not to be re®ned. All nodes found in these domains have D pointed to Node0 and the whole grid Mm is then indenti®ed. Note that the data structure described above, is quite similar, in principle, to that proposed by Rivara (see [15, 16]) while exploiting the features of the Clanguage. The structure Node represents, in fact, a direct implementation of the concept of `molecula' presented in that paper and the complex by-index relations required by its FORTRAN-based programming are completely avoided by the use of pointers. In particular, the introduction of pointer D creates a direct correspondence between the nodes of di€erent grids. In such way a binary search is not necessary. It must be noted, ®nally, that the memory required to store the adjacency pointers is quite small in comparison to that required by the ¯oating-point variables.

3.3. Re®nement procedure and error indicator The solution having already been obtained on grid Mm, the zones which need further re®ning are identi®ed using the a-posteriori estimate of discretization error proposed by Zienkiewicz and Zhu in [12]. This estimate is based on the assumption that a smooth and good enough approximation of the `exact' stress ®eld s*[x] is provided by its nodal value s*m through the same interpolation used for the displacement ®eld. We obtain: s*‰xŠ1Nm ‰xŠsm Z O

Qm ‰xŠT …s  ‰xŠ ÿ sm ‰xŠ†dO ˆ 0

Nm[x] being the displacement shape functions: um ‰xŠ ˆ Nm ‰xŠum ;

Fig. 3. Horizontal and vertical connectivity.

where um is the nodal displacement vector and Ql(x) a suitable projection operator. Using the energy norm for evaluating both the error and the solution:

630

kem ‰xŠk2 :ˆ

S. Lopez et al. / Computers and Structures 69 (1998) 625±637

Z O

fsm ‰xŠ ÿ s*‰xŠgT Dÿ1 fsm ‰xŠ ÿ s*‰xŠgdO …13†

kum ‰xŠk2 :ˆ

Z O

s*‰xŠTDÿ1 s*‰xŠdO

…14†

and assuming Dirac's delta function [13] for the operator Qm(x), we derive the following error test:  1=2 Oi kem ‰xŠki < E ; …15† …kum ‰xŠk2 ‡ kem ‰xŠk2 † O where E is the required relative tolerance on the error and kem(x)ki is obtained from (13) by integrating on the ith element with area Oi. The re®nement process is performed by scanning all elements of grid Gm in sequence. Elements where the test (15) fails are re®ned. During this process `®ne' nodes are generated, the `coarse' D pointer is properly set and new `®ne' nodal values are initialized by interpolation. Note that the resulting list of nodes consists, in general, of some contiguous segments. The relative multicolour ordering of nodes proves to be convenient within a Gauss±Seidel relaxation scheme (see [5, 10]). The absence of elements to be re®ned indicates that the required solution has been achieved and then the solution process can stop.

1 1 1 f1 …x† ˆ ÿ x ‡ x2 ; 8 2 2 3 2 f2 …x† ˆ ÿ x ; 4 1 1 1 f3 …x† ˆ ‡ x ‡ x2 : 8 2 2

Note that (17) represent the shape functions of the one-dimensional quadratic interpolation u(x) = S3i = 1 fi(x)wi of u(x) where the parameters wi do not coincide in general with the values of u(x) in the centre of the element. The (17) are obtained by forcing the inter-element continuity of u(x) and its derivatives (see Fig. 4). The extension of this to the two-dimensional case via (16) retains the C1 continuity in the whole domain using a minimal number of parameters. It is worth mentioning that such an interpolation scheme uses very simple algebra, which involves only small (3*3) matrices. In fact, the basic integrals required for the elemental matrices is expressed in the following simple form: Z 2 3 X @ u…x1 ; x2 † @2 u…x1 ; x2 † sq dO ˆ Arp p q ih Ajk wij whk @xr1 @xs2 @x1 @x2 Oi ijhk r; s; p; q ˆ 0; 1; 2; where

4. Numerical tests This section reports some numerical tests referring to plate vibration and buckling problems. In both cases a discretization based on the high continuity (HC) ®nite element interpolation proposed by Aristodemo in [17] (see also [15]) have been used. The element has been already used in [3, 9, 14, 18, 19], within a multigrid context, and in [17, 20], within a bucling analysis context. HC discretization assumes rectangular elements for the domain subdivision and interpolates the displacements ®eld by the following bi-quadratic shape functions:   3 X 3 X 1 1 u…x1 ; x2 † ˆ fi …x1 †fj …x2 †wij x1 ; x2 2 ÿ ; ; 2 2 iˆ1 jˆ1 …16† where interpolation parameters wij are located at the centre of the elements (see Fig. 4) and the basic shape functions fi(xk), are de®ned as follows:

…17†

Fig. 4. HC ®nite elements.

S. Lopez et al. / Computers and Structures 69 (1998) 625±637

631

Fig. 5. HC ®nite elements.

Ars ij :ˆ

Z

s dr fi …x† d fj …x† dx r dxs ÿ1=2 dx 1=2

r; s ˆ 0; 1; 2:

The interested reader can refer to [20] for further details. Using HC elements, the re®nement is performed by element bisection (a coarse element on Gh produces four ®ne elements on Gh + 1). Canonical transfer operators are directly de®ned by the interpolation law. By reference to Fig. 5, the interpolation (coarse-to-®ne) operator is expressed as: 2 3 F1 6 F2 7 6 7 2 2 3 36 F3 7 7 fa 1 3 0 3 9 0 0 0 0 6 6 7 6 0 3 1 0 9 3 0 0 0 76 F4 7 6 fb 7 1 6 7ˆ 6 76 7 4 fc 5 16 4 0 0 0 3 9 0 1 3 0 56 F5 7 6 F6 7 7 0 0 0 0 9 3 0 3 1 6 fd 6 F7 7 6 7 4 F8 5 F9

Two analyses have been performed using two di€erent tolerances (E = 10ÿ3 and E = 2.0  10ÿ4) for both the discretization and the solution error. The results are reported in Fig. 7 and in Tables 1 and 2. Fig. 7 shows the ®nal 4th level grid used in the ®rst analysis and the eigenmodes obtained. For an easier representation, the latter are shown using a quite coarse uniform visualization grid unrelated to that used in the analysis. Tables 1 and 2 report, for each level of discretization, the total number of nodes, the number of nodes actually considered in the local re®nement process, the CPU time (in seconds) required and the normalized values of the eigenfrequencies oi. The latter are related to the eigenfrequencies oi by the

the restriction (®ne-to-coarse) operator being its transposition. Numerical computations have been performed using a SUN (model Sparc 1+) work-station. For all tests, `reference' results have been obtained using extremely ®ne uniform grids. 4.1. Test 1 The ®rst test (see Fig. 6) refers to a vibration problem. A square Kirchho€ plate of uniform thickness, simply supported on two opposite edges and clamped on the other edges is considered and its ®rst three eigensolutions are computed.

Fig. 6. Test 1.

632

S. Lopez et al. / Computers and Structures 69 (1998) 625±637

formula

p D=m i oi ˆ o a2 a, m and D being the side length, the unit mass and the ¯exural rigidity of the plate, respectively. For the same problem, the reference values, corresponding to the 6th level of uniform re®nement, are  1 ˆ 28:951283; o

 2 ˆ 54:743908; o

 3 ˆ 69:329619: o

Note that while a reduction to less than 45% in the number of unknowns is provided by the local re®nement, the accuracy of oi of about 2  10ÿ3 for the ®rst analysis and 3  10ÿ4 for the second analysis agrees with the required tolerances (we expect a total relative error of the order of 2E). 4.2. Test 2 The second test refers to the free vibrations of a rectangular Kirchho€ plate simply-supported on the boundary ABDE, and free on the other boundary as shown in Fig. 8. Due to the presence of signi®cant discontinuous boundary conditions on sides ABC and

DEF, the problem presents singularities in the solution. The analysis has been performed with E = 5.5  10ÿ4 requiring the ®rst 4 eigensolutions. The results, shown in Figs. 9 and 10, report the ®nal 4th level grid and the four eigenmodes (as in the previous case, a quite coarse uniform grid has been used for the visualization). Table 3 reports the numerical results of the analysis in the same order as in the previous test. Reference values for the normalized values of the eigenfrequencies oi have been obtained using a uniform grid corresponding to the 5th level of re®nement:  2 ÿ1 p Eh3 lcr  2 b 12…1 ÿ n2 † The relative error of about 1  10ÿ3 for oi agrees with the assumed tolerance E. The highly appreciable reduction to 4% provided in this case by the local re®nement relates to the presence of singularities. 4.3. Test 3 The third test, shown in Fig. 11, refers to a buckling problem. A rectangular von Karman plate simply supported for out-of-plane displacements and with free

S. Lopez et al. / Computers and Structures 69 (1998) 625±637

633

Table 1 Test 1, E = 1.0  10ÿ3 Lev. 0 1 2 3 4

unif. 49 144 484 1764 6724

loc. ref. 49 144 484 1732 6644

CPU

o1

o2

o3

0.26 1.03 2.81 8.92 19.94

30.952815 29.407317 29.062962 29.984148 28.979659

58.451381 55.617783 54.960572 54.834859 54.832673

83.264869 72.164640 70.010862 69.575215 69.503072

CPU

o1

o2

o3

30.952815 29.407318 29.062554 28.978710 28.958038 28.954842

58.451381 55.617511 54.958869 54.797182 54.762361 54.758575

83.264869 72.163499 70.004378 69.495058 69.381476 69.355199

Table 2 Test 1, E = 2.0  10ÿ4 Lev. 0 1 2 3 4 5

unif. 49 144 484 1764 6724 26244

loc. ref. 49 144 484 1764 6644 11748

0.38 1.48 4.43 11.28 38.18 102.65

edge conditions for in-plane displacements is considered. The plate is compressed in the longitudinal direction by a normal load distribution of value (1  l) along the transversal edges. The Young modulus E = 2.1  106, Poisson coecient n = 0.25 and the thickness of the plate h = 0.1 are considered. The exact analytical solution of the critical load lcr=4  p2/b2 Eh3/12 (1 ÿ n2) is known from literature. The analysis has been performed using two di€erent tolerances E = 10ÿ3 and E = 2.5  10ÿ4.

In Fig. 12 the ®nal 4th level grid of the ®rst analysis is shown while the buckling mode is visualized using a coarse uniform grid. Tables 4 and 5 report the numerical results for both analyses. A reduction to less than 40% is provided by the local re®nement while the accuracy of the critical load of about 5  10ÿ4 for the ®rst analysis and 2  10ÿ4 for the second analysis, agrees with the assumed tolerances.

Fig. 8. Test 2.

Fig. 9. Test 2, E = 5.5  10ÿ4, Mm grid.

634

S. Lopez et al. / Computers and Structures 69 (1998) 625±637

Fig. 10. Test 2, E = 5.5  10ÿ4, eigensolutions.

4.4. Test 4

level of uniform re®nement. The relative error for the critical load of about 2  10ÿ3 agrees with the assumed tolerance E. Due to the presence of singularities and to the monomodal analysis, the reduction to 2% provided by the local re®nement is highly appreciable see Table 6.

The fourth test (see Fig. 13) refers to the buckling of a rectangular von Karman plate with the same load, the same mechanical features and the same conditions for in-plane displacements as the previous example. For out-of-plane displacement the discontinuous boundary conditions of Test 2 are considered. The results of the analysis, performed with E = 2.5  10ÿ4, are shown in Figs. 14 and 6. Fig. 14 shows the ®nal 5th level grid of the analysis and the buckling mode, visualized by means of a quite coarse uniform grid. The reference value lcr=0.177041 for critical load corresponds to the 6th

5. Conclusions A basic multigrid algorithm has been used for the ®nite element analysis of eigenvalue problems and sample applications have been made for plate vibration

Table 3 Test 2, E = 5.5  10ÿ4 Lev. 0 1 2 3 4

unif. 100 324 1156 4356 16 900

loc. ref. 100 324 1156 1460 710

CPU 12.21 18.02 29.07 50.76 87.39

o1 6.957895 6.780734 6.699536 6.662492 6.639291

o2 18.490673 17.957706 17.729533 17.630643 17.598406

o3 21.756361 21.638203 21.604183 21.599775 21.600098

o4 40.100699 39.309465 39.002852 38.890625 38.858792

S. Lopez et al. / Computers and Structures 69 (1998) 625±637

635

Fig. 11. Test 3.

and buckling problems. The algorithm has the following features: . The combination between the residual eigenvalue search and the multigrid solution provides a powerful tool for eigenvalue analysis. . The analysis is fast, in comparison with standard matrix methods, and allows very ®ne discretizations with accurate results for the eigensolutions. . The described data structure is simple, ¯exible and suitable for adaptive local re®nement.

Some further comments can be made: 1. While the algorithm has been described in relation to a linear eigenvalue problem, it can be easily extended to nonlinear eigenvalue searches such as that required by asymptotic post-buckling analysis [17]. Its ability to provide very ®ne solutions makes the multigrid approach highly suitable for this type of analysis [3]. The nonlinear extension only in¯uences the computation of the primary residual [14].

Table 4 Test 3, E = 1.0  10ÿ3

lcr  Lev. 0 1 2 3 4

un. 55 160 532 1924 7300

loc. ref. 55 160 532 1924 3012

 2 ÿ1 p Eh3 2 2 b 12…1 ÿ n †

CPU 0.73 1.41 2.15 5.64 8.55

4.192693 4.046320 4.011469 4.002859 4.002059

636

S. Lopez et al. / Computers and Structures 69 (1998) 625±637 Table 5 Test 3, E = 2.0  10ÿ4

lcr  Lev. 0 1 2 3 4 5

un. 55 160 532 1924 7300 28420

loc. ref. 55 160 532 1924 7252 9828

 2 ÿ1 p Eh3 2 2 b 12…1 ÿ n †

CPU 0.71 1.25 2.55 6.06 25.51 44.79

4.192693 4.046319 4.011463 4.002859 4.000796 4.000701

re®nement is transparent and does not in¯uence the complexity of the solution process. 3. Multigrid solution is based on repeated local computations performed at the element level. The simplicity of these computations, that is, a simple and fast algebra in the element description, plays an important role in the overall eciency of the process. Rough, small and simple elements tend to be more ecient than large, complex and accurate ones. This is just the opposite of what happens in standard matrix-based solution procedures where the computational burden is strictly related to the total number of unknowns. Fig. 13. Test 4.

2. The local re®nement option provides an appreciable reduction in both the computation time and storage demand. Obviously, the e€ect of the re®nement test (15), when extended to all the modes under consideration, tends to decrease with the number of modes to be analyzed. However, due to the unstructured, pointer-based organization of data, a partial

HC ®nite elements may be a good compromise for multigrid applications due to both the element simplicity and the relatively small number of variables required for accurate discretizations (at least when the domain to be discretized is regular enough). 4. Only regular domains are treated in the applications, in order to test the essential features of the algorithm in a simple context. The extension to general irregular domains would require some sort of

Fig. 14. Test 4, E = 2.5  10ÿ4, Mm grid and critical mode.

S. Lopez et al. / Computers and Structures 69 (1998) 625±637

637

Table 6 Test 4, E = 2.5  10ÿ4

lcr  Lev. 0 1 2 3 4 5

un. 60 180 612 2244 8580 33 540

loc. ref.

CPU

60 180 612 1852 1684 690

companion mesh-generation algorithm. Multigrid versions of generators based on coordinate transformations such as the Thompson±Thames algorithm [21] are possible candidates (see also [22]).

0.66 1.69 2.92 6.13 20.13 27.56

[12]

[13]

References [1] Lanzo A, Garcea G, Casciaro R. Asymptotic post-buckling analysis of rectangular plates by HC ®nite elements. Int J Numer Meth Engng 1995;38:2325±45. [2] Lanzo A, Garcea G. Koiter's post-buckling analysis of thin-walled structures by a ®nite element approach. Int J Numer Meth Engng 1996;39:3007±31. [3] Lopez S, Fortino S, Lanzo AD. Post-buckling analysis of plates by multigrid FEM approach. In: ECCOMAS `96 Conference. Wiley, Parigi, 1996. [4] Brandt A. Multi-level adaptive solutions to boundary value problems. Math Comp 1977;31:333±90. [5] StuÈben K, Trottenberg U. Multigrid methods: fundamental algorithms, model problem analysis and applications. In: Hackbusch W, Trottenberg U, editors. Multigrid Methods. Springer, Berlin, 1982:1±176. [6] Hackbusch W. Multi-grid eigenvalue problem. In: Advances in multi-grid methods. Proceedings Conference. Oberwolfach, Friedr. Vieweg-Sohn, 1984. [7] Bank R. Analysis of a multilevel inverse iteration procedure for eigenvalue problems. SIAM J Numer Anal 1982;19:886±93. [8] Hwang T, Parsons ID. Multigrid solution procedures for structural dynamics eigenvalue problems. Comput Mech 1992;10:247±62. [9] Lopez S, Casciaro R. Algorithmic aspects of adaptive multigrid ®nite element analysis. Int J Num Meth Engng 1997;40:919±36. [10] Lopez S. Aspetti algoritmici dei metodi multigrid in problemi di meccanica computazionale' (in italian). Ph.D. Thesis. Unical, Cosenza, Italy, 1994. [11] Casciaro R, Aristodemo M. Perturbation analysis of geometrically nonlinear structures. International Conference

 2 ÿ1 p Eh3 2 2 a 12…1 ÿ n †

[14]

[15]

[16]

[17]

[18]

[19]

[20]

[21]

[22]

0.190381 0.183172 0.179911 0.178388 0.177569 0.177510

on Finite elements nonlinear solid and structural mechanics, Gelio, Norway, 1977. Zienkiewicz OC, Zhu JZ. A simple error estimator and adaptive procedure for pratical engineering analysis. Int J Numer Meth Engng 1987;24:337±57. Zienkiewicz OC, Zhu JZ. Error estimates and adaptive re®nement for plate bending problems. Int J Numer Meth Engng 1989;28:2839±53. Lopez S, Fortino S. Soluzione adattativa multigrid di problemi bidimensionali elastici (in italian). Report No. 164. Dip. Strutture, Unical, Cosenza, Italy, 1995. Rivara MC. Design and data structure of fully adaptive, multigrid, ®nite element software. ACM Trans Math Soft 1984;10:242±64. Rivara MC. A grid generator based on 4-triangles conforming mesh-re®nement algorithms. Int J Numer Meth Engng 1987;28:1343±54. Storti M, Nigro N, Idelsohn S. Multigrid methods and adaptive re®nement techniques in elliptic problems by ®nite element methods. Comp Meth Appl Mech Engng 1991;93:13±30. Lopez S, Fortino S, Casciaro R. Multigrid analysis of two-dimensional eigenvalue problems. Report no. 169. Dip. Strutture, Unical, Cosenza, Italy, 1995. Lopez S, Fortino S, Casciaro S. Analisi multigrid di piastre alla Kirchho€ (in italian). Proceedings of XII Congresso AIMETA, Napoli, Italy, 1995. Aristodemo M. A high-continuity ®nite element model for two-dimensional elastic problems. Comput Struct 1985;21(5):987±93. Thompson JF, Thames FC, Waine Mastin C. Automatic numerical generation of body-®tted curvilinear coordinate system for ®els containing any number of arbitrary two-dimensional bodies. J Comp Phys 1974;15:299. Aristodemo M, Casciaro R, Vencia G. Analysis of twodimensional elastoplastic continua by high-continuity ®nite elements. In: Taylor et al., eds. Numerical methods for nonlinear problems. Pineridge Press, Swansea, 1982.