RR CM 0306

CSIR Centre for Mathematical Modelling and Computer Simulation

STUDIES ON VARIATIONAL CORRECTNESS OF FINITE ELEMENT ELASTODYNAMICS OF SOME PLATE ELEMENTS

R MURLIKRISHNA*, G PRATHAP ( * Engineering Analysis Center of Excellence Ltd,Bangalore

Research Report CM 0306

July 2003 Bangalore 560 037, India

STUDIES ON THE VARIATIONAL CORRECTNESS OF FINITE ELEMENT ELASTODYNAMICS OF SOME PLATE ELEMENTS

R. MURALIKRISHNA * and G. PRATHAP**

*Engineering Analysis Centre of Excellence Ltd., Bangalore 560 066, India, email: [email protected]

** CSIR Centre for Mathematical Modeling and Computer Simulation, Bangalore 560 037, India, email: [email protected]

Abstract

The finite element elastodynamics of structural elements like beams, plates and shells have been studied extensively, using both lumped and consistent mass matrix approaches. Recently, the variational bases for the derivation of the element mass matrices have become clearer, and in this paper, the performance of some popular C0 and C1 plate elements is closely examined from this point of view. It becomes clear that the use of non-conforming formulations or the use of lumped mass approaches is an extra -variational one, and introduces errors which are unpredictable in nature. It is also clearly evident that in some cases, the errors due to these extra-variational steps fortuitously compensate for errors in the stiffness matrix, leading one to the erroneous conclusion that a very accurate computational model has been achieved. This has important implications, especially where the aim is to derive error estimators for adaptive

mesh refinement using non-conforming elements and/or lumped mass approaches.

1.0 Introduction

The finite element method is considered to be one of the great innovations in engineering of the last century. Given the kind of wizardry that is now possible due to high performance computing, automated finite element computing based on adaptive optimal meshing is now practicable [1]. To reach this goal, it is necessary for rules to assess the error of computation in a given mesh configuration. An important trend in recent years is to perform a priori or a posteriori error estimation, and use this for adaptive mesh refinement, so that an optimal computational solution is arrived at automatically by the computational process. This approach is beginning to revolutionize the finite element industry in the past one and half decades [1,2]. Commercial finite element programs have started including adaptive mesh refinements as a part of their package.

Most error indicators (or error measures or estimates) are based on an a posteriori evaluation of errors as the computation proceeds [3-5]. These are introduced intuitively, based on an understanding of the energy terms involved, e.g. strain energy for elastostatics problems, strain and kinetic energy for elastodynamics, etc. The end-user of many black-box finite element packages is usually not fully aware of the theory behind the error measures and the strategies behind adaptive mesh refinement. And, often, a clear insight into the definition of the error estimates is wanting. At the same time, a large body of literature has been created by mathematical analysts using the concepts of variational calculus and functional analysis to derive projection theorems and energy error rules on an a priori basis.

These were first reported more than thirty years ago [6]. Unfortunately, these statements on error analysis of the finite element method are derived using rigorous mathematical abstractions which are difficult for the average engineer (who remains the biggest user) to grasp. It is indeed possible to re-derive these using the energy, virtual work and least action principles that the engineer or physicist is more familiar with, especially for linear elastostatics [10]. As an FEA solution is approximate, the errors in the strains and stresses have energies associated with them. It turns out that for linear elastostatics, it is possible to set up simple energy rules that show that a variationally correct FE statement will lead to proje ction rules that ensure that globally, for the problem, and in most cases, locally, in each element, the energy of the error is equal to the error of the energies [6,10]. This measure should be used in any adaptive strategy that tries to obtain an optimal mesh grading for a given mesh. In all these, it is assumed that all operations are based on variationally correct steps, e.g. that the basis functions are complete, and continuous within element and across element edges (continuity or conformity) and that all energy terms are exactly integrated. However, this has not prevented finite element researchers from experimenting with extra-variational procedures like the use of non-conforming elements, the use of reduced or selective integration, the use of lumping of mass terms at nodes to simplify the elastodynamical models, etc. While these variations are used with caution, they are attractive because quite frequently they produce very accurate results. So far, very few investigations have studied carefully, the consequences of the extra -variational procedures in the light of the projection theorems and energy-error rules of the analysis of the finite element method [6]. It is the aim of this work to bring to light these aspects. Though, publications exist which deal with the error measures of plate bending problems [12,13,14,15], and errors in their eigenfrequencies [add L&O here,16,17,18,19], a majority of these publications do not address the points discussed here.

In this paper, we investigate three familiar and widely used plate elements. Accurate structural elements for a variety of beam, plate and shell configurations and problems (elastostatics, elastodynamics, stability, etc.) have been available over the last three to four decades. The C0 and C1 plate elements which are based on the classical Kirchhoff-Love theory come in both the conforming and non-conforming types, and we use the well known ACM element which is nonconforming, and the BFS element, which is conforming in its rectangular form, to evaluate the role that the extra-variational aspect of non-conformity plays in determining the order or rate of convergence [20,21]. It is also useful to have along with these, the very popular C0 plate elements based on Reissner-Mindlin theory, for comparative studies [10]. All three elements are also investigated both from the consistent mass and lumped mass points of view. Our investigations will make clear that whilst for the sake of computational economy, lumped mass matrices are frequently used, adaptive refinement using these lumped matrices should be used with caution. Lumping is seen to be an extra-variational process which does not guarantee the boundedness of the results. The same is the case with non-conforming elements, wherein there is no definite boundedness of the solutions. Any error measures which involve such processes as lumping, non-conforming elements, and other extra-variational techniques are not a good measure of the actual error, and the refinements based on these errors must be used with care. 2.0 A non-conforming element (ACM), and two conforming elements (BFS, QUAD4)

The classical case of a non-conforming element can be quoted with the example of a thin plate bending problem. It has been known for nearly four decades that plate bending problems require C1 continuity. Shape functions insuring this C1

continuity for a rectangular element are easy to formulate, and many such formulations are available in literature, references [20,21]. Some of these elements are conforming (the BFS element) and others non-conforming (the ACM element). The rates of convergence of the results of these elements for both elastostatics and elastodynamics have also been reported. However, all these results were in the late 60s, when the concepts of error estimates and adaptive mesh refinement techniques were not fully established. Now, in the present context, it is a worthwhile effort to have a re-look at the behavior of these classical elements, with respect to their convergence for both elastostatics and elastodynamics. This is done in the light of recent work [22,23,24,25] to examine order of error and rates of convergence in terms of energy principles. It is worthwhile to compare this with the very popular C0 elements based on the Reissner-Mindlin theory. The convergence rates for elastostatics is examined
1 (but reported only for the C elements), and that for elastodynamics studied

carefully, for both beam and plate examples.

The ACM element has 3 degrees of freedom (displacement w, rotation w/ x, rotation w/ y) at each node. Thus we need 12 polynomial terms which would constitute the shape functions of the element. These 12 terms are selected so that we can maintain spatial isotropy. The following 12 terms are used [1, x, y, x2 , y2 , xy, x3 , x 2 y, xy2, y3 , x3 y, xy3 ]

The BFS element is an improvement on the ACM element. It has 4 degrees of freedom (displacement w, rotation w/ x, rotation w/ y, and a cross-twist 2w/ xy) at each node. In this case, we need to take 16 polynomial terms. [1, x, y, x2 , y2 , xy, x3 , x 2 y, xy2 , y3 , x3 y, xy3 , x 2y2 , x 3 y2 ,x 2y3 , x 3y3]

The inclusion of the additional 4 terms, as compared to the polynomial terms of ACM element, lends inter-element continuity up to the rotation degree of freedom. The interpolation function in this case can be expressed in terms of the Hermitian polynomials. The 4-noded element based on the shear-flexible Reissner-Mindlin theory, the QUAD4 is now the acknowledged work-horse of all finite element packages. Here, we adopt the field-consistency approach to remove the locking problems [10]. The element has 3 degrees of freedom (displacement w, rotations ?x and ?y ; but with the flexibility that these rotations are different from the slopes of the mid-surface w/ x and w/ y, respectively, these differences accounting for the transverse shear strains) at each node. Again, we need 12 degrees of freedom at the 4 nodes, but now as the transverse displacements and rotations are independently interpolated, we need only 4 polynomial terms to constitute the shape functions of the element for each of the displacement w, rotations ?x and ?y quantities: [1, x, y, xy ] 3.0 Error estimates and rates of convergence

Let us now review the projection theorem and error-energy rule for elastostatics, using an approach introduced in reference [22]. Following the nomenclature used in Strang and Fix [8], we write the weak form in terms of the energy inner product for the exact solution u to the problem.

a(u,u) = (f,u) a(u,u h) = (f,u h)

(1)

(2)

where a(u,u) is the bilinear symmetric functional and (f,u) is the integral

ufdV

over the system domain V. The vector u describes the degrees of freedom and the bilinear inner product a(u,u) is a measure of the total elastic strain energy, and (f,u) describes the potential to the applied loads. The first two virtual work statements refer to the exact solution of the elastostatic problem. In (1), the trial function and test function are taken as u and the virtual work argument establishes that (1) is truly satisfied only when u is the exact solution at the point of equilibrium. In (2), we take note of the fact that the test function uh (the Ritz or finite element solution) need not be the exact displacement function for the virtual work principle to be true. For convenience, we take this to be the discrete finite element displacement field, as long as it is admissible (that is, satisfies all the geometric boundary conditions). By using uh for both the trial and test function, we get the actual finite element equations, with the right hand side leading to the consistent load vector and the left hand side representing the stiffness matrix. a(u h,u h) = (f,u h)

(3)

This equation will now reflect the error due to the finite element discretisation. We are now in a position to see how the error can arrive at a(u,u h) = a(u h,uh ) e = u - uh can be assessed. Comparing (2) and (3) and noting that the energy inner product is bi-linear, we

and from this we obtain the projection theorem a(u - uh ,u h) = 0

(4)

The finite element solution is therefore seen to be a best-fit or best approximation solution. In most simple linear elastostatics cases, this would imply that the strains or stresses are obtained in a best-fit sense and that there would be points in the element domain where these stresses or strains are very accurately computed (superconvergence).

From the fact that the energy inner product is bi-linear, we can argue that a(u - uh,u-uh) = a(u,u) + a( u h,u h) - 2 a(u,uh ) = a(u,u) - a( u h,u h) - 2[ a(u,u h) - a( uh,u h) ] = a(u,u) - a( u h,u h) - 2[ a(u - uh,u h )] Introducing the result from (4), we get an energy error theorem which can be expressed as a(u - uh ,u-uh) = a(u,u) - a( u h,uh) i.e., Energy of the error = Error of the energy

(5)

This leads to a useful statement that as the left hand side of (5) is always positive definite, a( u h,uh) < a(u,u) (6)

Thus, in a variationally correct approach, the energy inner product of the approximate (Ritz or finite element) solution will always be a lower bound of the exact energy.

Equations (1) to (6) are a restatement of the equations covered by Theorem 1.1 of Strang and Fix [6]. We have started with three virtual work equations (equations (1) to (3)) and this led very easily to a projection theorem (equation (4)), and energy-error theorem (equation (5)), and a lower bound result (equation (6)). From these statements, we can set simple estimates for the error in a finite element discretisation and therefore, for the rate of convergence of the solution for an elastostatic problem. For the ACM and BFS elements, the functions for the displacement field are complete to cubic terms, i.e. [1, x, y, x2, y2, xy, x3, x2y, xy 2, y3] . . This would mean that the bending (including twisting) strains, which are second derivatives of the displacement field are complete to linear terms, i.e. [1, x, y]. It can be shown by the best-fit nature of the computed strains as described by the projection theorem, that the error of the energies, and hence, the theroretical rate of convergence is O(h4), where h is a characteristic element dimension, provided no variational crimes are committed [6]. However, the ACM element does not ensure the C1 continuity expected, and in plate problems, the order of convergence is degraded to O(h2). We shall see this later when we perform the numerical experiments. For the QUAD4 Reissner-Mindlin element, the domina nt energy is the flexural energy. The rotation fields are now complete to linear terms, i.e. [1, x, y]. As the bending (including twisting) strains are first derivatives of the rotation terms, this would mean that strains are accurate to a constant value w ithin elements. We can once again take recourse to the best-fit nature of the solution to argue that the error of the energies, and the rate of convergence are of O(h2). From the foregoing, it is possible to extend these results to a problem in elastodynamics. Now, the kinetic energy of motion enters into the picture, and a classical variational basis for this is already available through the Lagrangian or

Hamiltonian statements, where both potential energy and kinetic energy enter into the functional.

We can write the weak form in terms of the energy inner product for the exact solution u to the problem and replace the loading term f with the inertial force term 2 u. The earlier equations yield a(u,u) = 2.( u, u) a(u,u h) = 2.( u, u h)

(7)

(8)

where is the inertia density of the domain. Again, the first two virtual work statements refer to the exact solution of the elastodynamic problem. In (7), the trial function and test function are taken as u and the virtual work argument establishes that (7) is truly satisfied only when u is the exact eigenfunction and 2 is the exact eigenvalue. In (8), we take note of the fact that the test function uh (the Ritz or finite element solution) need not be the exact displacement function for the virtual work principle to be true, while 2 remains the exact eigenvalue. Again, we take uh to be the discrete finite element displacement field, as long as it is admissible (that is, satisfies all the geometric boundary conditions). By using uh for both the trial and test function, we get the actual finite element equations, with the right hand side leading to the consistent load vector and the left hand side representing the stiffness matrix. a(u h,uh) = (h)2.( uh, u h)

(9)

This equation will now reflect the error due to the finite element discretisation, appearing both in the eigenfunction, and in the eigenvalue. This makes the assessment of the errors a trifle more complicated than in the elastostatics case

10

earlier, as there is the error in the eigenfunction, u - uh as well as the error in the eigenvalue, (h)2 - 2 to be assessed. Comparing (8) and (9) and also noting that the energy inner product is bi-linear, we can arrive at

[a(u,uh) - [a(u h,u h)] - [ 2 .( u, uh) -(h) 2 .( uh, u h)] = 0 or a(u - uh,uh) - (2u - (h) 2 uh, uh ) = 0 (10)

becomes the projection theorem for elastodynamics.

The finite element solution is still seen to be a best-fit or best approximation solution. However, unlike the simple linear elastostatics cases, this would imply that the strains or stresses are only nearly obtained in a best-fit sense. It is no longer simple now to relate the error of the energy to the energy of the error, as was possible for the elastostatics case as in equation (5) earlier. We shall take up the error of the energies. From equation (7) and (9), we have, [a(u,u) - a(u h,uh)] - [2.( u, u) - (h) 2.( uh, u h )] = 0

(11)

We started with three virtual work equations (equations (7) to (9)) and this led very easily to a projection theorem (equation (10)), and the energy-error theorem (equation (11). Also, it is not easy to show that for a variationally correct formulation, the discretised eigenvalue would always be higher than the exact eigenvalue, whereas for the elastostatics problem, we did have an elegant lower bound result (equation (6)). This is because an eigenvalue problem does not give a unique displacement field u or uh . One can side-step this issue by

11

introducing into equation (11), the idea of a normalized generalized mass, where ( u, u) = ( uh, uh) = 1. This gives us, (h)2 - 2 = - [a(u,u) - a(u h,uh)] (12)

What is important to us now is that the error is still governed by the error in the strain energy, and therefore, as long as no variational crimes are committed, the orders of convergence discussed earlier for the elastostatic case should apply.

4.0 Numerical Experiments

In this section, we carefully choose a suite of numerical experiments to highlight that the error estimates projected in the analysis above are confirmed. An inhouse package is developed comprising the ACM, BFS and QUAD4 elements. The QUAD4 element including field -consistency and edge-consistency [ 10] is used. This QUAD4 formulation with field and edge consistency has been used successfully by Mukherjee and Krishnamurthy [15] , and has given excellent results. Numerical experiments are first carried out using the consistent mass matrices for all three element types. This is followed by experiments with lumped mass matrices. Lumping in all cases is done using HRZ lumping scheme [27]. It is important that we clearly define the measure of error to be used. For elastostatics, error can be estimated by comparing a typical deflection obtained by the finite element model, say wh with the analytical value, w. Thus, Error = wh/w - 1 is an appropriate measure. Note that if the model has been formulated in a variationally correct manner, the computed fem displacement will be lower than

12

the exact value, and the error will be negative in sign. Thus where logarithmic plots are made, it may be necessary to use |Error|.

For elastodynamics, where, a variationally correct formulation would have produced higher frequencies, one can define error as Error = (h )2/2 - 1 Again, if the problem is not formulated in a variationally correct way, boundedness is lost, and it may be possible that frequencies are lower than the analytical ones, in which case it may be necessary to use |Error| in the logarithmic plots. 4.1 Some plate elastostatics with ACM, BFS and QUAD4 elements 4.1.1 A beam elastostatic problem with ACM, BFS and QUAD4 elements

The first example chosen, was that of a simply supported beam, modeled as a single strip of plate elements, but with Poisson s ratio, ? = 0. The length of the beam is taken as 1000 mm, with the width and depth both as 10 mm. The Youngs Modulus used is 200000 N/mm 2. The beam is subjected to a uniformly distributed load of 1 unit. For this case, an exact solution for the central deflection is known. The problem has been modeled with the ACM, BFS and QUAD4 elements. All three elements produced stiffer solutions (wh < w). Figure 1 shows that the ACM and BFS element models, which are complete to cubic order, show identically, the h4 rate of convergence predicted theoretically. Here, in this example, the results for the nodal displacements are reported at the centre of the beam. For this example, the BFS and ACM models give exact displacements at the nodes. Therefore, f or the BFS and ACM element models, the beam was discretised into 1x1, 3x1, 5x1, 7x1, 9x1 and so on. Discretisation in this fashion ensures that no node comes at the centre of the beam, and the displacement at

13

the centre of the beam can be re-computed from the finite element displacements of the element on which the centre of the beam falls. This strategy enables us to arrive at the errors in the displacement of the centre of the beam. Again, as predicted, the QUAD4 element had h2 convergence. In this simple onedimensional example, where the exact solution is based on the Euler-Bernoulli theory, there is no twist term. The need for distinguish the two elements from this example. C 1 continuity is therefore

automatically met by both the ACM and BFS elements and so we could not

Variation of Error in Deflection of a simply supported beam

0 0 -0.5 -1 -1.5 0.5 1 1.5 2 2.5 3 3.5

log(|error|)

-2 -2.5 -3 -3.5 -4 -4.5 -5 log (h)

Mindlin Element O(h^2) BFS/ACM Element O(h^4)

Fig.1 Convergence of Displacement of a simply supported beam

4.1.2 A plate elastostatic problem with ACM, BFS and QUAD4 elements We now turn to a two -dimensional problem where the need for C1 continuity is important, especially as we anticipate that the ACM element does not have the necessary twist degree of freedom to ensure this degree of conformity. A square plate problem is chosen with dimensions 40mm x 40mm, 0.4 mm thick, and with elastic properties, =0.3, E=200000 N/mm2. Both simply-supported and clamped edge conditions were investigated and in turn, the point-load at centre and uniformly distributed loading cases were studied. For these cases, very accurate

14

solutions for the central deflection are known. The plate was divided into elements of equal size, and meshes ranged from 1x1 to 16x16 for one quarter of the plate. In all cases, the deflections from the BFS element model converged from below ( wh < w), showing that the condition expressed in Equation (6) was met, i.e. the energy inner product of the approximate solution will always be a lower bound to the exact energy, and hence deflections would converge from below. The results for the BFS element model showed h4 convergence, as expected. However, in all these cases, the ACM element showed larger deflections than predicted by theory ( wh > w). That is, the condition expressed in Equation (6) was violated, and the only explanation for this is that the nonconforming nature of the formulation has introduced extra-variational conditions. Also, the ACM element showed only a degraded h2 convergence. Figure 2 shows the convergence for the case of the simply-supported square plate under uniformly distributed loading conditions. Note that now, due to the difference in signs of the error (BFS and ACM give errors of different signs), the modulus of the error measure is used for graphical representation.

Variation of log ([e]) vs log(h), for SS plate with uniform load

0 0 -1 0.2 0.4 0.6 0.8 1 1.2 1.4

-2

-log(|e|)

ACM Element -3 BFS Element Order(h^4) Order(h^2) -4

-5

-6 log (h)

e = (wfem/wtheory ) - 1

Fig. 2 Effect of conformity (BFS) and non-conformity (ACM) on rates of convergence

15

4.2 A plate elastodynamic problem with ACM and BFS elements

We now turn to the problem of elastodynamics, which is the central theme of the present investigation. We shall use these elements for the study of the free vibration response of a simply sup ported plate. A plate of dimensions 40mm x 40mm, 0.4 mm thick, is chosen and with elastic properties, =0.3, E=200000 N/mm2, and mass density =7850 kg/m3. To study the rate of convergence of the errors, we follow an approach introduced in reference 18. We start with an accurate high density mesh of uniformly sized elements (NxN mesh) and then compute the error in the frequencies with increasing mode number (for a plate, m or n, varying one and keeping the other fixed). Consider a simply supported plate, where the exact solution involves simple trigonometric waves; for example, the (m,n) mode will be the eigenfunction sin(mpx/a)sin(npy/b). The ratio of the element length, h = a/N to the wave-length of the mth wave ? = a/m, is now r = h/? = m/N, when the other mode number (say n) is kept fixed. The errors will now be a function of r, and this will mean that the order of convergence with one of the mode numbers of the frequency (say m) will follow the exact trend as with h, if N and n are fixed for the problem. Thus, with one single eigenvalue

computation, we can sweep for the errors for varying mode numbers m and n.

Figure 3 shows the error trends in this fashion, for the BFS model of the free vibration of a simply-supported plate, when a consistent mass matrix is used. In Figure 3, we show only the error map for modes with wave numbers m > n. It is seen that the convergence is nearly consistent with the prediction that the eigenfrequencies will have a convergence of O(h4 ). However, what we

actually see when the full range (m,n) is surveyed is a complex map of errors with wave numbers m and n. To get a clearer picture, let us sweep for the errors keeping one mode number, say n a constant, while varying m. Figure 4 shows the errors in natural frequencies of a simply supported plate with consistent mass BFS elements with varying mode number m, for fixed mode number n = 5. It is seen that for m Y n, the errors are a minimum. This would

16

seem to suggest that errors are a minimum when the twisting effects in the plate are a minimum, and are compounded when twisting is significant, i.e. m @ n. Let us now examine the case of minimum twist, where the modes have m = n. Figure 5 depicts the errors in natural frequencies of the simply supported plate with BFS elements with varying mode numbers, m = n. It is seen that the convergence is almost exactly consistent with the prediction that the eigenfrequencies will have a convergence of O(h4). This was also the rate of convergence seen for the elastostatic one -dimensional beam example seen earlier, showing that where there is predominantly bending action, the optimal order of convergence is realized.

This is not the case when we go to the ACM model of the free vibration of the simply-supported plate, when a consistent mass matrix is used. Figure 6 shows the errors in the natural frequencies, when the ACM model is used. It first becomes apparent that the frequencies h computed by the ACM model are lower than the analytically predicted . This leads to the conclusion that unlike

the BFS model, which was a variationally correct formulation and computed higher frequencies, here the problem is no longer formulated in a variationally correct way, and boundedness is lost, with frequencies lower than the analytical ones. In fact, in Figure 6, it became necessary to use E | rror | in the logarithmic plots. Boundedness is lost because the lack of strict enforcement of the C1 continuity required across element edges obviously leads to an error in the tota l strain energy, and hence an error in the stiffness matrix. This leads to greater flexibility as seen by frequencies lower than the analytically correct ones. A consequence of this is that no definite order of convergence is seen in Figure 6. In fact, a mode where m = n, which is a predominantly bending mode, is not the most accurate one from the computational point of view. It is also very interesting to note that for m=n modes, the order of convergence has now degraded to O(h2), whereas it was O(h4) for the one-dimensional elastostatic case of the beam studied earlier. This seems to confirm that the loss of continuity

17

related to the twist term is the culprit that accounts for this degraded performance of the ACM element.

Variation of Log(Error) with Log(mode number) for a SS Plate

0 0 -1 n=1 n=2 -3 n=3 n=4 n=5 -4 O(h^4) 0.2 0.4 0.6 0.8 1 1.2

Log(Error)

-2

-5

-6 Log(m)

BFS, Consistent Mass

Fig. 3 Errors in natural frequencies of a simply supported plate with BFS elements

18

Variation of Log(Error) with Log(m) for a SS Plate, n=5

0 0 0.2 0.4 0.6 0.8 1 1.2

-0.5

-1

BFS, Consistent Mass

-1.5

Log(Error)

-2

n=5

-2.5

-3

-3.5

-4

Log(m)

Fig. 4 Errors in natural frequencies of a simply supported plate with BFS elements with varying mode number m, for fixed mode number n = 5

19

Variation of Log(Error) with Log(m) for a SS Plate, m=n

0 0 -1 0.2 0.4 0.6 0.8 1 1.2

BFS, Consistent Mass

-2

Log(Error)

-3

-4

-5

m=n
-6

O(h^4)
-7

Log(m)

Fig. 5 Errors in natural frequencies of a simply supported plate with BFS elements with varying mode numbers, m = n

20

Variation of Log(|Error|) with Log(mode number) for a SS Plate

0 0 -0.5 n=1 n=2 n=3 n=4 n=5 O(h^2) m=n modes 0.2 0.4 0.6 0.8 1 1.2

Log(|Error|)

-1

-1.5

-2

-2.5

-3 Log(m)

ACM Element, Consistent Mass

Fig. 6 Errors in natural frequencies of a simply supported plate with ACM elements

21

4.3 A plate elastodynamic example to show the boundedness aspects with with ACM and BFS elements

In this section, we examine carefully the question of boundedness of the eigenvalue predictions when conformity is not assured, and the mesh is not uniformly spaced. For this, we choose again, the simply supported plate but use only a 2x2 mesh. Also, we sweep the central nodal point (x, y) of the mesh of the plate, as shown in Figure 7 below, where x is the distance along the length and y is the distance along the width of the plate.

x Figure 7. A 2x2 mesh of the simply supported plate, with the central node (x,y) swept over the plate.

Figure 8 shows the eigenvalues obtained when the node is swept over the plate, so that we can have a highly-asymmetrical mesh, for the BFS element model.

22

Fig. 8 Fundamental Frequency of a Simply Supported Plate (BFS Element), as the central node in a 2x2 mesh is swept

Figure 8 shows that the frequencies are very neatly bounded from below, and the perfectly symmetrical mesh, which is also the mesh which is uniform, has the most accurate frequency. This is not unexpected, as it is known that globally, errors are minimized when the errors are equi-partitioned, and this is the case when the mesh is uniformly laid out. In Figure 9, we show what happens when the central node is swept along the diagonal (x = y). The boundedness aspect is seen very clearly.

23

Variation of Fundamental Frequency of a SS plate 2x2 Mesh - Full plate

247 246 245

Frequency (Hz)

244 243 242 241 240 239 238 237 0 5 10 15 20 25 30 35 40 Distance along diagonal of plate BFS-Consistent Exact Solution

Fig. 9 Fundamental Frequency of a Simply Supported Plate (BFS Element), when the node is swept along the diagonal

24

Fig.10 Fundamental Frequency of a Simply Supported Plate (ACM Element)

Variation of Fundamental Frequency of a SS Plate- Full plate 2x2 Mesh

240 235 230 225 220 215 210 205 0 10 20 30 40

Fundamental Frequency

ACM-Consistent Exact Solution

Distance along diagonal

25

Fig.11 Fundamental Frequency of a Simply Supported Plate (ACM Element), when the node is swept along the diagonal

Figures 10 and 11 show this exercise repeated with the ACM elements. Boundedness is now lost, and in fact, there are points where the mesh, if sufficiently non-symmetrical, gives the exact frequency. This can be interpreted to mean that the extra-variational errors introduced due to loss of conformity have exactly compensated for the errors due to the discretisation process described by Equations (10) to (12), along with the errors due to non-uniformity of mesh. It is easy to see that this has important implications where error estimates are used to achieve automatic mesh refinement. Thus, where compensation of errors is expected, it is possible that a bad mesh chosen deliberately, or unintentionally, can produce a solution with no error!

4.4 Plate elastodynamic examples to show the extra-variational aspect of lumping with ACM and BFS elements

So far, we have examined only formulations where consistent mass matrices were used. We next turn to discuss the effects of lumping of the mass matrices. The lumping scheme suggested in reference [21] is used in the present case. We first discuss the beam problem. Let us take the case of a simply supported beam (1000 mm length, 10mm width, 10mm depth, E=200000 N/mm2, = 7850 kg/m3. When we discretize the beam into two elements, and find the location of the middle node which makes the error in the fundamental frequency the least, we find that the middle node should be positioned right at the centre of the beam, for a model with consistent BFS elements. Also, all frequencies are bounded from below, consistent with what was predicted.

26

Variation of Fundamental Frequency w.r.t position of node

6.5

Frequency

5.5 5

BFS-Consistent Mass Theory BFS Lumped Mass (HRZ)

4.5

4 0 200 400 600 800 1000

Position of second node

Fig. 12 Fundamental Frequency of a simply supported beam effect of lumping

The same problem when repeated with lumped mass matrices, shows that boundedness is lost. In fact, it leads to the erroneous conclusion that the location of the middle node which makes the error in the fundamental frequency the least is not at the centre of the beam, but at either of the two locations shown in Fig. 12.

We now study the problem of the free vibrations of a simply supported plate , using the lumped mass approach (Figure 13). Again, a 16x16 uniform mesh is used, and error maps are drawn for varying m and n. First, we see that errors have changed sign, necessitating the use of kinetic energy, and this error is considerable. |Error| in the logarithmic plot. Lumping is therefore an extra-variational step that has introduced errors in the

27

Variation of Log(|Error|) with Log(mode number) for a SS Plate

0 0 -0.5 -1 n=1 n=2 -1.5 -2 -2.5 -3 -3.5 Log(m) n=3 n=4 n=5 O(h^2) m=n modes 0.2 0.4 0.6 0.8 1 1.2

Log(|Error|)

Lumped Mass,BFS

Fig. 13 Natural frequencies of a simply supported plate effect of lumping(BFS)

Note the rates of convergence have dropped to O(h2) as compared to the O(h4) achieved for the consistent mass case, reported in Fig. 3.

Variation of Log(|Error|) with Log(mode number) for a SS Plate

0 0 -0.5 n=1 n=2 n=3 n=4 n=5 O(h) m=n modes -2.5 O(h^2) 0.2 0.4 0.6 0.8 1 1.2

Log(|Error|)

-1

-1.5 -2

-3 Log(m)

ACM Element,Lumped

Fig. 14 Natural frequencies of a simply supported plate effect of lumping(ACM)

In Fig. 14, we show the effect of lumping on the natural frequencies for the ACM formulation. Here again, we observe the loss of rate of convergence. For m = n,

28

the rate is O(h2), but when there is considerable twist, e.g. for fixed n, and varying m, we see the rate of convergence is as poor as O(h)!

Fig.15 Fundamental Frequency of a Simply Supported Plate (BFS Element) effect of lumping

Fig. 15, emphasizes the boundedness aspect further from the lumping point of view. We use a 2x2 element mesh, with the central node being swept to cover the whole plate. Lumping causes the boundedness to be lost. The perfectly symmetrical mesh yields a frequency that is lower than the exact frequency. For a non-symmetrical mesh, the frequencies increase at first, and at some fortuitous configurations, the exact frequencies are obtained. At highly non-symmetrical situations, the trend changes and the frequencies begin to fall. The main conclusion from this is that it is unreliable to use lumped mass approaches in an automatic adaptive mesh refining approach, as the errors of lumping may actually be compensating for the errors of discretisation. Figure 16 shows what happens when the central node is swept along the diagonal of the plate.

29

Variation of Fundamental Frequency along Diagonal

360 340 320

Frequency (Hz)

300 280 260 240 220 200 0 10 20 30 40 BFS-Lumped Exact Solution

Distance along diagonal

Fig. 16 Fundamental Frequency of a simply supported plate effect of lumping with BFS elements

Figures 17 and 18 show the results obtained when the same exercise is repeated with lumped ACM elements.

Fig.17 Fundamental Frequency of a Simply Supported Plate (ACM Element) effect of lumping

30

Variation of Fundamental Frequency of a SS Plate- Full plate 2x2 Mesh

350 300 250 200 150 100 50 0 0 10 20 30 40

Fundamental Frequency

ACM-Lumped Exact Solution

Distance along diagonal

Fig.18 Fundamental Frequency of a simply supported plate effect of lumping with ACM elements

4.5 Plate elastodynamic examples with QUAD4 elements

We now see the results of the QUAD4 element for the errors in the frequencies of the simply supported plate.

Variation of Log(Error) with Log(mode number) for a SS Plate

0 0 -0.5 n=1 -1 n=2 n=3 -1.5 n=4 n=5 -2 O(h^2) m=n modes -2.5 0.2 0.4 0.6 0.8 1 1.2

-3 Log(m)

C0, Consistent Mass

Fig.19 Errors in Natural Frequencies of a Simply Supported Plate (QUAD4 Element)

31

Fig. 19 shows the error map for the simply supported plate when a uniform 16x16 mesh of QUAD4 elements is used with a consistent mass matrix. The results are similar to that obtained for the conforming C1 BFS element model, except that the order of convergence is nearly O(h 2 ). We can see next, if twist action alters the accuracy. Figure 20 shows how, for a fixed n, the errors vary as m varies. Again, we see that when m = n, the errors are least. In fact, in Fig. 21, we see that when m = n, the rate of convergence is almost exactly O(h2 ).

Variation of Log(Error) with Log(m) for a SS Plate, n=5

0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

-0.2

QUAD4, Consistent Mass

-0.4

Log(Error)

-0.6

n=5

-0.8

-1

-1.2

Log(m)

Fig. 20 Errors in natural frequencies of a simply supported plate with QUAD4 elements with varying mode number m, for fixed mode number n = 5

32

Fig.21 shows the boundedness of the frequency surface. The results of the diagonal sweep test are shown in Fig.22. The results for the QUAD4 element when lumped mass matrices are used are shown in Fig.23 and Fig.24.

Fig.21 Fundamental Frequency of a Simply Supported Plate (QUAD4 Element)

33

Variation of Fundamental Frequency of a Simply Supported Plate

QUAD4-Consistent Mass Theory 1000 800 600 400 200 0 0 5 10 15 20 25 30 35 40

Distance Along Diagonal

Fig.22 Fundamental Frequency of a Simply Supported Plate (QUAD4 Element)

Fig.23 Fundamental Frequency of a Simply Supported Plate (QUAD4 Element)

34

Variation of Fundamental Frequency of a Simply Supported Plate

Theory QUAD4-Lumped Mass 600 500 400 300 200 100 0 0 5 10 15 20 25 30 35 40

Distance Along Diagonal

Fig.24 Fundamental Frequency of a Simply Supported Plate (QUAD4 Element)

5.0 Conclusions :

A detailed study of the errors produced due to the use of extra-variational processes in the plate bending problems has been done. It is shown that whenever the variational principles are not adhered to, there is no guaranteed boundedness of the results. This lacuna needs to be seriously looked into before such error measures are used in the mesh refinement.

References

1. J.Mackerle, Error Estimates and adaptive finite element methods: A bibliography (1990-2000), Engineering Computations, 18, 5-6, (2001),802914. 2. L.Y.Li and P.Bettess, Adaptive finite element methods: A review, Applied Mechanics Review, Vol.50, No.10, pp.581-591, 1997.

35

3. I.Babuska and W. C. Rheinboldt, `A-posteriori error estimates for the finite element method', Internat. J. Numer. Methods Engrg., 12, (1978) 15971615. 4. O.C. Zienkiewicz and J.Z. Zhu "A Simple Error Estimator and Adaptive Procedure for Practical Engineering Analysis", Int. J. Num. Meth. Engineering, Vol.24, 337-357, 1987. 5. M.Ainsworth, and J.T.Oden, A posteriori Error Estimation in Finite Element Analysis, Comput. Methods Appl. Mech. Engrg., 142 (1997), 1-88. 6. G Strang, G J Fix, An Analysis of the Finite Element Method , Prentice-Hall Series in Automatic Computation, Prentice-Hall, Englewood Cliffs, NJ, 1973. 7. B M Fraejis de Veubeke, Displacement and equilibrium models in the finite element method, in Stress Analysis, ed. by O. C. Zienkiewicz and G. S. Hollister, Wiley, London, 145 -197, 1965, reprinted in Int. J. Num. Meth. Engrg ., 52, 287-342, 2001. 8. H C Hu, On some variational methods in the theory of elasticity and plasticity, Scientia Sinica , 4, 33-54, 1955. 9. K Washizu, On the variational principles of elasticity and plasticity, Aeroelastic and Structures Research Laboratories, Technical Report 2518, MIT, Cambridge, 1955. 10. G Prathap, The Finite Element Method in Structural Mechanics, Kluwer Academic Press, Dordrecht, 1993. 11. G Prathap, Barlow points and Gauss points and the aliasing and best fit paradigms, Computers and Structures , 58, 321-325, 1996. 12. J.S.Sandhu and H.Liebowitz, Mesh adaptation using a four-noded quadrilateral plate bending element, Engineering Fracture Mechanics, 50, 5/6, (1995), 737-758. 13. O.C. Zienkiewicz and J.Z. Zhu. Error estimates and adaptive refinement for plate bending problems, Int. J. Numer .Meth. Eng., 28, p. 2839-3854, 1989.

36

14. C.K.Lee and R.E.Hobbs, Automatic adaptive refinement for plate bending problems using Reissner-Mindlin plate bending elements, Int. J. Numer .Meth. Eng., 41, pp.1-63,1998 15. S. Mukherjee,S., and C.S. Krishnamoorthy, Adaptive FE analysis of plates by shear-flexible quadrilateral Reissner-Mindlin Elements, Finite Elements in Analysis and Design, 22, 329-366, 1996 16. C.Zhao and G.P.Steven, A-posteriori error estimator/corrector for natural frequencies of thin plate vibration problems, Computers and Structures, 59,5,(1996), 949-963. 17. D.B.Stephen and G.P.Steven, Error Estimation for natural frequency analysis using plate elements, Engineering Computations, 14,1, (1997), 119-134 18. F.J.Fuenmayor, J.L.Restrepo, J.E.Taranco ,L.Baeza, Error estimation and h -adaptive refinement in the analysis of natural frequencies, Finite Elements in Analysis and Design 38 (2001)137. 19. P.Hager, N.E. Wiberg, Error estimation and h-adaptivity for eigenfrequency analysis of plates in bending: numerical results, Computers and Structures 78 (2000) 1 -10. 20. F.K.Bogner, R.L.Fox. and L.A.Schmidt., The Generation of inter-element compatible stiffness and mass matrices by the use of interpolation formulas, in Conf. On Matrix Methods in Structural Mechanics, Wright Patterson Air Force Base, Ohio, October 1965.

21. R.J.Melosh, Basis for Derivation of Matrices for the Direct Stiffness Method, AIAA Journal, 1963, Vol.1, No.7, 1631-1637 22. G.Prathap and S.Mukherjee, The engineer grapples with Strang and Fixs Theorem 1.1 and Lemma 6.3, Submitted for Publication

23. G.Prathap and D.V.T.G.Pavan Kumar, Error Analysis of Timoshenko Beam Finite Element Dynamic Models, Int.J.of Computational Engineering Science, 2,1,(2001),1-10. 24. K M Liew and S Rajendran, New superconvergent points of the 8-node serendipity plane element for patch recovery, Int. J. Num. Meth. Engrg., 54, 1103-1130, 2002.

37

25. Tessler A. and Hughes T.J.R. An improved treatment of transverse shear in the Mindlin type four node quadrilateral element, Computational Methods in Applied Mechanics and Engineering, 39, 311-335, 1983. 26. G. Prathap., and B.R.Somashekar., Field- and Edge- Consistency Synthesis of a 4-noded quadrilateral plate bending element, International Journal for Numerical Methods in Engineering, 26, 1993,1692-1708. 27. E. Hinton., T.Rock and O.C.Zienkiewicz, A note on mass lumping and related processes in the finite element method, Earthquake Engrg. Struct. Dynamics, 4, (1976),245249

Add Lindberg & Olsons paper

38