- Email: [email protected]
Dynamic response of plate systems by combining finite differences, finite elements and laplace transform
004-7949 $3.00 + .oo Rrgamon press Ltd.
Conrpruers& Suucmres Vol. 19, No. 5/6, pp. 163-115, 1984 Printedin the U.S.A.
DYNAMIC RESPONSE OF PLATE SYSTEMS BY COMBINING FINITE DIFFERENCES, FINITE ELEMENTS AND LAPLACE TRANSFORM
D. E. BE~KOS~ Department of Civil Engineering, University of Patras, Patras, Greece and K. L. LEUNG$ Department of Civil and Mineral Engineering, University of Minnesota, Minneapolis, MN 55455, U.S.A. (Received 15 September 1983; received for publication 5 December
1983)
Abstract-A numerical method for the determination of the dynamic response of large rectangular plates or plate systems to lateral loads is proposed. The method is a combination of the finite difference method, the finite element method and the Laplace transform with respect to time. The plate system is considered as an assemblage of a small number of big rectangular superelements whose stiffness matrices are derived with the aid of the finite difference method in the Laplace transform domain. These superelements are then used to formulate and solve the problem by the finite element method in the transformed domain. The dynamic response is finally obtained by a numerical inversion of the transformed solution. External viscous or internal viscoelastic damping as well as the elastic foundation interaction effect can easily be taken into account. The method is illustrated and its merits demonstrated by means of numerical examples.
INTRODUCIION
The determination of the dynamic response of rectangular plate systems encountered in ship, aerospace and civil engineering structures is an important technological problem. The complex geometry and boundary and loading conditions of single large rectangular plates or plate systems such as one or two-way continuous plates on rigid or flexible supports do not permit an analytic solution of the problem but in very few simple cases[l-31. Static and dynamic analysis of complex plate systems is usually accomplished by numerical methods such as the Finite Difference Method (FDM) and the Finite Element Method (FEM)[3-51. However, both of these methods require a large computer memory for problems dealing with large plate systems due to the fact that a refined discretization is required for acceptable results. Besides, round-off error accumulation during the solution of large systems of equations usually leads to unsatisfactory results. In an effort to overcome these difficulties various special techniques have been employed in conjunction with the FEM for the analysis of plate systems such as the substructuring technique[6], the finite strip method[4, 71 and the finite panel method[8-lo]. With the exception of reference[l, Refs. [&lo] dealt with the static problem only. The finite strip and finite panel methods require a much smaller total number TProfessor and Visiting Professor at the Department of Civil and Mineral Engineering, University of Minnesota, Minneapolis, MN 55455, U.S.A. SGraduate Student.
of degrees of freedom for the whole system but they employ rather complicated displacement functions. It is apparent that a finite difference gridwork for analyzing a plate corresponds to a number of degrees of freedom which is much smaller than that for a finite element mesh with as many nodes as the gridwork points in the plate. On the other hand it is well known that treatment of the boundary conditions by the FEM is much simpler and easier than by the FDM. These two facts suggest that a combination of the two methods in such a way as to exploit their advantages and suppress their disadvantages might offer a solution to the problem of analyzing large plate systems. Indeed this is exactly what Ghali and Bathe[l 1, 121 did for the static analysis of plate systems under in-plane and lateral loading. In this work, the determination of the flexural dynamic response of large rectangular plates or plate systems to lateral loads is accomplished by a combination of the FDM, the FEM and the Laplace transform with respect to time. It is actually an extension to the dynamic case of the work of Ghali and Bathe[l2] with the aid of Laplace transform, which effectively reduces the dynamic problem to a static one. According to the proposed method the plate system is idealized as an assemblage of a small number of big rectangular superelements whose stiffness matrices are derived by the FDM on the basis of the Laplace transformed equation of motion. Use of these superelements in conjunction with the FEM leads to the formulation and solution of the problem in terms of nodal displacements in the transformed domain. Interior displacements and
763
764
D. E. BE~KOSand K. L. LEUNG
stresses for a superelement can then be obtained through relations connecting the boundary nodal displacements with the interior ones. The structural response can finally be obtained by a numerical inversion of the transformed solution. Use of Laplace transform has the advantages that the dynamic problem is effectively reduced to a static one, no need of natural frequencies and modal shapes is required as in modal analysis and other effects such as external viscous or internal viscoelastic damping as well as the foundation interaction can be easily taken into account. Dynamic forces of any time variation can be handled without any difficulty by employing a direct numerical Laplace transform. A comprehensive study on the direct and inverse numerical Laplace transform as applied to complex time-dependent linear problems of engineering science can be found in Narayanan and Beskos [ 131, while various applications of the Laplace transform to structural and soil dynamics problems can be found in Beskos et al.[14-231. The proposed method is illustrated by means of numerical examples and its merits are carefully assessed.
modulus of elasticity and Poi_sson’s ratio of the plate. The Laplace transform f(x, y, S) of a function f(x, y, t) with respect to time is defined by
fix,Y, s) =
sm
e-s’f(w,~) dr,
(3)
0
where s, in general a complex number, is the Laplace transform parameter. Application of the Laplace transform with respect to time on eqn (1) under zero initial conditions yields 4
g+2
a% a% _+_+,w=P. ax2ay2 ay4
-
d
N
(4)
It is very easy to see[24] that the force-displacement relations for the plate element of Fig. 1 in the Laplace transformed domain become
STIFFNESS MATRIX OF RECTANGULARPLATE SUPERELEMENT
Consider a rectangular plate bending element under distributed lateral load p(x, y, t) as shown in Fig. 1. The equation of its lateral flexural motion according to the linear theory of Kirchhoff-Love takes the form[l, 31 2+2-
adw saw M azw ax2ay2 +v+;dr”=$
p (1)
where w = w(x, y, t) is the lateral deflection, m is the plate mass per unit surface, t indicates time and the bending rigidity N is given by
where the moments k,, “y! &, and shear forces v, and v, are depicted in the positive sense in Fig. 1. Consider now a large rectangular plate bending element in the Laplace transformed domain with a typical finite difference gridwork as shown in Fig. 2. On the basis of that gridwork the finite difference form of eqn (4) takes ihe form[24] A,@,
+ A2(Gi+
I,j + wi_ Ij) + Aj(Gi,j+
1 + W,j_ 1)
+A4(~i+,,j+,+Wi+I,j-I+~i-I,j+I+~i-I,j-1)
N = Eh’/12(1 - v2),
(2) +
A.5(w&j+2
+
fi,j_2)
+
Ae(Wi+*,j
+ Wf_2J) = 0,
with h being the plate thickness and E and v being the
Fig. 1. A typical plate bending element with its corresponding forces.
(7)
lateral load and internal
Dynamic response of plate systems
765
and the geometric relationship
where
{a*} = [G]{G*}, A, =6+6a2+8a A,= -4a(l
+/?, A,= -4(l
+a), A,=2a,
As=a2,
a = (~,/~J*, j3 = (lQN)mS*.
A,= 1, (8)
Equation (7) is valid for an interior point (i,j) such as that of Fig. 2. For other points (i,j) close to or on the boundaries of the plate element, which are assumed to be free, eqn (4) can be written in finite difference form on the basis of the patterns shown in Fig. 3 with explicit expressions for the various coefficients A’s given in the Appendix. By writing eqn (4) in finite difference form for every point (node) of the gridwork (interior or near or on the boundary) of the rectangular plate superelement of Fig. 2 with free edges on the basis of eqn (7) and the patterns of Fig. 3 one can obtain a system of linear equations of the form WI{*} = {P},
= {P},
between the vector {6*} of the boundary nodal displacements and rotations (one vertical deflection and two rotations along the x and y axes per comer node and one vertical deflection and one normal slope per side node) and the vector {r?*} of the vertical nodal displacements (interior and on the boundary) corresponding to the spring supported superelement, is constructed. In deriving the matrix [G] the rotation of any boundary node is approximated by the slope at the middle of the first boundary interval at the node. Let now {P} and {P} be the two systems of forces corresponding to the two sets of displacements {D*} and {@*}, respectively. Application of Betti’s theorem on these two systems and use of (11) yield {p} = [G]r{E}. (12) Assuming also that the two systems are energy equivalent, i.e.
(9)
where {C} and {P} are the transformed nodal deflection and load vectors, respectively and [m is the square, symmetric and singular transformed stiffness matrix of the superelement. It is now desirable to replace the stiffness equation (9) by an energy equivalent one which connects boundary nodal displacements and rotations with corresponding boundary nodal forces as shown in Fig. 4 in effect eliminating all interior nodal displacements and forces appearing in (9) as it was done in [12] for the static case. In order to accomplish that, at first springs of arbitrary stiffness are introduced at three or more boundary nodes not lying along the same line so as to reduce (9) to [R*](P)
(11)
+a),
(10)
; {P}r[X*]-i{P}
=; {P}rlfC]{E},
(13)
where p] is the flexibility matrix of the spring supported superelement corresponding to boundary displacements and employing (12) one can find p]
= [G][R*] -‘[Cl?
(14)
Inversion of p] gives the stiffness matrix [$‘*I corresponding to the boundary degrees of freedom of the spring supported superelement. Subtracting the spring stiffness from [s*] finally provides the stiffness matrix [S] of the free, unsupported superelement such that
rwq
= {Q.
I I
Fig. 2. A typical finite difference plate gridwork for an interior point.
WI
D. E.
BESKOS
and K. L. LEUNG
Fig. 3. Finite difference patterns for points at or near a free plate edge.
ft should be noticed that the ap~roxima~on of the rotation during the construction of [Gf makes some of the elements of p] iess accurate than the other elements. To improve the accuracy of these elements a more accurate expression for the rotation at the boundary nodes is employed as explained in[24]. The vector of the nodal forces (Fj needed in (15) is not determined from (12) because [G] cannot be inverted as being, in general, a non-square matrix, Instead (F) is determined by observing that in view of (IO), (It) and (13) (a*) = [G][KK*]-‘(PJ = F]{F),
(16)
from which
of a large rectangular plate or ph& system can be obtained by an appropriate superposition of the stiffness matrices ]$ of the various superelements comprising the system. Thus the transformed equation of motion for the plate system takes the form rQ
%$= fQ%
where fd), and {P& are the system vectors of the nodal displacements and forces, respectively. After application of the transformed boundary conditions of the plate system one can solve (IS) to obtain the transformed nodal displacements. The interior deflections of every plate superelement due to the nodal forces {Pf are then obtained from {@I, = [R*]-‘[G]‘[~*](6),
FORMULATION AND SOI,UTION OF THE PROBLEM
In the previous section the procedure for obtaining the stiffness matrix [q of a free rectangular piate superelement in the Laplace transformed domain was described. Thus the transformed stiffness matrix IQ,
(1%
(19)
which is obtained by combining eqns (lo), (12) and (16). On the other band, the interior deflections of every plate superelement produced by the interior lateral forces {Is) and the sandal reactions ( - P> acting together can be expressed, with the aid of (10)
Dynamic response of plate systems
n
Fig. 4. A plate bending superelement with (a) Interior nodaf displacements; (bj boundary nodal dispiacements and rotations.
and (12) as (“}* = [x*]-‘({P}
- [C]T{F}).
(20)
Finally, the total transformed interior deflections of every plate superelement are obtained by summing up (w)~ and {w& of (19) and (20), i.e. from {+} = K*l-‘{Is}
I- [Glr([~*l{~}
- {F})}. (21)
Once the interior and boundary nodal transformed deAections of a plate superelement are known the transformed moments and forces &ZX,JZP Mm, V, and Vyskwn in Fig. 1 can be easily obtained from eqns (5) and (6) in finite difference form as given, e.g. in [4]. The above procedure is mainly recommended for dynamically analyzing rather complicated plate systems, such as one or two-way continuous plates on rigid supports_ Plates on ffexible supports can also be treated by the proposed method if one takes into account the Laplace transformed stiffnesses of the beams supporting the plate system. Explicit expressions for these stiffnesses can be found in Beskos et al. [ 15,l S-20]. For simple rectangular plates one can directly utilize eqn (9) to solve for the transformed interior deffectiuns. However, in establishing eqn (9) in this case the finite difference patterns of Fig. 3
corresponding to free edges have to be supplemented with those of Fig. 5 corresponding to other boundary conditions[4,24]. Another way is, of course, to model the plate by just one superelement and use eqns (15) and (21). Both of these two approaches are i&&rated in a numerical example, The transformed load vector {P) used in the analysis consists of nodal concentrated loadelements, which for a uniform spatial variation of the distributed load p, i.e. p =p(r), have the form ,FFi= F(S)&+,, Pi = pl(~)&/2 and Pi 3: @(~)lt,il,/4 for the node i being inside, at an edge or at a corner of the gridwork, respectively. In cases of nonuniform spatial variations of the load p, i.e. when p = p (x, y, t), mm complicated expressions are needed for computing the nodal concentrated loadelements of the vector {p}. These expressions can be found in [3,4]. According to the preceding discussion, eqns (9) or (15) and (21) can be used for the determination of the transformed displacement vector {$I of one plate superelement, while for a plate system one should employ eqns (18) and (21). A numerical inversion of the transformed solution will finally provide the time domain response. This implies that eqns (9) or (21) have to be solved numerically for a sequence of values of the transformed parameter s, It should be also noticed that in cases where the time va&tion of the
768
II. E. BESKOS and K. L.
bUNG
k--1x-h
A
I
T
I
hY
c
Implysb~#orlrd edgt
E
i? D 71”““““/“““.‘/“/‘/“““““““~‘fl /
b
;c
A-
(b)
WC=-W,
Fig. 5. Boundary conditions for built-in and simply supported plate edges.
applied dynamic load is too complicated for an analytic or through tables determination of its Laplace transform, a numerical direct Laplace transform is necessary. A very detailed account on the numerical direct and inverse Lapiace transform and its applications to various time-dependent linear problems can be found in Narayanan and Beskos[l3]. A comparative study in [13] proved that Durbin’s algorithm [2S) is the most accurate and most convenient to use for both inverse and direct Laplace transform ~rnpu~tio~ in structural dynamics. The only disadvantage of this algorithm which works with complex data (complex values of s) is that this is more costly than algorithms working with real data. Real data algorithms however, become useless for dynamic problems involving oscillatory responses[ I3]. The Durbin’s ~go~~[25] employed in this work combines both finite Fourier cosine and sine transforms to establish the following Laplace inversion formula for a function f(t) whose transform is 3’(s): f(r,) = 2(evjA’/T) I
where s, = y + in(2x/T), W = ei2*‘Q, A(n) = i
ReIi(y + i(n + lQ)f2nlT))},
I=0
B(n)
=
i
IrnCT’(r + i(n + 1(2)(2nlT))},
(23)
I=0 !,a
=
O,l,Z,
. . . .
i = VT
and f(t) is computed for Q equidistant points fi = jdr = jT/Q,j = 0,1,2, . . . , Q - 1. It is suggested that for L x Q ranging from 50-5000 one should select yT = 5-10 for good results, where T is the total time interval of interest. The computations involved in (22) are performed by employing the fast Fourier transform algorithm of Cooiey and Tukeyf26]. A modification of the above algorithm can lead to an algorithm for the direct numerical Laplace transform in the form[l3] BY + %Jd = - ff(l -t 5
n=l
= 0) (A(n) + mw’“,
(24)
Dynamic response of plate systems where f{(n + lQ)At}e-"("+Q)klAt,
A(n) = i I=0
B(n) = 0, W = e-2”i’Q, k, n = 1, 2, , . . , Q
(25)
~~ = 2nklT, T = Q x L x At.
EFFECTS OF DAMPING AND ELASTIC FOUNDATION INTERACTION
The effects of both the external viscous damping and the internal v&o-elastic damping on the dynamic response of the plate are considered in this section. In the case where there is external viscous damping, eqn (1) is replaced by
a4w a34 a9af+ay”+N
$+2-
m a%
c
-+Ndr=g’ at2
aw
(26)
where c is the viscous damping coefficient. Application of the Laplace transform with respect to time under zero initial conditions on (26) yields
a%
a%
2+&i+ a2ay
(WISP+ cs) _ N
p
w=N.
(27)
A comparison of eqns (4) and (27) reveals that external viscous damping can be taken into account very easily if (ms2) in (4) is replaced by (ms2 + cs) or equivalently if p in (8) is replaced by /3 = (~,4/N)(m.Y2+ 0).
i, j = 1,2,3,
(29)
where fl is the shear modulus, f is the internal damping coefficient and the deviatoric stress and strain tensors sI/ and ei/, respectively, are defined in terms of the stress and strain tensors eij and 4, respectively, by sii = bii - 6Jaii/3), eij = cii - &(cd3),
(30)
with 6, being Kronecker’s delta and repeated indices indicating summation, the correspondence principle can be stated as follows: the Laplace transform of the viscoelastic solution can be obtained from the Laplace transform of the elastic solution by replacing the elastic constants p and 1 (Lam6 constants) by[27] cc’= ~(1 +fi), 1’ = 1 + (2/3)b -$),
,_3v-(l-2vfi 3E(l +fi) 3 + (1 - 2v&’ v - 3 + (1 - 2vfi ’
4
$+2-
a%
ax2ay2
where 5 is the foundation modulus. Application of the Laplace transform with respect to time under zero initial conditions on (33) yields
a% a43 adiG (ms2+0 jg+2- ax2ay2 + ay4+ N
_ d w y.
(34)
A comparison of eqns (4) and (34) clearly indicates that the effect of the soil foundation interaction can be very easily taken into account in the analysis by simply replacing (ms’) in (4) by (mr2 + r) or equivalently by replacing /I in (8) by
B = (A4/W(~’ + 0
(35)
Obviously, the combined effects of damping and soil foundation interaction can also be easily taken into account.
NUMERICAL EXAMPLJB
In this section two numerical examples are treated in detail to illustrate the methodology proposed in this paper for dynamically analyzing plates and plate systems by the combined FDM-Laplace transform (based on eqn (9)) and the combined FDM-EMLaplace transform (based on eqns (15) or (18), and (21)) techniques. The computations were done on the CDC Cyber 730 computer of the University of Minnesota and the computer programs can be found in [241. Example 1 Consider a rectangular simply supported all around plate of plan size a x b = 60 x 40 in., uniform thickness h = 1 in., mass per unit area m = 0.00073 lb-&/in’ and elastic constants E = 30 x lo6 psi and v = 0.25, as shown in Fig. 6(a). This plate is subjected to a uniformly distributed dynamic load p(t) given by t 5 td t 2 td
(31)
or the elastic constants E and v by E, =
can be very easily included in this analysis by just replacing E and v in (2), (5) and (6) by E’ and v’ given by (32). Other more complicated viscoelastic models can also be employed without any particular difficulty[27J. In many practical situations the structural foundation consists of a simple or complicated mat which can be modeled for analysis purposes as a plate or plate system on elastic soil of the Winkler type. The equation of motion of a plate on an elastic foundation is
(28)
The effect of internal viscoelastic damping is treated with the aid of the correspondence principle, as described, e.g. in Boley and Weiner[27]. Thus, on the assumption that the plate material obeys the constitutive law of the Kelvin viscoelastic model sV= 2&e, + fde#/dt),
169
(32)
with intensity p. = 40 psi and time duration t, = 0.05 set, as shown in Fig. 6(b). By employing modal analysis one can easilv obtain the deflection at the center of the plate as[24] ’ m
where s is the Laplace transform parameter. Therefore, internal viscoelastic damping of the Kelvin type
(36)
m
.
we = c c A,sin$sin i=1 i=I
. y,
(37)
770
D. E.
BJSKOS and
K. L. LHJNG
Fig. 6. Simply supported all around plate uniformly loaded by a time dependent load.
where Aji =
16p,(DLF)/jin2mo,‘,
if i x j odd , ifi xjeven
0
The Laplace transform of p(r) of (36) needed for the analysis was obtained analytically as (38)
P(s) = P. +-(1-e-W) [
l-coso.l+-&-sinq,-J-, DLF= I
.
I
&[sin
t St,
nd
c0,t -sin
w,(t - td)]- cos
qt,
t 2
td
nd (39)
The bending moments M, and M, at the center of the plate can also be ohtained in closed form as[24]
(41)
The deflection wc of (37) and the bending moments M: and Mf, of (41) were computed for i, j = 1,2,3 from t = 0 to t = 0.08 set and it was found that the first mode solution differs from the four and nine mode solutions by 2.4 and 2.5%, resp. It was then decided to consider the four mode solution as the “exact solution” of the problem in the comparison studies which follow. Using the FDM-Laplace transform approach the same problem was solved for various mesh sizes and values of Q and for L = 1, y T = 6 and T = 0.08 sec.
d
1
(42)
The deflection wc and bending moments M: and M; at the center of the plate for a 12 x 8 finite difference mesh and Q = 80 are plotted versus time together with the exact solution and are shown in Figs. 7 and 8. Results obtained by the general purpose finite element program SAP IV[28] are also shown in Fig. 7 for 12 x 8 and 18 x 12 element discretizations. It is apparent from Figs. 7 and 8 that the proposed method shows an excellent to very good agreement for early to late times, respectively. The SAP IV results, however, as Fig. 7 clearly indicates, are very poor even for the 18 x 12 discretization which reduces the time shifting error but does not improve the response amplitude. The computer CP times spent were recorded as 90.070 set for the proposed method and 37.875 and 67.396sec for the two SAP IV discretizations. The influence of the mesh size on the results is shown in Fig. 9 which portrays the center deflection versus time for 6 x 4, 8 x 6 and 12 x 8 finite difference meshes and indicates how mesh refinement decreases the lag time in the solution at later times. The influence of the number of sample points Q needed for the numerical Laplace transform inversion on the response was found to be insignificant (differences less than 0.05%) by comparing the results for Q = 80, 160 and 240[24]. The response of the plate was also computed by incorporating two improvements in the proposed methodology, namely, taking into account the symmetry properties of the matrix [R] in (9) when solving
771
Dynamic response of plate systems -
1
fxact wlrlloa
~----
2
fDM-11
(12
-
I
SIP
IV
(16 x I21
--
4
SIP
IV
(12x II
x 6)
Fig. 7. Center deflection versus time by various methods for the plate of example 1.
-6
Fig. 8. Center internal moments versus time by various methods for the plate of example 1.
that equation in the transformed domain for complex values of s and computing the transformed solution G for only 2/3 of the 80 points (Q = 80) and providing the remaining l/3 of the points by linear interpolation. This is possible due to the fact that both the real and imaginary parts of the transformed function ii have a rea.Ily smooth variation for large values of s. This improved method provided results of the same high accuracy as before (Fig. 7) at a cost of only 63.563 set of CP computer time. The effect of internal viscoelastic damping on the dynamic response of the plate is illustrated in Fig. 10. It is shown there how the attenuating effect of damping increases for increasing values of f and time duration. The same problem was also solved by the combined
PDM-FEM-Laplace transform method. Actually, the plate was idealized both as an assemblage of two equal superelements of 30 x 40 in. size and as just one superelement. A 4 x 4 finite difference mesh for each superelement and employment of eqns (18) and (21) resulted in a response shown in Fig. 11 where it is obvious that the time shifting is unacceptable. The CP computer time spent, with improvements in the numerical inversion taken into account, was 1015.349 sec. The reason for this poor showing is, of course, the fact that the common nodes of the two superelements do not account for the tangential slope. Obviously, an addition of one more degree of freedom, namely the tangential slope, on all side nodes of the superelements could correct the result. However, this is not at all necessary in this example, because ex-
D. E.
BBSKOSand
K. L. LJWNG
--
I x6
mesh
------
IX4
,“rrh
Fig. 9. Influence of the mesh size on the center deflection history of the plate of example I.
Fig. 10. Effect of viscoelastic damping on the center deflection history of the plate of example 1.
cellent to very gocd results can be obtained much more easily by utilizing only one superelement for the whole
plate and employing eqns (15) and (21) for this purpose. The results are also shown in Fig. 11 for 4 x 4 and 8 X 6 finite difference meshes obtained at a cost of 479.072 and 3305.916 set of CP computer time, with the numerical improvements taken into account. In view of this large increase in computer time for a small increase in accuracy, one should be content with the results of the 4 x 4 discretization. In general, one can observe that results obtained by the proposed method, even when they are not considered acceptable, still provide very good estimates of the response amplitude and their only deficiency is with regard to the time shifting.
Example 2 Consider a simple rectangular one-way continuous plate system on rigid supports consisting of two equal panels as shown in Fig. 12. The two longitudinal sides of the structure have a length of 60 in. and are simply supported. The two short sides of the structure have a length of 40 in. and one of them is fixed while the other is free. The thickness of the plate system, its loading, as well as its material properties are the same with those of the plate of the previous example. The dynamic deflection wAof the point A at the middle of the free edge was determined by the FDM-FEMLaplace transform method on the basis of eqns (18) and (21) by considering one superlernent per panel
Dynamic respoese of plate syskms
Fig. 11. Center deflection versus time by the FDM-FEM-LT approach for the plate of example 1.
f
/ / / / / / / /
z / __-__-----3
fim*lc rrpprrl _____
-----_,
i
Fig. 12. Rectangular one-way continuous plate system of example 2.
with a 4 x 4 finite difference mesh and its time history is depicted in Fig. 13. An improved solution obtained by considering again two superelements but with three degrees of freedom per node is also shown in Fig. 13. Use of superelements with three degrees of freedom per node was expected to give better results here due to the fact that nodes along the free side of the plate system have both normai and tangential slopes. The CP computer times spent for the two versions, with improvements taken into account,
were 565.958 and 843.266sec, resp. Figure 13 also shows the results obtained by SAP IV for 8 x 12 and 12 x 18 finite element discretization of the whole structure by spending 65.530 and 71.248 sets of CP computer time, respectively, i.e. approx. l/IO of the computer time spent by the proposed method. To the authors’ knowledge, there is no analytic or other numerical solution available for this problem and any judgement about accuracy of the proposed method and SAP IV should be based on the detailed com-
114
D. E. BFZ,KOSand K. L. LEUNG
Fig. 13. Deflection at point A of plate of Fig. 12 versus time by various methods.
parison studies done in example 1. On the basis of that and the fact that the SAP IV results show a significant amplitude reduction with time, which is clearly incorrect, one should consider those results in significant error. It is apparent that much more refined discretizations are necessary for SAP IV to provide acceptable results, especially in view of the fact that improvements in response amplitude appear to be very small with discretization refinement. In summary, the proposed method might be time consuming in execution but requires very little time for data preparation and checking and provides results of acceptable accuracy. CONCLUSIONS
The main conclusions of the work described in this paper can be summarized as follows: (1) A general numerical method for determining the response of rectangular plates and plate systems to lateral dynamic loads has been proposed. This method is a combination of the FDM and the Laplace transform with respect to time for simple plates or a combination of the FDM, the FEM and the Laplace transform with respect to time for plate systems, where every panel is modeled by just one superelement. (2) Use of the Laplace transform reduces the dynamic problem to a static-like one which is formulated and solved in the transformed domain. A numerical inversion of the transformed solution provides the time domain response. Thus, the method does not require prior knowledge of natural frequencies and modal shapes. In general, the method provides results of acceptable accuracy by spending much less time for data preparation and checking and more computer execution time than conventional finite element techniques. (3) The importance of the method lies in that it combines the advantages of both the FDM and FEM and eliminates their respective disadvantages. Thus, very complex plates or plate systems with complicated boundary and loading conditions can be treated with acceptable accuracy and efficiency.
(4) The effect of damping, external viscous or internal viscoelastic, on the dynamic response can be very easily taken into account due to the presence of Laplace transform in the solution procedure. Actually damping can vary from superelement to superelement and it can be of any amount. Furthermore, the dynamic response of plates or plate systems on elastic foundation can be also easily determined by the proposed method again because of the employment of Laplace transform in the solution procedure.
Acknowledgements-The authors acknowledge the support provided by the University of Minnesota Computer Center. Thanks are also due to Dr. M.J. Frisch for his valuable suggestions in connection with the numerical computations.
RJWERENCES
1. S. P. Timoshenko, D. H. Young and W. Weaver, Jr., Vibration Problems in Engineering, 4th Edn. Wiley, New York (1974). 2. A. S. Veletsos and N. M. Newmark, Determination of natural frequencies of continuous plates hinged along two opposite edges. J. Appl. Mech. 23, 97-102 (1956). 3. R. Szilard, Theory and Analysis of Plates. Prentice Hall, Englewood Cliffs, New Jersey (1974). 4. A. Ghali and A. M. Neville, Structural Analysis: A Un$ed Classical and Matrix Approach. Intext Educ. Publ., Scranton, Penn. (1972). 5. R. H. Gallagher, Finite Element Analysis: Fur&mentals. Prentice Hall, Englewood Cliffs, New Jersey (1975). 6. J. S. Przemieniecki, Theory of Matrix Structural Analysis. McGraw-Hill, New York (1968). 7. C. L. Wu and Y. K. Cheung, Frequency analysis of rectangular plates continuous in one or two directions. Earth. Engng. Struct. Dyn. 3, 3-14 (1974). 8. S. French, A. P. Kabailaand V. A. Pulmano, Single element panel for flat plate structures. Proc. ASCE 101 (ST9), 18011812 (1975). 9. R. M. Gutkowski and C. K. Wang, Continuous plate analysis by finite panel method. Proc. ASCE 102 (ST3), 629643 (1976).
775
Dynamic response of plate systems 10. R. M. Gutkowski, Finite panel analysis of orthotropic plate systems. Proc. ASCE 103 (ST1 I), 2211-2224 (1977). 11. A. Ghali and K. J. Bathe, Analysis of plates subjected to in-plane forces using large finite elements. I.A.B.S.E. Publ. SO-I, 2940 (1970). 12. A. Ghali and K. J. Bathe, Analysis of plates in bending using large finite elements. Z.A.B.S.E. Publ. 30-11, 61-72 (1970). 13. G. V. Narayanan and D. E. Beskos, Numerical operational methods for time-dependent linear problems. Znt. .I. Num. Meth. Engng. 18, 1829-1854 (1982). 14. D. E. Beskos and B. A. Boley, Use of dynamic influence coefficients in forced vibration problems with the aid of Laplace transform. Comput. Structures 5, 263-269 (1975). 15. G. D. Manolis and D. E. Beskos, Thermally induced vibrations of beam structures. Comput. Meth. Appl. Mech. Engng 21, 337-355 (1980). 16. D. E. Beskos and J. B. Oates, Dynamic analysis of ring-stiffened circular cylindrical shells. J. Sound Vib. 75, 1-15 (1981). 11. G. D. Manolis and D. E. Beskos, Dynamic stress concentration studies by boundary integrals and Laplace transform. Int. J. Num. Meth. Engng 17, 573-599 (1981). 18. G. V. Narayanan and D. E. Beskos, Dynamic soilstructure interaction by numerical Laplace transform. Engng. Structures 4, 53-62 (1982). 19. G. D. Manolis and D. E. Beskos, Dynamic response of framed underground structures. Cotnput. Structures 15, 521-531 (1982). 20. D. E. Beskos and G. V. Narayanan, Dynamic response of frameworks by numerical Laplace transform. Comput. Meth. Appl. Mech. Engng 37, 289-307 (1983). 21. G. D. Manolis and D. E. Beskos, Dynamic response of lined tunnels by an isoparametric boundary element method. Comput. Meth. Appl. Mech. Engng 36, 291-307 (1983). 22. G. D. Manolis and D. E. Beskos, Dynamic stress field around a cavity embedded in a halfplane by the boundary element method. Earth. Engng Struct. Dyn., submitted. 23. D. E. Beskos and A. Y. Michael, Solution of plane transient elastodynamic problems by finite elements and Laplace transform. Comput. Structures 4, 695-701 (1984). 24. K. L. Leung, Dynamic analysis of plates using FDM-FEM and Laolace transform. M.S. Thesis. University of Minnesoti, Minneapolis, Minnesota (1982). 25. F. Durbin, Numerical inversion of Laplace transforms: an efficient improvement to Dubner’s and Abate’s method. Comput. J. 17, 371-376 (1974). 26. J. W. Cooley and J. W. Tukey, An algorithm for the machine calculation of complex Fourier series. Math. Comput. 19, 297-301 (1965): 27. B. A. Bolev and J. H. Weiner. Theorv of Thermal Stresses. Wiley, New York (196@. ’ 28. K. J. Bathe, E. L. Wilson and F. E. Peterson, SAP IV, a structural analysis program for static and dynamic response of linear systems. Rep. EERC 73-l 1, Univ. of California, Berkeley (June 1973).
APPENDIX Explicit expressions of the coefficients A’s appearing in the various finite difference patterns of Fig. 3 are given as follows[4,24]:
Al =6+6&+8a
+p
A5=a2
A2= -4(l+a)
A6=1
A3= -4a(l+a)
A7=5+6a2+8a+B
A4=2a
A9= -2(2a-va+l)
A8=a(2-v)
A10=1+4a(l-v)+3a2(l-v2)+B All=-2a(l-v+a(l-v’)) A12=ia2(l-v’)
(AlI
A13=5+5a2+8a+p Al4=2a(l Al5=
-v)
-2a(2-v+a)
Al6=1+4a(l-v)+ia’(l-v’)+/I
> Al8=2a(l-v)+i(l+a?(l-v2)+b
Al9=
-
a(, -v)+;(l
-v2) >
A20=;(1
-v’)
A21=6+5a2+8a+fl A22=a2+4a(l
-v)+3(1
-v2)+j
A23 = - 2(a(l - v) + (1 - v’)) A24=a2+4a(l-v)+i(l
-v*)+,9
Conrpruers& Suucmres Vol. 19, No. 5/6, pp. 163-115, 1984 Printedin the U.S.A.
DYNAMIC RESPONSE OF PLATE SYSTEMS BY COMBINING FINITE DIFFERENCES, FINITE ELEMENTS AND LAPLACE TRANSFORM
D. E. BE~KOS~ Department of Civil Engineering, University of Patras, Patras, Greece and K. L. LEUNG$ Department of Civil and Mineral Engineering, University of Minnesota, Minneapolis, MN 55455, U.S.A. (Received 15 September 1983; received for publication 5 December
1983)
Abstract-A numerical method for the determination of the dynamic response of large rectangular plates or plate systems to lateral loads is proposed. The method is a combination of the finite difference method, the finite element method and the Laplace transform with respect to time. The plate system is considered as an assemblage of a small number of big rectangular superelements whose stiffness matrices are derived with the aid of the finite difference method in the Laplace transform domain. These superelements are then used to formulate and solve the problem by the finite element method in the transformed domain. The dynamic response is finally obtained by a numerical inversion of the transformed solution. External viscous or internal viscoelastic damping as well as the elastic foundation interaction effect can easily be taken into account. The method is illustrated and its merits demonstrated by means of numerical examples.
INTRODUCIION
The determination of the dynamic response of rectangular plate systems encountered in ship, aerospace and civil engineering structures is an important technological problem. The complex geometry and boundary and loading conditions of single large rectangular plates or plate systems such as one or two-way continuous plates on rigid or flexible supports do not permit an analytic solution of the problem but in very few simple cases[l-31. Static and dynamic analysis of complex plate systems is usually accomplished by numerical methods such as the Finite Difference Method (FDM) and the Finite Element Method (FEM)[3-51. However, both of these methods require a large computer memory for problems dealing with large plate systems due to the fact that a refined discretization is required for acceptable results. Besides, round-off error accumulation during the solution of large systems of equations usually leads to unsatisfactory results. In an effort to overcome these difficulties various special techniques have been employed in conjunction with the FEM for the analysis of plate systems such as the substructuring technique[6], the finite strip method[4, 71 and the finite panel method[8-lo]. With the exception of reference[l, Refs. [&lo] dealt with the static problem only. The finite strip and finite panel methods require a much smaller total number TProfessor and Visiting Professor at the Department of Civil and Mineral Engineering, University of Minnesota, Minneapolis, MN 55455, U.S.A. SGraduate Student.
of degrees of freedom for the whole system but they employ rather complicated displacement functions. It is apparent that a finite difference gridwork for analyzing a plate corresponds to a number of degrees of freedom which is much smaller than that for a finite element mesh with as many nodes as the gridwork points in the plate. On the other hand it is well known that treatment of the boundary conditions by the FEM is much simpler and easier than by the FDM. These two facts suggest that a combination of the two methods in such a way as to exploit their advantages and suppress their disadvantages might offer a solution to the problem of analyzing large plate systems. Indeed this is exactly what Ghali and Bathe[l 1, 121 did for the static analysis of plate systems under in-plane and lateral loading. In this work, the determination of the flexural dynamic response of large rectangular plates or plate systems to lateral loads is accomplished by a combination of the FDM, the FEM and the Laplace transform with respect to time. It is actually an extension to the dynamic case of the work of Ghali and Bathe[l2] with the aid of Laplace transform, which effectively reduces the dynamic problem to a static one. According to the proposed method the plate system is idealized as an assemblage of a small number of big rectangular superelements whose stiffness matrices are derived by the FDM on the basis of the Laplace transformed equation of motion. Use of these superelements in conjunction with the FEM leads to the formulation and solution of the problem in terms of nodal displacements in the transformed domain. Interior displacements and
763
764
D. E. BE~KOSand K. L. LEUNG
stresses for a superelement can then be obtained through relations connecting the boundary nodal displacements with the interior ones. The structural response can finally be obtained by a numerical inversion of the transformed solution. Use of Laplace transform has the advantages that the dynamic problem is effectively reduced to a static one, no need of natural frequencies and modal shapes is required as in modal analysis and other effects such as external viscous or internal viscoelastic damping as well as the foundation interaction can be easily taken into account. Dynamic forces of any time variation can be handled without any difficulty by employing a direct numerical Laplace transform. A comprehensive study on the direct and inverse numerical Laplace transform as applied to complex time-dependent linear problems of engineering science can be found in Narayanan and Beskos [ 131, while various applications of the Laplace transform to structural and soil dynamics problems can be found in Beskos et al.[14-231. The proposed method is illustrated by means of numerical examples and its merits are carefully assessed.
modulus of elasticity and Poi_sson’s ratio of the plate. The Laplace transform f(x, y, S) of a function f(x, y, t) with respect to time is defined by
fix,Y, s) =
sm
e-s’f(w,~) dr,
(3)
0
where s, in general a complex number, is the Laplace transform parameter. Application of the Laplace transform with respect to time on eqn (1) under zero initial conditions yields 4
g+2
a% a% _+_+,w=P. ax2ay2 ay4
-
d
N
(4)
It is very easy to see[24] that the force-displacement relations for the plate element of Fig. 1 in the Laplace transformed domain become
STIFFNESS MATRIX OF RECTANGULARPLATE SUPERELEMENT
Consider a rectangular plate bending element under distributed lateral load p(x, y, t) as shown in Fig. 1. The equation of its lateral flexural motion according to the linear theory of Kirchhoff-Love takes the form[l, 31 2+2-
adw saw M azw ax2ay2 +v+;dr”=$
p (1)
where w = w(x, y, t) is the lateral deflection, m is the plate mass per unit surface, t indicates time and the bending rigidity N is given by
where the moments k,, “y! &, and shear forces v, and v, are depicted in the positive sense in Fig. 1. Consider now a large rectangular plate bending element in the Laplace transformed domain with a typical finite difference gridwork as shown in Fig. 2. On the basis of that gridwork the finite difference form of eqn (4) takes ihe form[24] A,@,
+ A2(Gi+
I,j + wi_ Ij) + Aj(Gi,j+
1 + W,j_ 1)
+A4(~i+,,j+,+Wi+I,j-I+~i-I,j+I+~i-I,j-1)
N = Eh’/12(1 - v2),
(2) +
A.5(w&j+2
+
fi,j_2)
+
Ae(Wi+*,j
+ Wf_2J) = 0,
with h being the plate thickness and E and v being the
Fig. 1. A typical plate bending element with its corresponding forces.
(7)
lateral load and internal
Dynamic response of plate systems
765
and the geometric relationship
where
{a*} = [G]{G*}, A, =6+6a2+8a A,= -4a(l
+/?, A,= -4(l
+a), A,=2a,
As=a2,
a = (~,/~J*, j3 = (lQN)mS*.
A,= 1, (8)
Equation (7) is valid for an interior point (i,j) such as that of Fig. 2. For other points (i,j) close to or on the boundaries of the plate element, which are assumed to be free, eqn (4) can be written in finite difference form on the basis of the patterns shown in Fig. 3 with explicit expressions for the various coefficients A’s given in the Appendix. By writing eqn (4) in finite difference form for every point (node) of the gridwork (interior or near or on the boundary) of the rectangular plate superelement of Fig. 2 with free edges on the basis of eqn (7) and the patterns of Fig. 3 one can obtain a system of linear equations of the form WI{*} = {P},
= {P},
between the vector {6*} of the boundary nodal displacements and rotations (one vertical deflection and two rotations along the x and y axes per comer node and one vertical deflection and one normal slope per side node) and the vector {r?*} of the vertical nodal displacements (interior and on the boundary) corresponding to the spring supported superelement, is constructed. In deriving the matrix [G] the rotation of any boundary node is approximated by the slope at the middle of the first boundary interval at the node. Let now {P} and {P} be the two systems of forces corresponding to the two sets of displacements {D*} and {@*}, respectively. Application of Betti’s theorem on these two systems and use of (11) yield {p} = [G]r{E}. (12) Assuming also that the two systems are energy equivalent, i.e.
(9)
where {C} and {P} are the transformed nodal deflection and load vectors, respectively and [m is the square, symmetric and singular transformed stiffness matrix of the superelement. It is now desirable to replace the stiffness equation (9) by an energy equivalent one which connects boundary nodal displacements and rotations with corresponding boundary nodal forces as shown in Fig. 4 in effect eliminating all interior nodal displacements and forces appearing in (9) as it was done in [12] for the static case. In order to accomplish that, at first springs of arbitrary stiffness are introduced at three or more boundary nodes not lying along the same line so as to reduce (9) to [R*](P)
(11)
+a),
(10)
; {P}r[X*]-i{P}
=; {P}rlfC]{E},
(13)
where p] is the flexibility matrix of the spring supported superelement corresponding to boundary displacements and employing (12) one can find p]
= [G][R*] -‘[Cl?
(14)
Inversion of p] gives the stiffness matrix [$‘*I corresponding to the boundary degrees of freedom of the spring supported superelement. Subtracting the spring stiffness from [s*] finally provides the stiffness matrix [S] of the free, unsupported superelement such that
rwq
= {Q.
I I
Fig. 2. A typical finite difference plate gridwork for an interior point.
WI
D. E.
BESKOS
and K. L. LEUNG
Fig. 3. Finite difference patterns for points at or near a free plate edge.
ft should be noticed that the ap~roxima~on of the rotation during the construction of [Gf makes some of the elements of p] iess accurate than the other elements. To improve the accuracy of these elements a more accurate expression for the rotation at the boundary nodes is employed as explained in[24]. The vector of the nodal forces (Fj needed in (15) is not determined from (12) because [G] cannot be inverted as being, in general, a non-square matrix, Instead (F) is determined by observing that in view of (IO), (It) and (13) (a*) = [G][KK*]-‘(PJ = F]{F),
(16)
from which
of a large rectangular plate or ph& system can be obtained by an appropriate superposition of the stiffness matrices ]$ of the various superelements comprising the system. Thus the transformed equation of motion for the plate system takes the form rQ
%$= fQ%
where fd), and {P& are the system vectors of the nodal displacements and forces, respectively. After application of the transformed boundary conditions of the plate system one can solve (IS) to obtain the transformed nodal displacements. The interior deflections of every plate superelement due to the nodal forces {Pf are then obtained from {@I, = [R*]-‘[G]‘[~*](6),
FORMULATION AND SOI,UTION OF THE PROBLEM
In the previous section the procedure for obtaining the stiffness matrix [q of a free rectangular piate superelement in the Laplace transformed domain was described. Thus the transformed stiffness matrix IQ,
(1%
(19)
which is obtained by combining eqns (lo), (12) and (16). On the other band, the interior deflections of every plate superelement produced by the interior lateral forces {Is) and the sandal reactions ( - P> acting together can be expressed, with the aid of (10)
Dynamic response of plate systems
n
Fig. 4. A plate bending superelement with (a) Interior nodaf displacements; (bj boundary nodal dispiacements and rotations.
and (12) as (“}* = [x*]-‘({P}
- [C]T{F}).
(20)
Finally, the total transformed interior deflections of every plate superelement are obtained by summing up (w)~ and {w& of (19) and (20), i.e. from {+} = K*l-‘{Is}
I- [Glr([~*l{~}
- {F})}. (21)
Once the interior and boundary nodal transformed deAections of a plate superelement are known the transformed moments and forces &ZX,JZP Mm, V, and Vyskwn in Fig. 1 can be easily obtained from eqns (5) and (6) in finite difference form as given, e.g. in [4]. The above procedure is mainly recommended for dynamically analyzing rather complicated plate systems, such as one or two-way continuous plates on rigid supports_ Plates on ffexible supports can also be treated by the proposed method if one takes into account the Laplace transformed stiffnesses of the beams supporting the plate system. Explicit expressions for these stiffnesses can be found in Beskos et al. [ 15,l S-20]. For simple rectangular plates one can directly utilize eqn (9) to solve for the transformed interior deffectiuns. However, in establishing eqn (9) in this case the finite difference patterns of Fig. 3
corresponding to free edges have to be supplemented with those of Fig. 5 corresponding to other boundary conditions[4,24]. Another way is, of course, to model the plate by just one superelement and use eqns (15) and (21). Both of these two approaches are i&&rated in a numerical example, The transformed load vector {P) used in the analysis consists of nodal concentrated loadelements, which for a uniform spatial variation of the distributed load p, i.e. p =p(r), have the form ,FFi= F(S)&+,, Pi = pl(~)&/2 and Pi 3: @(~)lt,il,/4 for the node i being inside, at an edge or at a corner of the gridwork, respectively. In cases of nonuniform spatial variations of the load p, i.e. when p = p (x, y, t), mm complicated expressions are needed for computing the nodal concentrated loadelements of the vector {p}. These expressions can be found in [3,4]. According to the preceding discussion, eqns (9) or (15) and (21) can be used for the determination of the transformed displacement vector {$I of one plate superelement, while for a plate system one should employ eqns (18) and (21). A numerical inversion of the transformed solution will finally provide the time domain response. This implies that eqns (9) or (21) have to be solved numerically for a sequence of values of the transformed parameter s, It should be also noticed that in cases where the time va&tion of the
768
II. E. BESKOS and K. L.
bUNG
k--1x-h
A
I
T
I
hY
c
Implysb~#orlrd edgt
E
i? D 71”““““/“““.‘/“/‘/“““““““~‘fl /
b
;c
A-
(b)
WC=-W,
Fig. 5. Boundary conditions for built-in and simply supported plate edges.
applied dynamic load is too complicated for an analytic or through tables determination of its Laplace transform, a numerical direct Laplace transform is necessary. A very detailed account on the numerical direct and inverse Lapiace transform and its applications to various time-dependent linear problems can be found in Narayanan and Beskos[l3]. A comparative study in [13] proved that Durbin’s algorithm [2S) is the most accurate and most convenient to use for both inverse and direct Laplace transform ~rnpu~tio~ in structural dynamics. The only disadvantage of this algorithm which works with complex data (complex values of s) is that this is more costly than algorithms working with real data. Real data algorithms however, become useless for dynamic problems involving oscillatory responses[ I3]. The Durbin’s ~go~~[25] employed in this work combines both finite Fourier cosine and sine transforms to establish the following Laplace inversion formula for a function f(t) whose transform is 3’(s): f(r,) = 2(evjA’/T) I
where s, = y + in(2x/T), W = ei2*‘Q, A(n) = i
ReIi(y + i(n + lQ)f2nlT))},
I=0
B(n)
=
i
IrnCT’(r + i(n + 1(2)(2nlT))},
(23)
I=0 !,a
=
O,l,Z,
. . . .
i = VT
and f(t) is computed for Q equidistant points fi = jdr = jT/Q,j = 0,1,2, . . . , Q - 1. It is suggested that for L x Q ranging from 50-5000 one should select yT = 5-10 for good results, where T is the total time interval of interest. The computations involved in (22) are performed by employing the fast Fourier transform algorithm of Cooiey and Tukeyf26]. A modification of the above algorithm can lead to an algorithm for the direct numerical Laplace transform in the form[l3] BY + %Jd = - ff(l -t 5
n=l
= 0) (A(n) + mw’“,
(24)
Dynamic response of plate systems where f{(n + lQ)At}e-"("+Q)klAt,
A(n) = i I=0
B(n) = 0, W = e-2”i’Q, k, n = 1, 2, , . . , Q
(25)
~~ = 2nklT, T = Q x L x At.
EFFECTS OF DAMPING AND ELASTIC FOUNDATION INTERACTION
The effects of both the external viscous damping and the internal v&o-elastic damping on the dynamic response of the plate are considered in this section. In the case where there is external viscous damping, eqn (1) is replaced by
a4w a34 a9af+ay”+N
$+2-
m a%
c
-+Ndr=g’ at2
aw
(26)
where c is the viscous damping coefficient. Application of the Laplace transform with respect to time under zero initial conditions on (26) yields
a%
a%
2+&i+ a2ay
(WISP+ cs) _ N
p
w=N.
(27)
A comparison of eqns (4) and (27) reveals that external viscous damping can be taken into account very easily if (ms2) in (4) is replaced by (ms2 + cs) or equivalently if p in (8) is replaced by /3 = (~,4/N)(m.Y2+ 0).
i, j = 1,2,3,
(29)
where fl is the shear modulus, f is the internal damping coefficient and the deviatoric stress and strain tensors sI/ and ei/, respectively, are defined in terms of the stress and strain tensors eij and 4, respectively, by sii = bii - 6Jaii/3), eij = cii - &(cd3),
(30)
with 6, being Kronecker’s delta and repeated indices indicating summation, the correspondence principle can be stated as follows: the Laplace transform of the viscoelastic solution can be obtained from the Laplace transform of the elastic solution by replacing the elastic constants p and 1 (Lam6 constants) by[27] cc’= ~(1 +fi), 1’ = 1 + (2/3)b -$),
,_3v-(l-2vfi 3E(l +fi) 3 + (1 - 2v&’ v - 3 + (1 - 2vfi ’
4
$+2-
a%
ax2ay2
where 5 is the foundation modulus. Application of the Laplace transform with respect to time under zero initial conditions on (33) yields
a% a43 adiG (ms2+0 jg+2- ax2ay2 + ay4+ N
_ d w y.
(34)
A comparison of eqns (4) and (34) clearly indicates that the effect of the soil foundation interaction can be very easily taken into account in the analysis by simply replacing (ms’) in (4) by (mr2 + r) or equivalently by replacing /I in (8) by
B = (A4/W(~’ + 0
(35)
Obviously, the combined effects of damping and soil foundation interaction can also be easily taken into account.
NUMERICAL EXAMPLJB
In this section two numerical examples are treated in detail to illustrate the methodology proposed in this paper for dynamically analyzing plates and plate systems by the combined FDM-Laplace transform (based on eqn (9)) and the combined FDM-EMLaplace transform (based on eqns (15) or (18), and (21)) techniques. The computations were done on the CDC Cyber 730 computer of the University of Minnesota and the computer programs can be found in [241. Example 1 Consider a rectangular simply supported all around plate of plan size a x b = 60 x 40 in., uniform thickness h = 1 in., mass per unit area m = 0.00073 lb-&/in’ and elastic constants E = 30 x lo6 psi and v = 0.25, as shown in Fig. 6(a). This plate is subjected to a uniformly distributed dynamic load p(t) given by t 5 td t 2 td
(31)
or the elastic constants E and v by E, =
can be very easily included in this analysis by just replacing E and v in (2), (5) and (6) by E’ and v’ given by (32). Other more complicated viscoelastic models can also be employed without any particular difficulty[27J. In many practical situations the structural foundation consists of a simple or complicated mat which can be modeled for analysis purposes as a plate or plate system on elastic soil of the Winkler type. The equation of motion of a plate on an elastic foundation is
(28)
The effect of internal viscoelastic damping is treated with the aid of the correspondence principle, as described, e.g. in Boley and Weiner[27]. Thus, on the assumption that the plate material obeys the constitutive law of the Kelvin viscoelastic model sV= 2&e, + fde#/dt),
169
(32)
with intensity p. = 40 psi and time duration t, = 0.05 set, as shown in Fig. 6(b). By employing modal analysis one can easilv obtain the deflection at the center of the plate as[24] ’ m
where s is the Laplace transform parameter. Therefore, internal viscoelastic damping of the Kelvin type
(36)
m
.
we = c c A,sin$sin i=1 i=I
. y,
(37)
770
D. E.
BJSKOS and
K. L. LHJNG
Fig. 6. Simply supported all around plate uniformly loaded by a time dependent load.
where Aji =
16p,(DLF)/jin2mo,‘,
if i x j odd , ifi xjeven
0
The Laplace transform of p(r) of (36) needed for the analysis was obtained analytically as (38)
P(s) = P. +-(1-e-W) [
l-coso.l+-&-sinq,-J-, DLF= I
.
I
&[sin
t St,
nd
c0,t -sin
w,(t - td)]- cos
qt,
t 2
td
nd (39)
The bending moments M, and M, at the center of the plate can also be ohtained in closed form as[24]
(41)
The deflection wc of (37) and the bending moments M: and Mf, of (41) were computed for i, j = 1,2,3 from t = 0 to t = 0.08 set and it was found that the first mode solution differs from the four and nine mode solutions by 2.4 and 2.5%, resp. It was then decided to consider the four mode solution as the “exact solution” of the problem in the comparison studies which follow. Using the FDM-Laplace transform approach the same problem was solved for various mesh sizes and values of Q and for L = 1, y T = 6 and T = 0.08 sec.
d
1
(42)
The deflection wc and bending moments M: and M; at the center of the plate for a 12 x 8 finite difference mesh and Q = 80 are plotted versus time together with the exact solution and are shown in Figs. 7 and 8. Results obtained by the general purpose finite element program SAP IV[28] are also shown in Fig. 7 for 12 x 8 and 18 x 12 element discretizations. It is apparent from Figs. 7 and 8 that the proposed method shows an excellent to very good agreement for early to late times, respectively. The SAP IV results, however, as Fig. 7 clearly indicates, are very poor even for the 18 x 12 discretization which reduces the time shifting error but does not improve the response amplitude. The computer CP times spent were recorded as 90.070 set for the proposed method and 37.875 and 67.396sec for the two SAP IV discretizations. The influence of the mesh size on the results is shown in Fig. 9 which portrays the center deflection versus time for 6 x 4, 8 x 6 and 12 x 8 finite difference meshes and indicates how mesh refinement decreases the lag time in the solution at later times. The influence of the number of sample points Q needed for the numerical Laplace transform inversion on the response was found to be insignificant (differences less than 0.05%) by comparing the results for Q = 80, 160 and 240[24]. The response of the plate was also computed by incorporating two improvements in the proposed methodology, namely, taking into account the symmetry properties of the matrix [R] in (9) when solving
771
Dynamic response of plate systems -
1
fxact wlrlloa
~----
2
fDM-11
(12
-
I
SIP
IV
(16 x I21
--
4
SIP
IV
(12x II
x 6)
Fig. 7. Center deflection versus time by various methods for the plate of example 1.
-6
Fig. 8. Center internal moments versus time by various methods for the plate of example 1.
that equation in the transformed domain for complex values of s and computing the transformed solution G for only 2/3 of the 80 points (Q = 80) and providing the remaining l/3 of the points by linear interpolation. This is possible due to the fact that both the real and imaginary parts of the transformed function ii have a rea.Ily smooth variation for large values of s. This improved method provided results of the same high accuracy as before (Fig. 7) at a cost of only 63.563 set of CP computer time. The effect of internal viscoelastic damping on the dynamic response of the plate is illustrated in Fig. 10. It is shown there how the attenuating effect of damping increases for increasing values of f and time duration. The same problem was also solved by the combined
PDM-FEM-Laplace transform method. Actually, the plate was idealized both as an assemblage of two equal superelements of 30 x 40 in. size and as just one superelement. A 4 x 4 finite difference mesh for each superelement and employment of eqns (18) and (21) resulted in a response shown in Fig. 11 where it is obvious that the time shifting is unacceptable. The CP computer time spent, with improvements in the numerical inversion taken into account, was 1015.349 sec. The reason for this poor showing is, of course, the fact that the common nodes of the two superelements do not account for the tangential slope. Obviously, an addition of one more degree of freedom, namely the tangential slope, on all side nodes of the superelements could correct the result. However, this is not at all necessary in this example, because ex-
D. E.
BBSKOSand
K. L. LJWNG
--
I x6
mesh
------
IX4
,“rrh
Fig. 9. Influence of the mesh size on the center deflection history of the plate of example I.
Fig. 10. Effect of viscoelastic damping on the center deflection history of the plate of example 1.
cellent to very gocd results can be obtained much more easily by utilizing only one superelement for the whole
plate and employing eqns (15) and (21) for this purpose. The results are also shown in Fig. 11 for 4 x 4 and 8 X 6 finite difference meshes obtained at a cost of 479.072 and 3305.916 set of CP computer time, with the numerical improvements taken into account. In view of this large increase in computer time for a small increase in accuracy, one should be content with the results of the 4 x 4 discretization. In general, one can observe that results obtained by the proposed method, even when they are not considered acceptable, still provide very good estimates of the response amplitude and their only deficiency is with regard to the time shifting.
Example 2 Consider a simple rectangular one-way continuous plate system on rigid supports consisting of two equal panels as shown in Fig. 12. The two longitudinal sides of the structure have a length of 60 in. and are simply supported. The two short sides of the structure have a length of 40 in. and one of them is fixed while the other is free. The thickness of the plate system, its loading, as well as its material properties are the same with those of the plate of the previous example. The dynamic deflection wAof the point A at the middle of the free edge was determined by the FDM-FEMLaplace transform method on the basis of eqns (18) and (21) by considering one superlernent per panel
Dynamic respoese of plate syskms
Fig. 11. Center deflection versus time by the FDM-FEM-LT approach for the plate of example 1.
f
/ / / / / / / /
z / __-__-----3
fim*lc rrpprrl _____
-----_,
i
Fig. 12. Rectangular one-way continuous plate system of example 2.
with a 4 x 4 finite difference mesh and its time history is depicted in Fig. 13. An improved solution obtained by considering again two superelements but with three degrees of freedom per node is also shown in Fig. 13. Use of superelements with three degrees of freedom per node was expected to give better results here due to the fact that nodes along the free side of the plate system have both normai and tangential slopes. The CP computer times spent for the two versions, with improvements taken into account,
were 565.958 and 843.266sec, resp. Figure 13 also shows the results obtained by SAP IV for 8 x 12 and 12 x 18 finite element discretization of the whole structure by spending 65.530 and 71.248 sets of CP computer time, respectively, i.e. approx. l/IO of the computer time spent by the proposed method. To the authors’ knowledge, there is no analytic or other numerical solution available for this problem and any judgement about accuracy of the proposed method and SAP IV should be based on the detailed com-
114
D. E. BFZ,KOSand K. L. LEUNG
Fig. 13. Deflection at point A of plate of Fig. 12 versus time by various methods.
parison studies done in example 1. On the basis of that and the fact that the SAP IV results show a significant amplitude reduction with time, which is clearly incorrect, one should consider those results in significant error. It is apparent that much more refined discretizations are necessary for SAP IV to provide acceptable results, especially in view of the fact that improvements in response amplitude appear to be very small with discretization refinement. In summary, the proposed method might be time consuming in execution but requires very little time for data preparation and checking and provides results of acceptable accuracy. CONCLUSIONS
The main conclusions of the work described in this paper can be summarized as follows: (1) A general numerical method for determining the response of rectangular plates and plate systems to lateral dynamic loads has been proposed. This method is a combination of the FDM and the Laplace transform with respect to time for simple plates or a combination of the FDM, the FEM and the Laplace transform with respect to time for plate systems, where every panel is modeled by just one superelement. (2) Use of the Laplace transform reduces the dynamic problem to a static-like one which is formulated and solved in the transformed domain. A numerical inversion of the transformed solution provides the time domain response. Thus, the method does not require prior knowledge of natural frequencies and modal shapes. In general, the method provides results of acceptable accuracy by spending much less time for data preparation and checking and more computer execution time than conventional finite element techniques. (3) The importance of the method lies in that it combines the advantages of both the FDM and FEM and eliminates their respective disadvantages. Thus, very complex plates or plate systems with complicated boundary and loading conditions can be treated with acceptable accuracy and efficiency.
(4) The effect of damping, external viscous or internal viscoelastic, on the dynamic response can be very easily taken into account due to the presence of Laplace transform in the solution procedure. Actually damping can vary from superelement to superelement and it can be of any amount. Furthermore, the dynamic response of plates or plate systems on elastic foundation can be also easily determined by the proposed method again because of the employment of Laplace transform in the solution procedure.
Acknowledgements-The authors acknowledge the support provided by the University of Minnesota Computer Center. Thanks are also due to Dr. M.J. Frisch for his valuable suggestions in connection with the numerical computations.
RJWERENCES
1. S. P. Timoshenko, D. H. Young and W. Weaver, Jr., Vibration Problems in Engineering, 4th Edn. Wiley, New York (1974). 2. A. S. Veletsos and N. M. Newmark, Determination of natural frequencies of continuous plates hinged along two opposite edges. J. Appl. Mech. 23, 97-102 (1956). 3. R. Szilard, Theory and Analysis of Plates. Prentice Hall, Englewood Cliffs, New Jersey (1974). 4. A. Ghali and A. M. Neville, Structural Analysis: A Un$ed Classical and Matrix Approach. Intext Educ. Publ., Scranton, Penn. (1972). 5. R. H. Gallagher, Finite Element Analysis: Fur&mentals. Prentice Hall, Englewood Cliffs, New Jersey (1975). 6. J. S. Przemieniecki, Theory of Matrix Structural Analysis. McGraw-Hill, New York (1968). 7. C. L. Wu and Y. K. Cheung, Frequency analysis of rectangular plates continuous in one or two directions. Earth. Engng. Struct. Dyn. 3, 3-14 (1974). 8. S. French, A. P. Kabailaand V. A. Pulmano, Single element panel for flat plate structures. Proc. ASCE 101 (ST9), 18011812 (1975). 9. R. M. Gutkowski and C. K. Wang, Continuous plate analysis by finite panel method. Proc. ASCE 102 (ST3), 629643 (1976).
775
Dynamic response of plate systems 10. R. M. Gutkowski, Finite panel analysis of orthotropic plate systems. Proc. ASCE 103 (ST1 I), 2211-2224 (1977). 11. A. Ghali and K. J. Bathe, Analysis of plates subjected to in-plane forces using large finite elements. I.A.B.S.E. Publ. SO-I, 2940 (1970). 12. A. Ghali and K. J. Bathe, Analysis of plates in bending using large finite elements. Z.A.B.S.E. Publ. 30-11, 61-72 (1970). 13. G. V. Narayanan and D. E. Beskos, Numerical operational methods for time-dependent linear problems. Znt. .I. Num. Meth. Engng. 18, 1829-1854 (1982). 14. D. E. Beskos and B. A. Boley, Use of dynamic influence coefficients in forced vibration problems with the aid of Laplace transform. Comput. Structures 5, 263-269 (1975). 15. G. D. Manolis and D. E. Beskos, Thermally induced vibrations of beam structures. Comput. Meth. Appl. Mech. Engng 21, 337-355 (1980). 16. D. E. Beskos and J. B. Oates, Dynamic analysis of ring-stiffened circular cylindrical shells. J. Sound Vib. 75, 1-15 (1981). 11. G. D. Manolis and D. E. Beskos, Dynamic stress concentration studies by boundary integrals and Laplace transform. Int. J. Num. Meth. Engng 17, 573-599 (1981). 18. G. V. Narayanan and D. E. Beskos, Dynamic soilstructure interaction by numerical Laplace transform. Engng. Structures 4, 53-62 (1982). 19. G. D. Manolis and D. E. Beskos, Dynamic response of framed underground structures. Cotnput. Structures 15, 521-531 (1982). 20. D. E. Beskos and G. V. Narayanan, Dynamic response of frameworks by numerical Laplace transform. Comput. Meth. Appl. Mech. Engng 37, 289-307 (1983). 21. G. D. Manolis and D. E. Beskos, Dynamic response of lined tunnels by an isoparametric boundary element method. Comput. Meth. Appl. Mech. Engng 36, 291-307 (1983). 22. G. D. Manolis and D. E. Beskos, Dynamic stress field around a cavity embedded in a halfplane by the boundary element method. Earth. Engng Struct. Dyn., submitted. 23. D. E. Beskos and A. Y. Michael, Solution of plane transient elastodynamic problems by finite elements and Laplace transform. Comput. Structures 4, 695-701 (1984). 24. K. L. Leung, Dynamic analysis of plates using FDM-FEM and Laolace transform. M.S. Thesis. University of Minnesoti, Minneapolis, Minnesota (1982). 25. F. Durbin, Numerical inversion of Laplace transforms: an efficient improvement to Dubner’s and Abate’s method. Comput. J. 17, 371-376 (1974). 26. J. W. Cooley and J. W. Tukey, An algorithm for the machine calculation of complex Fourier series. Math. Comput. 19, 297-301 (1965): 27. B. A. Bolev and J. H. Weiner. Theorv of Thermal Stresses. Wiley, New York (196@. ’ 28. K. J. Bathe, E. L. Wilson and F. E. Peterson, SAP IV, a structural analysis program for static and dynamic response of linear systems. Rep. EERC 73-l 1, Univ. of California, Berkeley (June 1973).
APPENDIX Explicit expressions of the coefficients A’s appearing in the various finite difference patterns of Fig. 3 are given as follows[4,24]:
Al =6+6&+8a
+p
A5=a2
A2= -4(l+a)
A6=1
A3= -4a(l+a)
A7=5+6a2+8a+B
A4=2a
A9= -2(2a-va+l)
A8=a(2-v)
A10=1+4a(l-v)+3a2(l-v2)+B All=-2a(l-v+a(l-v’)) A12=ia2(l-v’)
(AlI
A13=5+5a2+8a+p Al4=2a(l Al5=
-v)
-2a(2-v+a)
Al6=1+4a(l-v)+ia’(l-v’)+/I
> Al8=2a(l-v)+i(l+a?(l-v2)+b
Al9=
-
a(, -v)+;(l
-v2) >
A20=;(1
-v’)
A21=6+5a2+8a+fl A22=a2+4a(l
-v)+3(1
-v2)+j
A23 = - 2(a(l - v) + (1 - v’)) A24=a2+4a(l-v)+i(l
-v*)+,9