$3.00 + 0.00 0045-7949/88 0 1988 Pergamon Pressplc
Compurers & S~ucrures Vol. 29, No. 6. pp. I 119-l 130, 1988 Printed in Great Britain.
PROCREATED
C’ DISPLACEMENT GENERIC VIBRATION MODEL
PLATE
AYSE ALAYLIQCLU TPA, Structures, D-605, P/Bag X197, Pretoria 0001, South Africa (Received
I December 1987)
Abstract-Procreation of the C’ continuous displacement FE efficient vibration model, which can be attained by exploiting the unparalleled conciseness and cost-objectiveness of a computer automated exact analytical integrator, is considered. The proposed Kirchhoff-conforming, rectangular four dof/node plate FE vibration simulator is demonstrated to achieve high accuracy with low CPU expenditure in a straightforward and systematic framework without calling for any ad hoc mechanisms. Due to idiosyncrasies of the formulation, the code has the potential for providing high-speed computational power to aid accurate and reliable identification of a wide spectrum of vibration modes that faithfully captures the physical behavior of thin plates, frequently encountered in real-life engineering applications.
BASIC QUERIES ON CURRENT C’ DISPLACEMENT PLATE CODES Many engineering structures are comprised of predominantly thin plates, and thus simulations of their vibration characteristics necessitate the development of thin plate finite element codes. The conforming thin FE requires C’ continuity, as dictated by the Kirchhoff-Love hypothesis defining the kinematic schemata. This has been however, ab initio, an arduous and expensive number crunching practice dependent upon exact quadrature execution to maintain the intrinsic accuracy/reliability of the C’ continuum representation, since a high number of total degrees of freedom accumulates in the assembled mesh. Although many of the difficulties besetting the implementations of the classical thin plate theory are circumvented by problem discretization decisions based on Co continuity approach, the CO-type thin plate description has often been rather unreliable if not inaccurate for nonacademical, real-life engineering practice. To alleviate this problem, while many important applications await thin plates, some proponents of the FEM have decided to deceive the plate with regard to its thickness. This has involved a mutually motivated procedure with the optimism of creating a dual-purpose plate generator to serve well with Co continuity for every thick/thin plate situation. However, since a physical plate element cannot concurrently be thick and thin, so, not surprisingly, neither can the numerically oriented theoretical plate FE formulation, in order not to defy the long-established plate theories. Ironically, suppression of the well-known theoretical/physical differences in the thick and thin plate behaviors has been tried through the Mindlin [I] thick plate theory for the sake of generality, if not for the high quadrature expenses. Thereby, a slave thin plate behavior has been emulated from its thick plate
master, upon evoking some heuristic numerate integrators which possibly can be manipulated on an ad hoc basis by FE experts for attaining good results under selected scientific circumstances. The related publications of such results have grown at a prodigious rate, albeit a state of confusion and frustration created among the already bewildered average FE users. The thin plate simulation capability of the conventional thick plate FE formalism, ameliorated by encompassing some painstaking heuristic features, such as uniform/selective reduced integration with a priori or a posteriori anti-hourglassing controls, while purporting to be efficient, cannot be taken at face value. With hindsight, the effectiveness of an FE refers to the faithful representation of the overall physical phenomenon rather than the instantaneous numerical performance. Indeed, for the last several years, the FE analysts have investigated the possibility of obtaining accurate solutions from thick plate FE formalism, as a means of significantly improving computational efficiency. These solutions, however, have often exhibited overstiff plate behavior, the so-called locking phenomenon, as the thin plate region is entered. Ad hoc algorithms, although useful to cure locking in some applications, can be costly in terms of analysts’ time and uncertainty about the solution. Such ad hoc codes are not user-friendly and can be used effectively only by FE experts. With this in mind, the state-of-the-art plate problem of today is the rigorous design of accurate and cost-objective elements maintaining spurious mechanisms/locking-free behavior in all situations. In view of the fact that the economic scale of global analysis is manipulating the users to favor relatively simpler elements, a need for efficiently implemented, adequately robust and user-oriented elements is in evidence. This is possible through a combination of
AYSE ALAYLKJGLU
1120
good planning of computer requirements and proper quality assurance of a comprehensive description of an appropriate plate generator. In this setting, it is important to recognize the success of the displacement based FE due to the simplicity of its formulation and to identify the real time spenders in order to make them efficient. As the well-known slogan ‘keep it simple to keep it fast’ sums up the sentiments very well, many efforts in the direction of simplicity have been undertaken. These efforts have been mainly exercised in the form of underintegration, which in turn produces spurious mechanisms and loss of accuracy. The fact remains that the behavior of an FE with full quadrature to maintain accuracy is suboptimal from the simplicity point of view, and can be improved by the unparalleled simplicity of a computerized exact analytical treatment. With the automated analytic integration, the computational speed-up in the context of the costly integrand evaluations is equivalent to the simplicity of one-point Gauss quadrature interaction. The onepoint rule, however, often cannot be used because of its accuracy and stability handicaps (hourglassing; undesirable oscillations due to spurious modes in the kernel of the underintegrated element matrices). The present thin plate vibration capability, therefore, reverts to the Cl continuous, Kirchhoff conforming plate modelling with the simplicity and efficiency provided by exact analytic integral routine invocation, not allowing quadrature or allied ad hoc numerate procedures as a part of design optimization. Elements of this type have apparently not been implemented in other displacement FE codes heretofore. With the present rigorous treatment, it is possible to reflect the genuine behavior of the thin plate structure more closely and comprehensively. The contents of this paper are arranged as follows. In the forthcoming section, the Analytically Integrated Displacement ‘AID’ version of the dynamic FE formulation is presented. Simple expressions for the FE matrices, which can easily be coded in a high level language, viz. Fortran, are provided. This is followed by the rectangular vibration element representation. The fourth section identifies the skeletal computer code that paces the progress of reliable and cost-objective FE modelling. The fifth section describes the benchmark test problems and provides a brief review of the annotated plate FE and a catalogue of their results by identifying their major characteristics. The accuracy profiles of the numerical solutions are reviewed in a tabular format. Conclusive remarks regarding the experiences with the present treatment constitute the penultimate section.
A in terms of the generalized nodal displacements A as
where i refers to the elemental displacement component index, S stands for the nodal displacement component index cumulating over the sequential nodes of the element (the intrinsic element’s degree of freedom index) and S, identifies the assumed displacement shape function which may be represented in a genera1 cast of bivariate form
(2)
Herein, xi/k and Y$,are the shape function coefficients, Gilk and Y,, are the natural coordinate functions @i/k= rl*lfk)
Y,, = ( Ql,
(3)
expressed within the intrinsic plane ~-6 A, and B,, are the exponents, and the indices k and 1 refer to the entity counters within each shape function entry. Since the strain vector Kis related to the displacement vector A through the linear operator matrix L, that is, rci= LjiAi,
(4)
any of the strain components ~~ can be deduced in a format similar to that of Ai in eqns (lH3). Henceforth, the format for uj becomes
with
wherej is the strain component index, Nj,k and Z, are the coefficients, while C,, and Oj,, identify the exponents. Based on the above definitions, equations of dynamic equilibrium can be brought in through the use of Hamilton’s principle for the energy functional II,, within the time interval [t,-t2], requiring 6ll,, = 0, i.e. 6 ‘2 1 ,12Ve S{S[
P
C&&-CC(KLE~~K~ i k
I
(7) DYNAMIC
SCHEMATA DESCRWTION ‘AID’ MODEL
The present dynamic
FOR
displacement finite element ‘AID’ modelling assumes the elemental displacements
Here, Ehi specifies the elasticity tensor, pi and e, are the prescribed forces and strains, V’ is the elemental volume, p is the material density, and & is the
1121
Procreated C’ displacement velocity field vector based on the same shape functions as those of the displacements. Ai = 1 &/A,, /
(8)
whilst /i, is describing the generalized nodal velocities. Upon inserting the definitions (l)--(6) into the condition 6l7,= 0, and performing the required operations for fulfilling the Hamilton’s principle, one obtains
In the above algorithmic representations of the inertial and stiffness matrices and the force vector elements, the requisite integrations can be executed analytically within the computer automated anatomy of these algorithms. This type of symbolic exact integration procedure is simple if one recalls that d V’ = a b t drl dc dc for an element with characteristic lengths a, b and 1, and the natural coordinate variables q, 5 and [ are scaled to vary in the interval [@l]. If eqn (11) is processed for the volume integral of the inertial matrix element Mti according to the definitions of @ and Y in eqn (3), one obtains the following simple integrator-generator algorithm:
imnkl
X
Y;/(Aigm
+
Ajfk
+
1))’ (Big, + Bl/l+ 1))‘.
(15)
The integrator-generator algorithm for the stiffness matrix element Kg, in eqn (12) can be formed via a similar procedure. Henceforth, by consulting eqn (6) for the definitions of r and Q the stiffness matrix element Kg, integrator-generator algorithm is constructed as follows: x pi
dV’dt
=O,
(9)
Kg/ = cbt 2 ccc
and the indicesareg=f, m=n=k=land h=j. Since 6Ar in general can be other than zero, the stationary value of l7” for an element at dynamic equilibrium can be maintained only if
M,,ji, + K,,A,
- F$ - F; = 0,
(10)
cc&,
hjmnk:
I ’
N,,kZj,,(Chgm+C,~k+l)-’
Nkgm
zhgn
tDhgn+Dfl/+
I)-‘.
(16)
Similarly, the integrations for generating the force vector elements in eqns (13) and (14) are performed with the following results:
where c EhjIzj NhgmZ,,,
F? = abt ~~~ h
’ X
Fgn’J’ignX,@,
F,/Yt,,dV’
tchg,n
jmn
+
Dhg,
+
l)-’
(17)
(11)
is the inertial matrix element;
s
Kgf=
“I
x (A,, + B,,, + 11-l.
(18)
CCCCCCNhgmrh.v h
,
m
n
k
/
x Z hgn0 hrnEhl.N.tik r.J/kZ.J/, Q.J/l dl”
(12)
is the stiffness matrix element;
is the force vector element due to cj; and (14)
is the force vector element due to I’,
For computer implementation of the foregoing integrator-generator algorithms depicted in (15)-(18), which consist essentially of a series summation, the simple do-looping technique is used within the related generator-routine structure, and these generators operate only on the non-zero values of the array groups (2, r A, B), (N, Z, C, D), E, I? and P which are declared within the related routines according to their array description format. These algorithms are constructed such that the element generation is independent of its geometry and can be bolted-on to any finite element software. The algorithms are monitored to account for the formation of the lower triangles of the symmetric
1122
AWE ALAYLDGLU
inertial and stiffness matrices within the reduced ranges of the outer do-loops, which sweep through the pair of indices g andS, whilst a matrix element is being generated by performing summations over the nested inner loopings typified by the computing spectrum of the related indices.
computer implementation. In view of the physical boundary conditions, the constraints on the directional changes of W (for restraining 0, and/or 0,) and on the directional changes of 0, or 0, (for restraining K) can be easily implemented and coded for accurate and rapid computer processing. In the above interpolation, the bivariate functions G,, R,, T, and Q, are G,=&‘,Z,,
AID-C'PLATE FE REPRESENTATION Itis the intention here to specify a simple plate FE which can straightforwardly reflect the true thin plate behavior, is amenable to simple and rapid FE processing and can be handled confidently by nonspecialized FE practicers. Simplicity of the FE geometry is one of the major cost reducing factors, and simple elements can potentially be used for a numerically efficient FE discretization of large-scale engineering structures. Accordingly, the plate FE defined here is of a simple rectangle with four corner nodes (Fig. 1) accommodating the C’ continuous field representation by the following interpolation for the flexural displacement W in terms of the degrees of freedom (DOF) at each node r
R,=&‘,Z3
&(q) = &(q) = q - 2q2 + q3 x$(q)
= &4(q) = -q* +q3
T2=X4%,,
Q,=.%$,
G3=H2El,
R3=X2E.,
T3 = X4&,
Q3 = X4&
G4=X,Zz,
R4=X,E4
T4=.X3E2,
Q4=Z3E4,
(22)
(21)
while W= i G,WI,-bR,B,I,-oT,BiI,-~P,~l,, r=l
(19)
while the nodal out-of-plane rotations O,, 0, and twist K (curvature strain component) of the middle surface are defined as
iw, (q)
= S,(q) = 1 - 3q* + 2q-’
x;(q)
= Z*(q) = 3q* - 2q3
X3(q)=&(q)=q X4(q)=Z4(q)=
It is noteworthy that the direct inclusion of the Q,,, 0, and K DOFs in this form is quite effective for the
-2q2+q3 -q2+q3
(22)
denote the first order Hermite polynomials, with q = q for X”s and q = 5 for Es respectively. The stiffness matrix generator algorithm (16) is capable of
Fig. 1. Element topology and notations.
1123
Procreated C’ displacement handling a completely anisotropic material behavior with six independent elastic constants for Ehj as in E,,
[El =
E,z
43
41
E23
1 ’
E33
[ (Symm.)
-I/a*S’O;8,
K,, =
W, + b/a2.?f’P;E30,,,
+ l/a&‘ED;8,Qc, + b/2a &‘yE3K, - l/a**;E,
(23)
W, + b/a23EP;Z3Q,,
+ I/a &:8,Q,,
In the case of orthotropic plates (principle orthotropic directions coinciding with the elemental coordinates), the number of independent constants reduces to four, since E,3 = Eg3 = 0. For numerical performance evaluations of the present algorithm, however, isotropic material behavior is assumed here (to provide a wider range of comparison), retaining only two independent elastic constants; the Young’s modulus E and the Poisson’s ratio v, which define the elasticity matrix elements by
+ b/2a fl;B,
- l/a2&~S2 + l/a *;z,
K*
W3 + b/a*S’~8,0,,
Qc3 + b/2a &‘; z4 ~~
- I/a*XO;E, W, + b/a2S~E4Q44 + l/a~~B2QE4+b/2a~“;~~K,.
(27)
Similarly KC =
l/b*&‘,Z; W, + l/b#,4F;Q,,
-
E,, = E,, = Et3/[12(l -v*)]
+ a/b*.X;%;Q,, 42
=
+ a/2bJf’3E;K,
~4,
- l/b*X;E”; W,+ E33 =
(4,
-
(24)
E,,)P.
+ a/b2ti4E;
With these elements in hand, the stiffness matrix integrator-generator algorithm (16) can be written more specificially as Krf = a& ccc
{E,,N,,,
c
mnkl
’
(&m
+
c,,
+
+
E,z[N,,,Z,,,N2~kZ2,,(C,gm
Z,,,
N,
1)-‘(~,,
x (C2gm +
c,,k
+
E22 N2w
z,,
N2,k z2,, cc,,,,,
+
D2,
+
x tC3gm +
c3,k
+
+
+
X (D2gn
D,,
+
+
@I,
I)-’
+
+
+
X4E;~*
1/b2X2Z; W, + l/b x;Z-;Q,,,
+ a/b2&“8;Qc3
+ a/2b&4%!,‘K3
- l/b*X,Z;W,+
l/bX,E:Q,,,
+a/b2~3~“;Q,,+a/2b~38;K,,
l)-’
c,k
Q,, + a/26
(28)
and
I)-’
K,,(= -2/abX0;3;
W,+2/aX{EGQ,,,
-2/abX”;H;
W2+2/aXiB;Q,,,
N,,,Z,N,,.&,
I)-‘(47,
+
1) - ’ +
-
Z,,
+
D2,
x
l/bS2E;Q,,2
D,,
+
+
c2,k
E3, N,,
Z3,
I)-‘]
+
1) - ’
N3/k Z3[,
- 2/ab %‘;E; W, + 2/aZ;H;Q,, I)-’
P3,
+
D,,,
+
I)-‘},
+ 2/b &‘;B;Q;,
+ c%;E-;K,
(25)
where the shape function coefficients Ns and ZS, and CS and Ds of the strain vector K are related to those of the displacement W through the strain-displacement relation the exponents
a*w = -
{
w9
a*w
2
a*w T
-b*at*' abatlar
1 (26) .
Accordingly, the strains are evaluated explicitly from eqns (19), (21) and (22) as follows:
(29)
In the above representations, single (‘) and double (“) primes stand for the first and the second derivatives of the Hermite polynomials &‘(q)s and Z(<)s as depicted in (22), with respect to the related natural coordinate variables q or 5. The subscript numbers (IA) used for the nodal DOFs ( W, Q,, Q,, K) indicate the node numbers for which these DOFs are associated. The strain components which are explicitly written in eqns (27)-(29) have the same structure as prescribed by algorithm (5) with j = 1, 2, 3 re-
1124
AYSE ALAYLICGLU
ferring to the components K?, K;, K,( respectively. Hereunder, the computer coding for the shape function coefficients and exponents (i.e. Ns, Zs and Cs, Ds respectively) is written in the stiffness matrix generation routine merely for the non-zero entries, and the plate element stiffness matrix is generated according to algorithm (25). For generating the inertial matrix of the plate element, the flexural displacement W, which is specified through eqns (19~(22), is first perceived in conformity with the algorithmic definitions (1)43), with i = 1 referring to the W displacement. Thereby, the shape function coefficients xs and Ys, and the exponent As and Bs are coded into the generator routine for their non-zero entries, and the plate element inertial matrix is constructed according to the general form of the inertial matrix generation algorithm (15). In these stiffness and inertial matrix generators, considerable time saving is achieved by performing the computations once, as this avoids the large number of repetitive numerical evaluations which, in the case of fine meshes, may amount to a sizable proportion of the total running cost.
A(l,
11, l)=B(l,
x(1,11,2)=
11, 1)=2
--a,
A(l, ll,2)=B(l,
r(l, 11,2)=
-2
l1,2)=3
~(1, 12, 1) = ~b/2,
r(l, 12, I) = - I
A(l, 12, 1) = B(1, 12, I) = 2 x(1, 12,2) = -ab/2,
r-(1, 12,2) = 1
A(1, 12,2) = B(1, 12,2) = 3.
With a view to illustrating the computer codes for the strain components given by eqns (27)-(29), the K; component is sampled (j = 2) and coded according to the algorithmic descriptions (5) and (6). Due to space considerations, only the entries related to the fourth node (r = 4) are sampled below, wherein f = l3( I) I6 for the fourth node, k = l(l)3 for the 2 and I = I, 2 for the 8 entries of this node, as N(2, 13, 1) = - l/b2
SKELETAL
COMPUTER
CODING
Z(2, 13, 1) = 6 The computer ready integrator-generator algorithms for inertial and stiffness matrices can be processed in the related generator routines, if the computer codes for the interpolated FE displacement and strain components are available. Computer implementations of the FE variables are sampled hereinafter. The accounts of the W displacement interpolation are present in eqns (19x22). These are genetically implemented through the algorithms (1)<3). For instance, if the entries related to the third node (r = 3) are sampled, the code descriptions can be written as follows, with i = 1 for W, f= 9(1)12 for the third node having a nodal DOF sequence of { W, O,, O,, K}, and k = 1 = 1, 2 since the maximum number of entities in 2 and S of this node is 2. x(1,9,1)= A(l,9,
l)=B(l,9,
x(1,9,2)=
A(l,
A(l,
r(1, 10, 1) = - 1
r(l,lO,2)=
IO, 2) = B(l, 10,2) = 3
x(1,11, l)=a,
D(2, l3,2) = I
N(2, 13, 3) = -2/b2 C(2, 13, 3) = 3 N(2, 14, 1) = l/b Z(2, 14, I) = -2 N(2, l4,2) = -3/b
D(2, l4,2) = 1
C(2, 14, 3) = 3
Y(l, 11, I)=3
N(2, 15, I) = a/b2 Z(2, IS, I) = 6
IO, 1) = B(1, IO, I) = 2
x(1,10,2)=2&
C(2, 13,2) = 2,
N(2, 14,3) = 2/b
-2
A(1,9,2)=B(1,9,2)=3 x(1, IO, 1) = -36,
Z(2, 13,2) = - I2
C(2, l4,2) = 2,
1)=2
Y(l,9,2)=
N(2, 13, 2) = 3/b2
Z(2, 14,2) = 6
Y(l,9,1)=3
I
(30)
C(2, 15, I) = 1 N(2, l&2) = -k/b2 Z(2, 15,2) = - 12
Procreated C’ displacement C(2,15,2)=2,
0(2,15,2)=
1
Case (iii). A square plate clamped at one edge, and the rest are free (CFFF). Case (iv). A square plate with opposite edges simply supported and clamped (SCSC).
N(2, l&3) = a/b2 C(2, 15,3) = 3 N(2, 16, 1) = a/(2b) Z(2, 16, I) = -2 C(2, 16, I) = 1 N(2, 16,2) = --a/b Z(2, 16,2) = 6 C(2, 16,2) = 2,
O(2, 16,2) = I
N(2, 16, 3) = n/(26) C(2, 16, 3) = 3. SOME EXAMPLES FOR ACCURACY OF AID-C’ PLATE
1125
(31) ASSESSMENT
As an outgrowth of the FE activity and because of the need for verifying accuracy by independent element testing, a consensus for code validation efforts is evident. However, a set of &facto tests that can expose element errors and weaknesses in dynamic situations is, at present, not available to include all of the FE requirements which have important effects on vibration capabilities. While standards are sought after as tests of fitness for dynamic situations in the spirit of MacNeal and Harder’s [2] static benchmarks, some typical applications to square plate vibration problems are considered herein for the quantitative estimates of the errors in vibratory capabilities, establishing the pattern of root convergence of the present model. The performance of the present AID-C’ element under various boundary conditions is compared with the acclaimed computational results. In some test problems, the solutions are not available to a large number of digits, nor are they provided for all the practical vibration modes of interest. For practical problems, no one class of elements in vibrating plate analysis has yet been found to be pre-eminent. Convergence studies of the present element are made in comparison with the annotated plate vibration solutions, wherein considerable disagreements appear among the existing results. The proposed AID-C’ plate element is examined through four typical test examples. Test cases
Case (i). A square plate simply supported on all edges (SSSS). Case (ii). A square plate clamped on all edges (CCCC).
The FE representations are assimilated to their respective physical boundary conditions. For instance, the kinematic scheme along the simply supported boundaries accounts for restrained lateral displacement together with its tangential derivative, while the normal and cross derivatives are free. In the case of clamped boundaries all DOFs are restrained. If the symmetricity conditions are imposed (as in the quarter plate analysis), the lateral displacement and its tangential derivative are set free, whereas the normal and cross derivatives are constrained. The plates are modelled using different mesh grids of up to 5 x 5 or 6 x 6, and the resulting eigen-problems are solved for the dominant eigenvalues by using the updated subspace iteration based on the simultaneous inverse vector iteration followed by the Sturm sequence check for verification. The references cited in Table 1 are selected from the vast literature available for illustrating the significant contributions to the plate vibration simulation developments that contain numerical examples to demonstrate the effectiveness of the related computational scheme. Primary information for the catalogue of the elements participated in the comparative convergence analyses is presented in Table 1, where the referred elements are typified by the first named author. The recent publications [3,4] by the present author also propose two hybrid FE models connotated by Dl7-Cl and Dll-C? for plates based on the automated analytical integration technique by using algorithms similar to those of the present AID-C’ model. The D17-C’ element is fully conforming with 17 independent stress parameters and C’ continuous displacements, whilst the Dll-Co element is patch test conforming with 11 independent stress parameters and Co continuous displacements. With a view to enabling a direct comparison on the basis of total DOF per problem, only the results of the D17-C’ high precision element are contrasted in the tables. The results of the numerical studies are depicted in Tables 2-5 in terms of the nondimensional frequency parameters i,, defined as I,, = w,pu2(pt/D)“2,
(32)
where o is the circular frequency and D is the flexural rigidity with Poisson’s ratio v = 0.3. The subscripted indices CLand /I refer to the normal modes of the vibrating plate according to the number of half waves parallel to the respective edges. (i) SSSSplute. The results of the simply supported plate analysis are illustrated in Table 2 for the first five distinct vibration modes. The remarkable performance of the AID-C’ model is seen to exist with
Thick quadrilateral
Thin triangular
12
12
12
12
9
4
4
4
4
3
Disp
Disp
Disp
1181
1191
WI
Pll
Negm
Papenfuss
Tessler
Zienkiewicz
Neale
Disp
nonconforming nonconforming
Thin rectangular
Thin rectangular
nonconforming
isoparametric
nonconforming
plate; constraints
Thin rectangular
Thin semiloof
quadrilateral
Hybr
USI
Martins
Thick
11 parameter
shape
rule
terms
operators
and inertial
exact Gauss
integrations
fields
exact numerical
rule for all terms
(p = I) and manual
factor;
functions
rule
matrices
1171 type displacement
to 16, then 2 x 2 Gauss moment-Birkhoff
DOFs
for shear
[I] plate with shear correction
plate; exact Gauss
Mindlin
plate via beam functions
plate based
stiffness
displacements
(c = 7) to avoid locking;
boundary
multiplier
w.r.t.
with displacement
integration
to reduce
order
solution
dependent
on parametric
terms
plate (HCR):
on shear
shell; reduced
is averaged
integrations
integrations
DOFs
exact analytical
(8.2.2) with abstract
fields; computerized
plates
polynomials
displacement
for vibrating
slope and displacement
continuous
by a shear stiffness
frequency
analyses
Brief description
solution
plate; manual matrix
plate; closed form
involving
balancing
plate; inertial
nonconforming
element
isoparametric
conforming
Thin rectangular
1161
u41
MacNeal
dynamic
manual
W distribution
(ACM);
quintic
plate
plate: energy constraints
Plane rectangular
27
1131
Lindberg
Disp
18
1121
Gupta
Thick (3-D) rectangular with discrete Kirchhoff
Disp
Disp
s 11
Fried
nonconforming
Thin triangular
Hybr
20
UOl
Dungat
isoparametric
Thin rectangular
test
stress and 2nd order
on the patch
plate;
Disp
4
191
Dawe
conforming
Thin triangular
Disp
Disp
]gl
Cowper
nonconforming
Thin rectangular
Disp
16
]71
Clough
plate: 8th order
Thin rectangular series expansions.
24
Hybr
4
PI
Brandt
flat shell based
Thin triangular
18
Disp
Disp
PI
Argyris
stress-C’
in the comparative
plate: 17 parameter
conforming
Thin rectangular
DOF
16
Node
participated
Hybr
Method
of the elements
8
141
D17-C’
1. Catalogue
4
Ref.
Element
Table
integrations
and matrix
integrations
z : r 2 x z:
Table 2. Comparative convergence profiles for frequency parameters ks of an SSSS plate Mesh
Element
11,
221= A,*
1x1
AID-C’ Brandt
20.976 23.349
58.992 71.690
92.564 149.677
2x2
AIDC’ D17-C’ Brandt Cowper Dawe Fried Negm
19.782 19.740 20.163 19.744 18.02 20.100 19.218
52.99 1 51.698 51.963 50.215
83.905 80.490 93.395 80.882
3x3
AID-C’ D17-C’ Brandt
19.748 19.740 19.905
49.737 49.683 51.605
79.486 78.970 82.372
107.68 107.54
135.76 129.27
4x4
AID-C’ D17-C’ Brandt Cowper Dawe Fried Negm
19.742 19.739 19.827 19.739 19.15 19.802 19.537
49.476 49.459 50.533 49.360
79.127 78.958 80.652 78.985
100.186 99.323
129.57 129.02
98.948 96.21
128.66
5x5
AIDC’ D17-C’
19.740 19.739
49.401 49.346
79.027 78.957
99.340 98.711
128.84 128.46
6x6
Cowper Dawe Fried
19.739 19.46 19.764
49.349
78.959
98.718 96.87
128.34
8x8
Argyris MacNeal Tessler
19.778 19.952 19.675
49.496
19.739
49.348
Classical solution [22]
122
A,, = A,,
118.70 117.61
146.94 129.75
100.76 124.42
142.80
124.92
98.933
110.65 99.728 78.957
98.696
128.30
Table 3. Comparative convergence profiles for frequency parameters I, CCCC plate Mesh (a)
A,, = A23
of a
Element
i II
I 31
I 13
133
AID-C’ D17-C’ Clough Dawe Fried Neale Negm
36.057 36.015 34.309 34.31 34.843 36.240 35.488
134.55 133.31
135.05 133.90
227.39 217.54
133.27
134.19
215.30
3x3
AID-C’ D17-C’ Dawe Dungar Neale
36.001 35.993 35.09 36.800 36.129
132.47 132.09 129.70 146.85 132.82
133.07 132.73
221.79 217.95
133.20
221.78
4x4
AID-C’ D17-C’ Clough Dawe Dungar Fried Martins Neale Negm
35.991 35.988 35.454 35.45 36.410 35.755 35.941 36.058 35.823
131.87 131.69
132.49 132.31
220.62 219.30
119.45 132.16
132.66
221.93
AID-C’ D17-C’ Fried Neale
35.987 35.986 35.947 36.026
131.70 131.64
132.32 132.26
220.28 219.54
131.89
132.45
221.36
6x6
AID-C’ Clough Negm
35.986 35.739 35.907
131.64
132.26
220.16
8x8
Argyris
35.983
Classical solution [22]
35.986
131.64
132.24
219.93
2x2
5x5
1127
129.62
130.30 139.27
1128
APSE ALAYLIOGLU Table 4. Comparative convergence profiles for frequency parameters I,, of a CFFF plate Mesh
Element
A,,
1 12
1 21
I 13
1x1
AID-C’ Dl7-C’ Brandt Dawe Gupta
3.518 3.510 3.711 3.32 3.915
8.953 8.894 9.449 8.93 7.841
28.89 27.29 25.90 30.32
37.41 31.07 32.68 35.86
44.24 39.66 36.50 46.39
2x2
AID-C’ Dl7-C’ Brandt Dawe Gupta Martins
3.486 3.481 3.524 3.46 3.732 3.328
8.567 8.557 8.739 8.55 8.136 8.133
21.48 21.46 24.36 21.77
27.41 27.33 29.39 26.33
31.53 31.42 35.13 30.23
19.80
27.24
29.27
3x3
AID-C’ Dl7-C’ Brandt Dawe Gupta Martins
3.478 3.475 3.496 3.47 3.616 3.404
8.528 8.524 8.610 8.53 8.228 8.336
21.37 21.37 22.34 21.67
27.33 27.31 27.84 26.85
31.16 31.12 32.29 30.80
20.17
26.43
30.20
4x4
AID-C’ Dl7-C’ Brandt Gupta Martins
3.475 3.473 3.486 3.561 3.432
8.518 8.516 8.565 8.266 8.403
21.33 21.32 21.85
27.25 27.25 27.55
31.04 31.03 31.66
20.56
26.72
30.63
5x5
AID-C’ Dl7-C’ Gupta Zienkiewicz
3.474 3.473 3.532 3.469
8.514 8.512 8.284 8.535
21.31 21.30
27.22 27.22
31.00 30.99
54.32 54.30
21.45
27.06
3.472
8.525
21.29
27.20
31.10
54.30
Classical solution [23]
42
I 23
66.47
55.80 54.96
54.85 54.67
54.46 54.40
Table 5. Comparative convergence profiles for frequency parameters A,, of an SCSC plate Mesh
Element
2x2
AID-C’ Dl7-C’ Clough Fried Lindberg Negm Papenfuss
1 IL 29.49 29.34 25.83 26.86 29.48 28.14 29.50
3x3
AIDC’ Dl7-C’
29.06 29.03
55.33 54.46
70.73 69.97
96.66 93.76
111.30 109.21
149.25 132.98
4x4
AID<’ Dl7-C’ Clough Fried Lindberg Negm Papenfuss
28.99 28.98 27.74 28.83 28.99 28.61 29.70
54.95 54.58
69.89 69.60
95.31 93.82
103.83 102.38 98.61
131.84 130.99 125.89
102.08 105.74
129.64 132.18
5x5
AID-C’ Dl7-C’
28.97 28.96
54.83 54.66
69.56 6944
94.89 94.15
102.93 102.30
130.74 130.38
6x6
AID-C’ Dl7-C’ Fried Lindberg
28.96 28.95 28.95 28.96
54.79 54.70
6944 69.38
94.73 94.34
102.57 102.24
129.93 129.75
8x8
Clough Negm Papenfuss
28.61 28.84 29.45
100.36 101.68 103.84
127.63 128.80 130.60
102.22
129.10
Classical solution [22] a minimum
number
28.95
121
i 12
1 22
I 31
I I3
58.91 54.06
79.02 76.31
104.54 92.26
122.01 118.29 117.87
161.48 135.50
119.86 122.82
54.74
of mesh grid. For a moderate
mesh size of 5 x 5, the overall response of the element is impressive and is seen to converge to the classical solution [22]. The performance of the proposed AID-C’ element is much better than the pertinent
69.33
94.59
results of the other models, except the author’s earlier element D17-C’, which uses potentially more powerful hybrid technique. The slightly better performance of the more elaborate and higher order triangular displacement model of Cowper [8], that embodies six
Procreated C’ displacement DOF per node, however, contributes
influence for academical in lieu of practical applications. (ii) CCCC plate. The comparative convergence profiles of the frequency parameters for one quadrant analysis of the clamped plate is reported in Table 3 for the first four vibration modes. Due to the more accurate predictions throughout, the present element is advocated to afford greatly improved performance than all of the compared elements, except that of the high precision Dl7-C’ hybrid treatment. (iii) CFFF plate. Table 4 highlights the numerical potentials of the referred elements, in comparison with the accurate and efficient AID-C] model for the first six vibration modes of this plate. The present element is observed to yield better accuracy than the other elements, except in a few instances of Dawe’s [9] element and the slightly more accurate Dl7-C’ model. (iv) SCSC plate. The results given in Table 5 demonstrate the convergence characteristics of the AID-C’ model to the high precision hybrid FE solutions by the Dl7-C’ element for this plate analysis, for which relatively few results exist in the literature. The AID-C’ model is seen to provide better accuracy than the other elements, except in a few instances of Negm’s [18] solution. CLOSURE
A bona Jide thin, C’ plate simulation numerate software for determining the dynamic response of rectangular plates has been proposed. This method does not require prior knowledge of natural frequencies and modal shapes. In general, the method provides results of high accuracy by spending less time than the conventional C’ biased, ad hoc, thick plate FE techniques. The importance of the methodology lies in that it combines the advantages of both the computerized exact analytical integration and the symbolic treatment of the required FE properties. The proposed AID-C’ plate vibration treatment is straightforward and user-oriented in order that it can be administered easily by nonexpert analysts. The practical efficiency of the present FE model, which can be implemented without any difficulty into the libraries of existing FE software, has been recognized with reference to the following desirable element characteristics: (i) cost-objectiveness; (ii) rapid convergence; (iii) freedom from zero energy mechanisms and locking; (iv) invariance with respect to the direction of reference coordinates; (v) suitability to problems with different boundary conditions; and (vi) robustness and reliability. When the proposed AID-C’ plate vibration displacement model is contrasted with the present author’s earlier Dl7-C’ hybrid model, both of which are based on the same displacement functions encompassing C’ continuity and automated exact analytical integration algorithms, the following
1129
observations are made: The Dl7-Cl model provides marginally higher accuracy in vibration mode predictions than the AID-C’ model. However, the AID-C’ model works with marginally lower CPU time than its hybrid FE counterpart. The desirability of these performance characteristics will dictate which of the two models is more suitable for the analyst’s situation; the analyst will then place more weight on that performance to reflect its importance. In consequence, the present AID-C’ displacement model suggests a direction for improving productivity that appears to have significant importance in the context of large-scale dynamic FE investigations. REFERENCES
1. R. D. Mindlin, Influence of rotary inertia and shear on
flexural motions of isotropic, elastic plates. ASME, J. appl. Mech. 18, 31-38 (1951).
2. R. H. MacNeal and R. L. Harder, A proposed standard set of problems to test finite element accuracy. FE Anal. Des. 1, 3-20 (1985). 3. A. Alaylioglu and H. Alaylioglu, Effective hybrid hierarchical element generator for vibration analysis of plates. Int. J. Numer. Meth. Engng 24, 16151628 (1987). 4. A. Alaylioglu and H. Alaylioglu, Hybrid plate vibration models for coarse-mesh analysis. Compur. Struct. 28, 789-796 (1988). 5. J. Argyris, M. Haase, H. P. Mlejnek and P. K. Schmolz. Trunc for shells-an element possibly to the taste of Bruce Irons. Inr. J. Numer. Meth. Engng 22, 93-115 (1986).
6. K. Brandt, Calculation of vibration frequencies by a hybrid element method based on a generalized complementary energy principle. Int. J. Numer. Meth. Engng 11, 231-246 (1977). I. R. W. Clough and J. L. Tocher, Finite element stiffness matrices for analysis of plate bending. Proc. Conf Matrix Meth. Siruci. Mech. DD. 515-545. WriahtPatterson Air Force Base, OH ii966). 8. G. R. Cowper, E. Kosko, G. M. Lindberg and M. D. Olson, Static and dynamic applications of a highprecision triangular plate bending element. AIAA Jnl7, 1957-1965 (1969).
9. D. J. Dawe, A finite element approach to plate vibration problems. J. Mech. Engng Sci. 7, 28-32 (1965). 10. R. Dungar, R. T. Sevem and P. R. Taylor, Vibration of plate and shell structures using triangular finite elements. J. Strain Anal. 2, 73-83 (1967). 11. I. Fried, Residual energy balancing technique in the generation of plate bending finite elements. Compur. Strucr. 4, 771-778 (1974).
12. K. K. Gupta, Development of a finite dynamic element for free vibration analysis of two-dimensional structures. Int. J. Numer. Meth. Engng 12, 13 I l-1327 (1978). 13. G. M. Lindberg and M. D. Olson, Convergence studies ^. ot eigenvalue solutions using two finite plate bending elements. Inf. J. Numer. Mefh. Engng 2, 99-116 (1970). 14. R. H. MacNeal, A simple quadrilateral shell element. Compuf. Sfruc~. 8, 175-183 (1978).
15. R. A. F. Martins and D. R. J. Owen, Thin plate semiloof element for structural analysis-including stability and natural vibrations. Int. J. Numer. Meth. Engng 12, 1667-1676 (1978). 16. B. K. Neale, R. D. Henshell and G. Edwards, Hybrid plate bending elements. J. Sound Vibr. 23, 101-112 (1972). 17. G. Birkhoff and H. L. Garabedien, Smooth surface interpolation. J. Math. Phys. 39, 353-368 (1960).
1130
Avsn ALAYLIOGLU
18. H. M. Negm and E. A. Armanios, Calculation of the natural frequencies and steady state response of thin plates in bending by an improved rectangular element. Comput. Strucf. 17, 139-147 (1983). 19. S. W. Papenfuss, Lateral plate deflection by stiffness matrix methods with application to a marquee. MSc. Thesis, Dept of Civil Engng, University of Washington, Seattle, WA (1959). 20. A. Tessler and T. J. R. Hughes, An improved treatment
of transverse shear in the Mindlin-type four-node quadrilateral element. Comput. Merh. appl. Mech. Engng 39, 311-335 (1983). 21. 0. C. Zienkiewicz, The Finite Element Method, 3rd Edn. McGraw-Hill. Maidenhead. U.K. (1977). 22. A. W. Leissa,‘The free vibration of rectangular plates. J. Sound Vibr. 31, 257-293 (1973). 23. A. W. Leissa, Vibration of Plates. NASA SP-160, U.S. Govt Printing Ofice (1969).