High precision finite elements

FINITE ELEMENTS IN ANALYSIS A N D DESIGN

ELSEVIER

Finite Elements in Analysis and Design 26 (1997) 97-114

High precision finite elements P.H. Kulla RESULT, Hochholz 2, D-83556 Griestaett, Germany

Abstract

The finite elements method has established itself as a synonym for Structural Analysis over the years. Nevertheless, the accuracy of FE-results continues to be unsatisfactory in many cases. The error is related to the fact that finite elements only approximate the partial differential equations describing the structure, where the quality of the approximation depends on the elementation. In this paper a new type of finite element is presented, which is based on exact solutions of the PDEs. This element is in fact the dynamic stiffness matrix of a continuous substructure that has not been discretized. Its elements are transcendental functions of the frequency. Special tools have to be used to find the eigensolutions. On the other hand, the matrix accurately describes the behaviour of a possibly very large structural element up to high frequencies. The elementation of the structure is solely determined by its geometry and not by precision constraints. In the past, only one-dimensional continuous FEs have been available under different names including "continuous elements" and "Exact Finite Elements". Recently, the method has been modified to include two-dimensional elements. The modification does introduce an error, but still the performance compares favourably with the traditional FEM. Accordingly, the technique is termed high precision finite elements method.

Keywords." Continuum methods; Mechanical impedance; Dynamic stiffness matrices; HPFEM

1. Introduction Mechanical systems are continuous in three spatial coordinates and in time. For the purpose of analysis the actual system is replaced by a theoretical model. The complexity chosen for this model depends on the accuracy desired for the results and on the resources available for the investigation. The spectrum reaches from three-dimensional continuous models to simple lumped-parameter approximations. After transformation to the frequency domain, continuous models are described by partial differential equations (PDEs) plus boundary conditions, whereas discrete models are described by sets of algebraic equations. To avoid the solution of PDEs, discrete models (e.g. traditional finite elements) are preferred in general. The direct treatment of continuous systems, though historically older, has E-mail: [email protected]. S0168-874X/97/$17.00 @ 1997 Elsevier Science B.V. All rights reserved

PII S 0 1 6 8 - 8 7 4 X ( 9 6 ) 0 0 0 7 3 - X

98

P.H. Kulla / Finite Elements in Analysis and Design 26 (1997) 97-114

never reached the degree of popularity of the FE method. The FE method is often considered to be more flexible than methods based on the direct solution of continuous systems. The key idea behind the method to be presented is to combine the flexibility of the FE approach with the increased accuracy of continuous methods. Continuous finite elements are not new but have been in use under different names for quite some time [1, 2]. However, the scope of the method has been restricted to one-dimensional continua, e.g. beams. With a few exceptions, the method has never been automated and cast into program modules, that would allow an average user to profit from it without going into the details of the approach. Rather, it was generally used for the individual solution of special problems by specialized people. For some classes of problems, however, the method has been furnished with more or less user-friendly surfaces and is in practical use with a few institutions. Among these are Chalmers University, Sweden [3], University of Wales [4], NASA/LRC [2], ESA/ESTEC [5, 6], STCAN, France [7] and RESULT, Germany. All these program packages give "exact" results, the last one also for non-conservative systems. Their main drawback is their lack of two-dimensional elements. This shortcoming has been a major obstacle to a broader acceptance of the otherwise powerful method. Two-dimensional continuous elements have been introduced to overcome this shortcoming [8]. This generalisation, however, also includes one marked drawback: For one-dimensional continua the boundary conditions can always be satisfied in an exact way following standard procedures. For multidimensional continua the boundary conditions have to be approximated in general, which introduces an error and abandons the label "exact". In our case, the continuous functions representing displacements and forces along the edges are decomposed into Fourier-type series. Truncating these series introduces an error, quite in contrast to the one-dimensional case. On the other hand, this drawback is not too serious from a practical point of view, since the truncation point can be easily controlled by the user and the effect of the truncation on the result is readily evaluated without re-elementation of the structure. The method may be classified as an approximation method with high (and easily tunable) precision. These features are reflected by the name high precision finite element method (HPFEM). The precision of the method is one of its strongest points. As a consequence, arbitrarily and even infinitely large elements may be used without loss of accuracy. The elementation of the structure is determined by its geometric properties rather than by a certain a priori knowledge of the result to be obtained. Geometrically, simple structures need not be cut down to pieces in order to get reasonable precision. There are no hidden internal contradictions that may lead to sudden blocking eventually. Precision is a topic that has traditionally been treated in a somewhat casual way. It is sometimes claimed that the errors of the analysis are completely masked by uncertainties inherent in the model or in the assumed loads. This may be true, but still it is mandatory to scrutinize the accuracy in order to eliminate cases, where many small errors combine in a way which gives completely useless results. Typically, these error-sensitive situations are encountered when evaluating stress concentrations, higher eigenforms or stability limits. The direct solution (without discretization) of PDEs has been critisised as being too complicated, academic and of little use for the solution of practical problems. In addition, it was claimed that the method is based on a separation-of-variables technique and hence ignores all solutions which do not fit this specific product form. Without going into details here, it is stated that all these arguments miss the point: The general solution of (linear) PDEs itself does not pose any serious problem, as

P.H. Kulla/Finite Elements in Analysis and Design 26 (1997) 97-114

99

will be shown. The coding of the subroutines delivering the element matrices is straightforward. Travelling wave effects are covered by HPFEM just as well as by traditional FEM, but with much higher precision. The main problem to be solved is the satisfaction of the boundary conditions. On the element level, the boundary conditions are not yet decided on. Consequently, the solutions of the PDEs used to compose continuous finite elements must be unspecific w.r.t, boundary conditions. This excludes all specific choices such as "admissible", "natural", "free" or "clamped" functions. In the literature [9-11], in contrast to the above application, exact solutions of PDEs are generally used to generate solutions for systems with well-defined boundary conditions. In this case the solutions of the PDE should a priori incorporate as many of the system boundary conditions as possible. Our objectiw~ is different, however, since we want to use the elements for arbitrary boundary conditions. Nevertheless, the generation of high precision finite elements presented here is based on these earlier works and may be considered a generalization of the results given there. Using continuou,; elements, the behaviour of a system is described by a dynamic stiffness matrix, the elements of which are transcendental functions of the frequency. Special tools are necessary to evaluate the related eigensolutions [4, 12, 13]. On the other hand, the lack of discretization errors will not only give more accurate numerical results, but may also lead to a broader understanding of the theory of (structural) dynamic systems. Continuous elements are close to boundary elements since both methods fulfill the PDEs accurately. The methods use different types of base solutions and different procedures of discretizing the boundary conditions. Both methods claim high precision.

2. One-dimensional continuous elements

A typical one-&~mensional continuous substructure is a beam; see Fig. 1. Its deflection w is described as a function of locus x and time t by (EIw")" + # ~ = p ( x , t )

with - a<~x<<.a

(1)

where E1 is the bending stiffness, # is the mass per unit length, p(x, t) is the external load and 2a is the length of the beam. Restricting ourselves to the homogeneous case (p(x, t) = 0) and assuming that the coefficients. E1 and /t are constant over the beam length, the solutions of (1) can be given (without loss of generality) under the form w(x, t) = ~epx e st

=

ff;epx+st.

(2)

Inserting (2) into (1) and splitting off the terms contained in (2) we obtain the characteristic equation or dispersion equation of the beam: EIp 4 +

fls 2 =

0.

(3)

The ansatz (2) has reduced the PDE (1) to the algebraic equation (3). Having in mind the harmonic motion with circular frequency ~o (which is standard in the frequency domain) we write co2 = - s 2.

(4)

P.H. Kulla / Finite Elements in Analysis and Design 26 '1997) 9 ~ 1 1 4

100

Q

_k_ I

Q ---~

~w 0

I~

x

I d1 do

d2

ol 3

Fig. 1. Beam element.

Solving (3) for p we get four solutions P0-3: P0,1 = :~0~, P2,3 =-+-i~

with i = x/L-i-

(5)

where 4/~

(6/

Inserting the four solutions (5) into (2) gives four base solutions of (1). Omitting the (harmonic) dependency on time and grouping conjugate complex terms in x to real expressions, the four base solutions are

Wo(X) = ~0 cosh ~x, wj(x) = #l sinh ~x,

w2(x) = ~2 cos ~x, w3(x) = w3 sin ~x.

(7)

If the parameters in (1) do depend on x, four (linearly independent) base solutions equivalent to (7) can always be generated by numeric integration of (1) with four linearly independent starting values at x = 0. The only difference to the analytical solution (7) is the runtime. The deformations and forces at the two boundaries of the beam are now expressed in terms of the base solutions:

Db = d,

(8)

Fb = f ,

(9)

where

D=

cosh - ~ a ~sinh - ~ a cosh ~a ~sinh ~a

sinh-~a cosh - ~ a sinh ~a cosh ~a

cos-c~a -c~ sin - ~ a cos ~a - ~ sin ~a

sin-~a | ~ cos -c~a sin ~a ' ~ cos ~a

(10)

P.H. Kulla/Finite Elements in Analysis and Design 26 (1997) 97-114

sinh -c~a - cosh -eta - ~ sinh cca cosh eta

F = EIe 2

c~cosh - c a - sinh -eta -c~ cosh ~a sinh ~a

~ sin -c~a - cos -eta -c~ sin c~a - cos ~a

-ct cos -cta ] / sin - ~ a / ~ cos ~ta | - sin c~a t.l

101

(11)

1,~0 b =

w2

'

(12)

~3

d =

f =El

w(-a) w'(-a) w(a) w'(a)

w"'(-a) -w"(--a) -w'"(a) w"(a)

(13)

(14)

The signs in (13) have been chosen in accordance with the DOFs of the end-points introduced in Fig. 1. Note that the same sense of rotation has arbitrarily been chosen for rotations at the fight and the left ends of the beam. The signs for the forces in (14) follow from the rule that the forces are adjoint to the displacements or, in other words, that d T f is the "dynamic potential" or the energy pumped into the system at the boundaries. Premultiplying (81) with D -1 and inserting b into (9) gives F D - t d =- Z d =:f ,

(15)

where Z is the dynamic stiffness matrix of the beam element. It may be used as an exact FE. For conservative systems it is symmetric. The symmetry of Z is due to the self-adjointness of the underlying differential equation plus boundary conditions. The derivative of Z w.r.t, s2(= -092) is positive definite. Due to the inversion performed in (15) it has poles. These poles are the natural frequencies of the element clamped at the boundaries. The eigensolutions of (1 5) (zero right-hand side) represent the natural frequencies and eigenforms of the free element. Since no discretization has been used, (15) is exact and contains an infinity of eigensolutions. It can therefore not be split up into frequency-independent mass and stiffness matrices. Evaluating (15) for one frequency implies a call to a subroutine which composes F and D and processes them according to (15). For every individual frequency, the subroutine has to be called again. For more elaborate presentations of dynamic stiffness matrices of beams see, e.g., [1]. The assembly of elements to a more complex structure is essentially identical for all methods of structural analysis based on stiffness matrices, i.e. the assembly follows the lines known from traditional FEM. The resulting matrix represents the system as assembled of elements. This system matrix has the specific properties mentioned above with the element matrices. The system matrix may depend on certain parameters, e.g. frequency or prestress. The specific values for these parameters,

102

P.H. Kulla/Finite Elements in Analysis and Design 26 (1997) 97-114

which turn the system matrix singular, are eigenvalues. Prominent types of eigenvalues are natural frequencies or bucklin9 loads. Severe numerical problems will occur if the eigenvalues are located by simply scanning the determinant, which possibly has further hindered a broad acceptance of the method. Special tools are available for the solution of this transcendental eigenvalue problem [13, 12, 4]. After solution of the system equations, the complete state of the system can be determined for every point of every element by back-substitution of the solution of the system equations into the DEs of the beam elements.

3. Two-dimensional high precision finite elements Continuous elements as described above are accurate. This means that only numerical (and no methodical) errors are present in the result. In other words, the error is only due to the finite word length of the compiler in use. (Note, however, that poor conditioning still may turn the result completely useless.) Once the merits of one-dimensional continuous elements have been realized, it is natural to ask whether the method can be generalized to include two-dimensional elements. Obviously, the availability of such elements would enormously broaden the scope of applicability and would make the method a true general-purpose tool. The concept of composing the system stiffness matrix from element submatrices implies that the interface and boundary conditions are imposed only at the system level, i.e. the solutions for the elements have to cope with all potentially arising boundary conditions. To make the point quite clear: Our base solutions must not fulfill any well-defined set of boundary conditions. In contrast to one-dimensional substructures (having a finite number of boundary conditions), two-dimensional substructures have continuous boundary conditions. Fulfilling continuous boundary conditions is comparable to fulfil discrete boundary conditions on an infinity of points. This is a major obstacle to an exact solution of the problem. Only for very specific boundary conditions can exact solutions be given. In general, the continuous boundary conditions have to be approximated by a series of discrete conditions. The discretization, performed by whatever method, invariably introduces an error. In contrast to the one-dimensional case, the method can no more be called "exact", and the question arises, whether the increased effort at the element level (as compared to traditional finite elements) is still worthwhile. Ultimately, the answer can only be given by comparing the results in terms of precision and runtime. The error incurred by discretization may slightly destroy the symmetry of the dynamic stiffness matrix and the positive definiteness of its derivative w.r.t, s 2. This is serious, since it alters basic physical properties of the system, e.g. conservativeness. As a consequence, "phantom modes" as reported in the literature may be detected. Here we cut the discussion of this point by stating that both typical properties (symmetry and positive definiteness of the derivative w.r.t, s 2) can be conserved by using special techniques. This will safely avoid erraneous modes. The resuiting modes will be approximations, though, but the precision of the modes (and inspecific the related fields for bending moments and transversal force) will compare favourably with competing methods.

P.H. Kulla/Finite Elements in Analysis and Design 26 (1997) 97-114

103

4. Differential equations Here we restrict ourselves to Kirchhoff plates. The principles presented here are directly applicable to other types, e.g. Mindlin plates. The PDE describing transversial displacement w is D A A w + It# = p(x, y, t)

(16)

with -a<~x<.a

and

-b<.y<<.b,

(17)

where D is the plate's bending stiffness D =

Eh 3

12(1 - v2)

(18)

with Young's modulus E, plate thickness h, Poisson's ratio v, mass distribution # and surface load p. Dots stand for the derivation w.r.t, time and the Laplace operator is defined in Cartesian coordinates as A( ) = ( ),xx+( ),yy.

(19)

Eq. (16) is a PDE with three independent variables x, y and t. With respect to x and y it describes a boundary value problem and w.r.t, t it describes an initial value problem. In the following, only the homogeneous pilate (with p(x, y , t ) = 0) is considered. The generalization to include the loaded plate is straightforward [14]. The equations reilating the deflection w and its partial derivatives to the forces and moments in the plate domain (and on its boundaries) are given in standard text books, e.g. [15]. The inclination of the surface is defined by ~x = -w~,

(20)

fly = -w,y.

(21 )

The moments for bending and twist (per unit length) are related to w by m~ = - D ( w ~ x + YW,yy),

(22)

mxy = - D ( 1 - v)Way,

(23)

my = -D(w,yy + VWax).

(24)

Due to Kirchhoff's approximation (no shear deformation), the shear force is combined with the derivative of the twist moment to form the transversal force v: + (2 - V)W yy),~,

(25)

t)y = - D ( w , y y -.[- (2 - v)W~x),y.

(26)

vx = - D ( w ~

The combination of the derivative of the twist moment with the shear force also gives rise to singular point forces in the four comer points. The point force in the comer (x = a, y = b) is Fc = - 2 m x y .

(27)

104

P.H. KullalFinite Elements in Analysis and Design 26 (1997) 97-114

The approximative character of Kirhoff's plate theory is also reflected by the fact that the transversal force as given by (25) and (26) is not a vector field. This means that the field is not independent of the coordinate system chosen to resolve the vectors into components.

5. Generation of base solutions

Since PDE (16) is linear with coefficients constant w.r.t, time and locus, the general solution to the homogeneous case (p(x, y, t) = 0) is of the form w(x, y, t) = ~epXeqYe "t = l~e px+qy+st,

(28)

where s is the frequency parameter and, analogously, p and q are wave parameters. Inserting (28) into (16), prescribing zero for the right-hand side and splitting off the terms contained in (28) we get the characteristic equation of the plate. D(p2 + q2)2

+ ]2s 2 :

0.

(29)

Taking the square root and replacing s 2 by -092 we get the characteristic equation of the form p2 d- q2 = ±CO,/ # .

(30)

VEI

The two parallel inclined lines in Fig. 2 represent the solutions of (30) in the pZ_q2 plane for plate parameters equalling unity. Any point on either line is also a solution of (29) and hence represents a solution of the PDE (16). An infinite series of base solutions is easily generated by introducing for p2 the series p~ = - (/a~)2 = -c~,

i = 0,1,2...

(31)

Series (31) is indicated in Fig. 2 by the vertical dashed lines with index i. Taking the square root of (31 ) we get Pi = +i~i

with

i = x/~-l.

(32)

Inserting (32) into the factor of (28) containing x and recombining the resulting two exponential terms gives two series of harmonic functions in x: 1

~(exp ieix + exp i7ix) = cos ctix, 1

~ ( e x p ic~ix - exp i~ix) = sin ctix.

(33) (34)

Thus, every term p~ of (31) represents a pair of functions in x, one symmetric and one antisymmetric. Further, for every term p~ of (31) we get two terms from (30), q~0 and q~l; see also Fig. 2. The square terms q2 may be negative or positive, having conjugate complex or real roots. The related

P.H. Kulla/Finite Elements in Analysis and Desiyn 26 (1997) 97-114

105

Cl ~ I2

I I

1 I i

=

0

p~

2

Fig. 2. Solutions of characteristic equation.

functions in y may be brought to real forms as shown above. Depending on the sign of q~, the resulting functions will be cyclic or hyperbolic cosine and sine functions in y. A second series of base solutions is generated by introducing for q2 in (30) the series q} = - ( b ) 2 = -fir,

j = 0,1,2...

(35)

Relenting for a moment, we state that every term p2 in (31 ) corresponds to a vertical line in Fig. 2 containing two pairs of solutions p~, q~ of the characteristic equation (29). Every pair represents four solutions of the PDE (16), where every solution represents one out of the four possible symmetry cases w.r.t, x and y: sym-sym, sym-ant, ant-sym and ant-ant. The same is true for the terms q~ of series (35), where each term corresponds to a horizontal line in Fig. 2. A minimum practical set of base solutions is obtained if i and j are restricted to the first terms of the Fourier series (31) and (35), i.e. i = 0 and j = 0. This results in a minimum set of 16 base solutions. The deflection of the plate can now be written in terms of the base solutions. We restrict the presentation to the double symmetric case; the other symmetry cases follow from analogy. Omitting the (harmonic) dependence on time we get

m--1 Wss(X' Y) = E l (l~iO c o s h ~ioY ~- l~il i=0

ircx cos/~il

Y) cos

-a

+ IY'~(~.o cosh ~,oX + ~), cos ~,ly) cos J;Y. .i=0

(36)

P.H. Kulla/Finite Elements in Analysis and Design 26 (1997) 97-114

106

In (36) it was assumed that one of the functions contained in parantheses was hyperbolic while the other was cyclic. Inspecting Fig. 2 it may be inferred that both functions may become hyperbolic, while both never become cyclic simultanously. The inclinations, moments and forces related to every base solution of series (36) are derived according to (20)-(26). Restricting the presentation to the first term of (36) and omitting the indices "ss" and "07, we have iT~x

wi(x, y) = wi cosh [3iy cos - -a ,

(37)

ire iTzx Gi(x, y) = wi cosh fiiY--a sin a ,

(38)

~ i ( x, Y ) = --Wi~i sinh fliY COS

mxi(X,y):-]~2i(-(i~) 2 mxyi(X, y) = wi(1

-

-

(39)

,

a

+ vfl~)cosh fliY cos --7 ircx a ire

izcx

a

a

V)fl i sinh f l i Y - - sin

myi(X'Y)=ffq(fl~--v(i-~)2) ((i~) 2 Vxi(X,y) = --Wi

iT~x

-- (2 -

1)yi(X,Y) = -wi (fl2i - (2

(40)

,

(41)

c°shfliyc°sizc~xa '

v)fl

~)

iTz coshfliy--

a

(42)

iTzx

sin -a ,

(43)

(44)

-

The technique of generating base solutions shown above was also used in [10, 11]. In these works the solutions were used to investigate the dynamics of plates with specific boundary conditions. Here, in contrast, the solutions are used to generate a high precision finite element (which is capable of coping with arbitrary boundary conditions).

6. Series development of deflections and forces on boundaries The deflections and forces on the boundaries are discretized by development into Fourier-type series. For example, the deflection on the boundary (x,y = b) follows by substituting b for y in (37): iTzx

wi(x, y = b) = wi cosh flib cos - -a According to series (31), this

is already a Fourier series in x.

(45)

P.H. Kulla/Finite Elements in Analysis and Design 26 (1997) 9~114

107

The deflection along the boundary (x = a, y) follows from (37) if x is replaced by a.

wi(x -----a, y) = 1~ i cosh fliY COS ire =

wi

cosh fliY

(-- 1 )i.

(46)

Developing (46) into a Fourier series in y gives Wi(X ~-

a,y)

= w i ( - 1 ) i ~ 1,~ik COS k=0

kzcy b

(47) '

where the Fourier coefficients ffik are determined by =

1 J~_

cosh fliY

krcy COS

dy

(48)

Carrying out the intergration one gets

ffik = ( - 1 )~

2 flib

(/ jb)2 + (k )2

sinh flib.

(49)

The development of rotations (inclinations), forces and moments along the boundaries is performed in an analogous way; see [14]. Since only trigonometric functions are contained in (37)-(44), the formula work is modest. The infinite series (47) is truncated, which introduces an error. The user can determine the truncation point and hence the truncation error. Truncating series (47) after the first term will result in a minimum displacement vector of dimension 4 per side: symmetric and antimetric deflections and rotations. This results in a displacement vector of 16 for the plate, the same is true for the (adjoint) force vector. The minimum dimension 16 for displacement and force vectors corresponds to the minimum number of base solutions introduced earlier. Having developed the deformations and forces of the base solutions into series along the boundaries, the problem i,; reduced to the form of the one-dimensional case given by (8)-(15). In contrast to the one-dimensional case, however, the length of the vectors d, f and b is not only determined by the order of the differential equation, but also by the point of truncation of the series development, i.e. by the desired accuracy. The number of terms in the series development of the displacements and forces according to (47) is chosen to correspond to the truncation of the base solutions (36) in order to make the matrices D and F square. The sequence of lhe elements in the displacement (column) matrix d is arbitrary, as is the sequence of the elements of ~:he matrix b. For the ease of assembly it is advantageous to group the elements of d (and hence f ) according to the four sides of the plate. The rest of the procedure to develop the element matrix Z is identical to the procedure outlined for the one-dimensional case.

7. Example: cantilever plate consisting of two elements With the plate element developed above, arbitrary planar combinations of rectangular plates can be assembled. The assembly technique is exactly as used with traditional finite elements. Instead of using a grid of small elements, HPFEM uses only a few large elements. The HPFEM stiffness matrix relating to a single element is fully populated, as is a traditional FEM element matrix. Since

108

P.H. KullalFinite Elements in Analysis and Design 26 (1997) 97-114

r

22

L)

Fig. 3. Cantilever plate with line load.

typically only a few elements are used with HPFEM, the banded structure of the system matrix known with traditional FEs is less prevailing. We demonstrate the performance of the HPFEM by analyzing a school-book example: A square cantilever plate is loaded by a constant line load along one of its sides; see Fig. 3. The parameters of the system are: 2a = 2b = 1.0; v = 0.30; p = 1.0. The plate is loaded by a line load along one of its edges. The load is P = 1 where P = 2bp with p being the load per unit length. The stiffness of the plate is D = 1 in the right-half of the plate (x > 0) and D = 2 in the left-half. Accordingly, two plate elements were used to model the plate, one for each half-plate. The plate is excited at 3 different frequencies, co = 1, co = 100 and co = 1000. For excitation at the lowest frequency, the program was run with two different levels of accuracy. The results shown in Fig. 4 were obtained using 4 Fourier terms (resulting in a 64 x 64 dynamic stiffness matrix for one element). The results presented in Fig. 5 for the same case were obtained using 6 Fourier terms (96 x 96 matrix). The medium- and high-frequency cases (Figs. 6 and 7) were also run with 6 Fourier terms. The runtime per figure including the generation of all related pictures was of the order of tens of seconds on a desktop.

8. Results and conclusions The figures show from lower left to upper right contour plots for the following scalar fields: displacement w, inclination ~bx, inclination ~by, moment rex, moment mxy, moment my, transversal force vx and transversal force Vy. The contour lines are presented in 10 shades, where the brighter shades correspond to higher function values. For the low-frequency excitation, the scaling for a group of 10 shades is 0.1 for the deflection, 0.2 for the inclinations, 0.5 for the moments and 1.0 for the forces. Different scalings are used for medium- and high-frequency excitation. The results presented here were not subject to any smoothing algorithms. Note the discontinuities along the element interface line in the functions mxy, my and Vy. The solutions fulfill the differential

P.H. Kulla/Finite Elements in Analysis and Design 26 (1997) 97-114

~R

E~

Fig. 4. Cantilever plate, ~o = 1, two 64 × 64 HPFEs.

109

110

P.H. Kulla / Finite Elements in Analysis" and Design 26 (1997) 97-114

E~

Fig. 5. Cantilever plate, co = 1, two 96 × 96 HPFEs.

P.H. Kullal Finite Elements in Analysis and Desiyn 26 (1997) 97-114

._>,

~R

.x

E~

Fig. 6. Cantilever plate, o~ = 100, two 96 × 96 HPFEs.

111

112

P.H. Kulla/Finite Elements in Analysis and Design 26 (1997) 97-114

X'~ >0.

E~

Fig. 7. Cantilever plate, co -- 1000, two 96 x 96 HPFEs.

P.H. Kulla/Finite Elements in Analysis and Design 26 (1997) 97-114

113

equation o f the plate exactly. The boundary and interface conditions, however, are only fulfilled in an approximate way. The deflection, for instance, oscillates between small positive and negative values along the clamped edge. Also, the transversal force Vx oscillates around the value 1.0 along the loaded free edge. The error observed along the boundaries quickly decays towards the inner o f the element domain. We conclude the presentation o f the HPFEM by listing a few practical aspects of the method: • The HPFEM giw~s results that are unrivaled by traditional techniques in terms of precision and runtime. • Potential fields of application are high frequency and/or stability analysis. Also, the method can be used to verify FEM results without using FEM. • The method seems ideal for fast analysis of systems in the design phase. • The method enables interactive comparison of analysis and test (based on full field laser measurement techniques). • Presently, (two-dimensional) HPFEs are available for Kirchhoff- and Mindlin-type plates. Plates undergoing in-plane deformations are also available. • Provided with a user-friendly surface, the method may serve as a supplement to commercially available FEM packages.

Acknowledgements The author developed the HPFEM in the sequel o f ESA/ESTEC Study Contract No. 9562/91/NL/ JG(SC). He gratefu![ly appreciates the cooperation of the ESTEC staff involved in the contract.

References [1] V. Kolou6ek, Dynamics in Engineering Structures, Butterworths, London, 1961. [2] T.H. Broome, "The VIBIC1 computer program", Presented at the 49th Annual Convention of the National Technical Association, Hampton, VA, 1977. [3] B.A. Akesson, "PFVIBAT - a computer program for plane frame vibration analysis by an exact method", Int. J. Numer. Method~ Eng. 10, 1976. [4] F.W. Williams, "Computation of natural frequencies and initial buckling stresses of prismatic plate assemblies", J. Sound Vib. 21, 1972. [5] D. Poelaert, "DIST?EL, a distributed element program for dynamic modelling and response analysis of flexible structures", Proc. 4th VPI & SU/AIAA Symp. on Dynamics and Control of Large Space Structures, Blacksburg, VA. 1983. [6] F. Janssens and E. Crellin, "Analytical modal analysis of flexible spacecraft", Paper Presented at the 8th Seminar on Modal Analysis, KU. Leuven, Belgium, September 1983. [7] C. Duforet, H. Grungier and J.P. Quintin, "Correlation calculus-mesnre de la dynamique d'une structure tubulaire jusqu'/~ des frbquerLces~lev6es", STRUCOME'88. [8] P.H. Kulla, "Two-dimensional continuous elements", European Forum on Aeroelastic and Structural Dynamics, Aachen, 17-19 April 1989, DGLR-Bericht 89-01. [9] M. Levy, "Sur 1,6quilibre 61astique d'une plaque rectangulaire", Comptes Rendus l'Acadkmie Sci. 129, 1899. [10] S. Iguchi, "Die Eigenschwingungen und Klangfiguren der vierseitig freien rechteckigen Platte", Ingenieur Archiv 21, pp. 303-323, 1953. [11] D.J. Gorman, "Free vibration analysis of the completely free rectangular plate by the method of superposition", J. Sound Vib. 57, pp. 437-447, 1978.

114

P.H. Kulla/Finite Elements in Analysis and Design 26 (1997) 97-114

[12] P.H. Kulla, "A solution to the implicit eigenvalue problem", Proc. 2nd Int. Modal Analysis Conf., Orlando, FL, 1984. [13] F. Derksen, "On the transcendental eigenvalue problem", ESA Workshop on Continuum Methods, 15-16 June 1989, ESTEC, Noordwijk, Netherlands. [14] P.H. Kulla, "Dynamic stiffness of rectangular plates and its spectral dcomposition", ESA/ESTEC Contract No, 9562/91/NL/JG(SC), Final Report, March 1993. [15] A.W. Leissa, "Vibration of plates", NASA Report Sp-160, 1969. [16] T.H. Richards and A.Y.T. Leung, "An accurate method in structural vibration analysis", J. Sound Vib. 55, 1977.