- Email: [email protected]
Nonlinear dynamic analysis of stiffened plates
Compulers & S~ucrures Vol. 29, No. 6, pp, 929-941, Printed in Great Britain.
0045.7949/a 13.00 + 0.00 Pergamon Press plc
1988
NONLINEAR DYNAMIC ANALYSIS OF STIFFENED PLATES M. R. KHALIL, M. D. OLSON and D. L. ANDERSON Department of Civil Engineering, University of British Columbia, Vancouver, B.C., Canada V6T lW5 (Receic& 10 September 1987) Abstract-A new finite strip formulation for the nonlinear analysis of stiffened plate structures subjected to transient pressure loadings is presented. The effects of large deflections, and strain rate sensitive yielding material properties are included. An explicit central difference/diagonal mass matrix time stepping method is adopted. Example results are presented for an I-beam, an isotropic plate and a five-bay stiffened panel and compared with other predictions and/or experimental results. It is observed that design level accuracy can be obtained for practical structures for a fraction of the cost of full finite element analyses.
shown in Fig. l(b), are treated much as large finite elements with degrees of freedom located at certain Over the years, numerous studies have been directed nodal displacements. It is the choice of the displacetowards predicting the nonlinear dynamic response of ment functions (which will vary depending on plate-like structures subjected to intense pressure so as to enable one to write the displacements anyloading. The problem is time dependent and involves where in the middle plane of the strip in terms of the large deflections, high strain rates and material nodal displacements. It is in the choice of the disnonlinearity. Finite element [ 1,2] and finite differ- placement functions (which will vary depending on ence [3,4] programs have been developed for the boundary conditions at the ends of the strips) predicting the static, elasto-plastic, large deflection based on a priori knowledge of the displacement behaviour of beams, plates and stiffened panels [S, 61. shape, that allows the use of so few degrees of There are also finite element [7] and finite difference [8] freedom for such a large element. programs for predicting the dynamic response of members subjected to surface pressure loadings of 2.1. Displacement functions varying time histories. However, these finite element Finite strip displacement functions vary across the and difference methods were found to be very time strip in the same manner as finite elements, but the consuming and expensive to run, and require a large amount of input data unless mesh generators are Y available. The finite strip method can be used to approximate the response of plates and plate assemblages under blast type loads at a cost and effort significantly less than that required for finite element or finite difference analyses [9]. An added advantage is that such problems can be run on smaller computers. This paper describes the development of such a finite strip approach and the associated computer program. Many test examples have been run and some of the results are presented herein. More details (al Typical Strip Division for a Plate are available [lo]. The present work relies heavily on the static investigations by Abayakoon [l 11,which will be reported in a companion paper [ 121. I. INTRODUCTION
,r
2. THE FINITE STRIP
APPROACH
The finite strip approach [13] is generally used for members that have constant cross-sectional properties along one axis. Although it is not restricted to such that restriction is made in the following analyses. The member is divided into a number of parallel strips as shown in Fig. l(a). The individual strips,
(b)
Strip
Parameters
Fig. 1. Finite strip modelling of a plate 929
M. R. KHALIL et al.
930
variation along the length is given by continuous functions selected to satisfy the boundary conditions at the ends of the strips as well as provide a reasonably accurate description of the displacement. In addition they should satisfy compatibility requirements between adjacent strips. Letting the degrees of freedom be the U, v, w displacements and the rotation (0) about the x axis (i.e. 0 = c?w/dy) at each of the two nodal positions, the displacements for a strip are given by (see Fig.
2. Clamped ends, including axial restraint: ah(5 ) g:(r) =
m=l,3,5
at i + added mode sin 4~5
gX[)=sinnnt;
n=l,3,5,...
(3)
g;(5) = 4,,(r) = aP (sinh a,,5 - sin a,,lj)
1(b))
+ cash a,,5 - cos ap(;
u = c [(I - rl)uhl+ su*mlg:(5) m
p=l,3,5..., lJ =
c w - rl)U,”+ W,“k~(i;)
(1)
n
It’ = c [(1 - 3$ + 2tj3)w,, P
where c(,,= (cos ap - cash a,)/(sinh ap - sin a,,). The a,, are solutions of the transcendental equation, cash ap = set ap, giving a, = 4.7300407448, a3 = 10.995607838.
+ (q - 2q2 + ‘lf)bB,, + (3q2 - 2?/3)w,,
(ii) Members plates).
+ +I3 - v2)b~,,k;(r)?
with stiffeners (I-beams
1. Simply supported axial restraint:
where
Typically a,,,, represents the contribution to the U, degree of freedom of the mth mode. In many cases only one mode is sufficient, but for some displacements in some cases more than one mode is needed. In these situations
4 =
c UlW m
variable
m = 1,3,5
+ added mode sin 27~5
g:(t) = sin nnc;
n = 1, 3, 5
g;(s) = sin&;
p = 1,3,5.
(4)
2. Clamped ends, including axial restraint:
r
w,(r)
(W
or stiffened
ends for bending,
cosma{;
g:(5) =
.
-;
g:(5) =
m=l,3,5...
at
[ + added mode sin 471{
The variation of the displacements across the strip is chosen to satisfy compatibility, i.e. the variation is linear for u and v and cubic for w. Along the length the displacements vary according to the gL,gL, gi functions, which for different boundary conditions and for the assumption of loads symmetrical about the middle of strip are as follows. (i) Members without stiffeners (beams of rectangular cross-section or uniform plates). 1. Simply supported restrained:
g:(t)
= sin mn<;
ends for bending,
axially
m = 2,4,6 . .
,&(5)=sinnn~;
n = 1,3,5...
g;(r)=sinPt<;
p=l,3,5...
gXO=9n(5);
n=l,3,5...
g;(5)=+,(5);
p=l,3,5...,
(5)
where the 4, functions are the same as in the clamped case above, and are the linear vibration mode shapes for a uniform clamped beam. Stiffeners, such as shown in Fig. 2, can be modelled by considering the strip turned on its edge. In the same way I beams can be built up using a series of strips as shown later. For these cases the v displacement in the web stiffener must match the w displacement in connecting strips, and thus they must have the same functional form and magnitude. In general the displacements can be expressed as
(2) =N6 IWJ
(6)
931
Nonlinear dynamic analysis of stiffened plates
Fig. 3. Linear strain hardening. 2.3. Constitutive relations
Fig. 2.
Plate-stiffener assemblage.
by combining eqns (2)-(S), as appropriate, with eqn (1). 6 is the strip nodal displacement vector made up ofalltheu,,,v,,,w,,,e,, ,..., U*,,,,... terms,andas a minimum 6 will contain eight terms. The N matrix is made up of the appropriate g:(t), g:(t) and g:(c) functions in conjunction with the functions in q. 2.2. Strain-displacement relations For strips used to model rectangular beams or uniform plates the w displacement will be much larger than the u or u displacements. However, for finite strips used as web members in beams or stiffened plates the v displacement will be larger than the u or w displacements. The minimum number of nonlinear terms that adequately account for in-plane stretching are thus
The material behaviour is considered to be isotropic, elastic-plastic, and strain rate sensitive. For simplicity a bilinear model has been assumed with linear strain hardening as shown in Fig. 3. The effect of strain rate is included by using the Cowper-Symonds relation [ 141,
to calculate the effective yield stress I?,,,where o,, is the static yield stress. Typical parameters for mild steel are p = 5 and D = 40 sec. The von Mises yield criteria and associated flow rule have been adopted, along with an isotropic hardening model. Figures 4 and 5 show the
Ex=--zg+;[(q+(g)‘]
E.v_;_z$+#2J+(gl (a 1 Perfectly
azw awaw au au axay +aX& +Z$
ya = @ + fi - 22 -
ay
In those cases corresponding Substituting into the above
ax
(7)
where the u or w terms are small, the nonlinear terms may be dropped. the displacement functions eqn (1) gives, in symbolic form
=B6 +C,@),
plastic
surface (b)
Isotropic
strain
hardening
( C)
Kinematic
strain
hardening
(8)
where B represents the linear strain-displacement terms, while C,, which is a function of the S parameters, includes the contribution of the nonlinear terms in the strain-displacement relations eqns (7).
Fig. 4. Hardening models.
M. R. KHALIL et al.
932
increased. Upon unloading and stress reversal the yield again depends upon the strain rate but also upon the isotropic increase in the size due to the hardening effects. Following the above assumptions it is possible to determine the stress increment for a strain increment. This is written as da = D, dc,
(10)
where D, is the elastic-plastic constitutive matrix. It depends on the stress state, which determines the direction of the possible strain increments, the previous strain history which determines the amount of strain hardening, and the current strain rate. ABCDE ABCJK ABCHI ABCFG
-
Kinematic hardening Actual case (experimets) Independent yield Isotropic hordenlng
Fig. 5. Bauschinger effect representations difference in the isotropic and kinematic hardening models, with Fig. 4 depicting the yield stress in the usual principal stress space representation, while Fig. 5 shows the stress-strain relation for either a uniaxial test case or for the equivalent stress and strain definitions of more complex stress states. It is seen that isotropic hardening assumes that the yield surface remains centered at the origin but gets larger as the material hardens, i.e. is strained beyond the elastic limit. Independent hardening on the other hand assumes that the yield surface on reversal remains the same size as the original. For monotonic loading both models give the same result, while on reversal of loading a difference occurs. Figure 5 also shows an idealized experimental result and it is seen that the yield stress on reversal usually falls between the yield stress predicted by the kinematic and independent models (Bauschinger effect). However, the isotropic model is used here for simplicity. The von Mises yield criteria and associated flow law are well documented for strain rate insensitive or for quasi-static loading [15, 161. Figure 6(a) shows the equivalent stress-strain relation for a strain rate insensitive material. If at some time the strain and stress are given by point A, and during a time increment the strain is incremented by dc, then the stress must be increased by da so that the stress strain point falls on point B. For a strain rate sensitive material under high strain rates the strain rate model schematically shown in Fig. 6(b) has been assumed. Assume at some time the strain and strain rate are such that the stress is u2, (Fig. 6(b)), i.e. the stress strain point falls on point A. If during an increment of time the strain increases by the amount dc and the strain rate increases from iZ to i,, then the stress is taken to be bj corresponding to the stress point B. On the other hand if during the increment the strain rate decreases, say to i,, then the stress would be assumed to decrease to 0, (point C) even though the strain has
2.4. Virtual work The dynamic equilibrium equations are obtained via the principle of virtual work. Defining 8 as a virtual displacement of the nodal degrees of freedom vector, the resulting virtual strains P are given by 3 = [B + C(S)]&
(11)
where
is a linear function
of 6.
(a) Strain
rate
insensitive
Stroin
rate
sensitive
(b)
material
material
Fig, 6. Idealization of strain-rate behaviour for elasticplastic strain hardening material.
933
Nonlinear dynamic analysis of stiffened plates Using d’Alembert’s principle to include the inertia forces, the virtual work equation is
where p = mass per unit area V = strip volume ii = acceleration vector p = applied load vector A = strip area = uxb. Substituting
2.5. Integration points
eqns (6), (11) and (12) gives
pNrNdV%+ sY
NrpdA, IA
[B+C]‘adl’= sY
(13)
where 8 is the second time derivative of the nodal degrees of freedom. Defining M =
Equation (16) as derived forms the equation of motion for one strip. The strip equations can then be assembled into a set of equations for the entire structure, that would have the same form as eqn (16) except that now 6 and P would represent all the generalized degrees of freedoms and forces, M would be the assembled mass matrix and V would represent the entire volume of the structure. In the following, eqn (16) will be taken to represent the assembled structure equations.
pNrN dV = consistent mass matrix jY
(14)
The volume integrals in eqns (14) and (16) are evaluated numerically using a Gauss integration scheme. The number of integration points depends on the mode shape and whether the material behaviour is plastic or remains elastic. For the simplest elastic problem with a single half sine wave displacement, it was found that a scheme with 5-2-2 integration points in the x-y-z directions, respectively, was adequate. Figure 7 shows the location and the weights of the S-2-2 integration points. For other problems the number of integration points was increased appropriately. 2.6. Loading characteristics
Nrp dA = consistent
P=
load vector
s”
(15)
gives M8+
[B+Cj’udI’=P.
(16)
s”
Equation (16) represents the equation of motion in terms of the stresses in the strip. Since the stress-strain relations are nonlinear, these equations are solved using an incremental time step procedure. The C matrix, which is a function of the displacements, and the stresses u, must be updated at each time step.
For the analysis herein, the loading is considered to be a pressure uniformly applied to each strip, i.e. spatial variations are not considered, nor are forces in the x and y directions. The program is designed to read a pressure-time file so that any arbitrary time history can be considered. In addition, many blast pressure histories can be modelled fairly accurately by the following expression, which can be called by an option in the program. p(t)
=
p,,,(l
= 0
-
t/T)[“‘B’-‘l
O,
t >
T,
0.
0.04
0.19390 1 Y-
0.2113b
0.2113b
JJyiTq-
Fig. 7. Gauss numerical
integration
points, S-2-2 (two points weights: 0.4786, 0.5689, 0.4786, 0.2369).
1.0, 1.0; five points weights: 0.2369.
(17)
934
M. R. KHALIL et al.
where
been shown [ 191 that for stability P,,, = peak pressure
At < 2/w,,,,
duration of load p = decay parameter.
(21)
7 =
2.7. Time stepping procedure Equation (16) are solved using a numerical time stepping procedure. The method selected here is an explicit central difference scheme, which has been shown [17, 181 to be a very efficient method for computing transient response if the mass matrix is diagonalized. In the central difference scheme, the displacement at time t + At is given by
where w,,. is the largest frequency. For most of the problems considered here the dynamic motion of interest is in the w direction; the inertia terms associated with the u and u displacements are small and can be ignored. Since it also turns out that the maximum frequencies are associated with the u and u displacements, it has been found possible to assign large masses to the u and u degrees of freedom which increases the smallest periods and allows for a larger time step. 3. NUMERICAL APPLICATIONS
6 , + A,
=
26,- 4 -A,
+
(At)*&
(18)
where At = time step, and the subscripts denote time. Knowing S, and d,, conditions at time t can be calculated giving, from eqn (16) [B+C]ra,dI’ ”
1
(19)
which then gives 6, + &,from eqn (18). If the material is strain rate sensitive the velocity 8, +A, must also be determined. This is given by 6I+ Ari ,+A!=-.
4
3.1. Symmetrical
At
In the program the consistent mass matrix is diagonalized by ignoring the off-diagonal terms while multiplying the diagonal terms by a factor. This factor was determined by trial and error such as to give the correct frequency for the fundamental mode of vibration. To be computationally stable the time step must be a fraction of the smallest period of vibration. It has
495lb
A new finite strip program based on the foregoing theoretical formulation has recently been developed and tested on many examples. Some of these examples and the results obtained are described in the following. The first two applications are to relatively simple beam and plate problems for which other analytical or numerical results were available, and are intended to illustrate the basic capabilities of the method. The last two applications are to test panels for which experiments were performed by the Defence Research Establishment in Suffield, Alberta (DRES). I-beam
Five strips were used to model a symmetrical I-beam with clamped ends as shown in Fig. 8. All the displacements were free at the nodal line mid-points. The calculated fundamental frequency values are shown in Table 1, where the lumped mass matrix mentioned was taken as equal to 1.3 times the diagonal terms of the consistent mass matrix. The exact frequency was calculated from simple beam theory and while the finite strip predictions are quite
/in.
R
E = 30 x 106psi cO= 36 000
psi
Y = 0.0 p = 0.733x
10-31b.sec2/in4
At = 20 microseconds Fig. 8. Symmetric
I-beam
with clamped ends
configuration
and data.
935
Nonlinear dynamic analysis of stiffened plates Table 1. Fundamental frequency of clamped I-beam (rad/sec)
F.S.
F.S.
Exact
(Consistent M)
(Lumped M)
122.73
120.9
120.55
close, they are slightly on the low side. All calculations use the same exact beam vibration mode shape, so the differences must be due to shear deformation and rotary inertia effects which are included in the finite strip modelling and also, perhaps to a lesser extent, to numerical integration. This has now been confirmed by a beam finite element calculation including shear deformation and rotary inertia which gave a fundamental frequency of 120.0 rad/sec. The dynamic response of this I-beam was obtained for a step load in time assuming elastic-perfectly plastic material behaviour and nonlinear geometry. The load intensity and yield stress used are shown in Fig. 8. The central displacement response history obtained is shown in Fig. 9 along with that from the FENTAB program [20]. The second w mode referred to in the figure was sin 47rl as given in Sec. 2.1. The lumped mass matrix was used, and the time step for the finite strip calculation was 0.02msec. The two results compare quite well for most of the response, but the finite strip is stiff compared to the FENTAB model as the deflections are less and the peak displacement is reached sooner. However, the displacement levels predicted here are very extreme being about four times the beam depth. At these levels, the beam is acting essentially as a plastic string and it is known that the clamped beam mode used in the finite strip model is not optimum for this case [12]. Hence considering all this, an error of only 15% in the maximum displacement is really exceptionally good. 3.2. Simply supported square plate Four strips were used to model the simply supported square plate [called plate (I)] shown in Fig. 10. The calculated fundamental frequencies are shown in Table 2 along with the exact result. The lumped mass terms were taken as pab/4 and pab2/96 for w and .0, respectively. The predicted frequencies are very close to the exact one with the consistent mass one higher as expected.
E = p
205000
= 7.9x
10’
FENTAB\,/--, I
36
/
‘;; 32 2 .: 28 & 24 E % 20 z ,; 16
TWO U-MODES 10-2-2 GAUSS
4
0
0
12
24 Time,t
36
48
(msec
60
72
1
Fig. 9. Symmetric I-beam response to step load. For the transient calculations, the mass terms for the u and u degrees of freedom were multiplied by 100 to increase the minimum period and allow the use of a relatively large time step. The linear response to a step load of 2310 N/m2 was calculated with a time step of 0.025 msec and the result is shown in Fig. 11. The result is almost the same as the exact solution (cf. the peak displacement and its location shown Fig. 11). Figure 12 shows a comparison of the finite strip results with other numerical results [7] for a nonlinear, elastic-plastic response of a simply supported plate with edges constrained against in-plane motion (plate (2)) under a transient uniform pressure. This example is one of very large deflections with a central deflection of approximately 30 thicknesses or 0.4 times the span. Increasing the number of strips, and using two modes instead of one for each of the displacement functions reduced the stiffness of the finite strip model and made it closer to the other predicted results shown in Fig. 12.
MPo
kg/m3
v = 0.3 Fig. 10. Square plate configuration and data
936
KHALIL er al
M. R.
Peak-Elosttc,Smoll Osflcciion Theory
/2 Pre6ent
6”“‘$/ 0
/
,
0.5
Analysis
/ /,
1.0 Time,t
1.5 (msec
)
2.0
1
Fig. 11.Linear response of simply supported square plate to a step load. 0
’
’
0.4
’
’
0.8
Time,t(msec
3.3. DRES square plate
Fig. 12. Response
A square steel plate of constant thickness was subjected to an air blast loading in the DRES blast chamber facility [21]. The panel was instrumented with pressure, acceleration and strain transducers. The pressure measurements showed that the assumption of a spatially uniform pressure distribution over the plate was reasonable. Figure 13 shows the average pressure-time history and the plate particulars. The plate was modelled by four finite strips using a one mode approximation for each of the u, v and w displacements. The in-plane boundary conditions
F.S. (Consistent M)
F.S. (Lumped M)
3042.9
3079.9
3039.5
’
1
of simply supported load.
plate to intense blast
190
150
0
I
2
3 Time,
Fig. 13. Average
pressure
1
I.6
were fully constrained, while two different bending conditions were considered, namely, simply supported and clamped. The plate vibration frequencies obtained from the eigenvalue problems are shown in Tables 3 and 4, for the two different bending boundary conditions, respectively. The numbers in the mode column refer to the number of half waves across the plate in the direction perpendicular to the strips. The finite strip model only gives results for modes with one half wave in the strip direction because of the mode shape assumptions. Overall the predicted frequencies compare well with the theoretical ones. It was found from these eigenvalue problem results that with the use of large masses for the u and v degrees of freedom, a time step of 0.2 msec was adequate for the dynamic analysis for either boundary condition.
Table 2. Fundamental frequency for S.S. square plate (Plate (1)) Exact
’
12
for
4 5 t (msec 1
6
7
DRES plate (Shot No. 7).
8
937
Nonlinear dynamic analysis of stiffened plates Table 3. Frequencies of DRES simply supported square plate (rad/sec) Mode
F.S. (Consistent M)
F.S. (Lumped M)
Exact?
(l,l) (1.2) (1.3) (134) (IS)
402.66 1007.00 2016.80 4118.26 5753.46
397.41 952.36 1757.16 3362.60 3998.79
402.69 1006.73 2013.46 3422.89 5235.00
t v-21.
-III,
0
The nonlinear elastic-plastic response of the above model to the pressure load shown in Fig. 13 was calculated for both boundary condition cases. Hopefully, the experimental boundary conditions should lie somewhere between these two. Figure 14 shows a comparison of the experimental central displacement response with the two calculated results. The results for the simply supported case are reasonably close to the experiment. However, it appears that higher modes are contributing to the experimental results. To check this, the simply supported case was analysed again using a two mode approximation for each of the U, v and w displacements and using eight finite strips instead of four. The numerical integration was increased to 10-2-2 for the second model, while the normal 5-2-2 scheme was used for the first. Also the time step had to be reduced to 0.05 msec to avoid instability. The results of the two mode analysis are shown and compared with the experimental ones in Fig. 15. A modest improvement in the correlation of the results with experiment was achieved by the use of the second mode and more strips. No yielding was observed in any of the above calculations so the numerical results shown here are all strictly large deflection elastic.
r
I
08
1.6
2.4
3.2
I
4.0
T1me.t
48 (msec
56
6.4
7 2
I
8.0
)
Fig. 14. Response to DRES square plate (four strips, one mode). 3.4.1. Trial No. 314. The stiffened panel was suspended over a 96” x 180” opening, with the edges bolted to the foundations such that a relatively rigid connection was made between the plate and the foundation, simulating a clamped boundary. The panel was subjected to a blast resulting from the detonation of TNT explosives above the panel at the Height-of-Burst (HOB) site located at DRES [21]. The pressure transducers were used to establish the average pressure-time history of the air blast which is shown in Fig. 17. The square symbols on the plot show the data points used for the piecewise linear representation used in the program. At first 11 finite strips were used to model half of the stiffened panel, with symmetry boundary conditions used at the centre line. As shown in Fig. 18
13r
3.4. DRES stiffened panel Several 8’ x 15’ T-beam stiffened steel panels with l/4” plating and 3” x 6” T-beam stiffeners were tested at DRES. Figure 16 shows the details for these panels. Several different load cases were considered at DRES but only two are presented in the following along with their corresponding finite strip analyses. Note that the same time step of 5 psec was used for all the runs described in the following.
-III0
’
08
’
1.6
’
24
’
3.2
’
4.0
’
5.6
’
64
’
7.2
’
8.0
Tlme,t (msec 1 Fig. 15. Response of DRES square plate with extra modes and strips (simply supported in bending).
Table 4. Frequencies of DRES clamped square plate (rad/sec) Mode
Upper bound?
Lower bound?
F.S. (Consistent M)
F.S. (Lumped M)
(171) (192) (1,3)--(3,l) (1.3) + (3.1) (1.4) (175)- (591) (ITS)+ (531)
734.13 1497.28 2684.33 2697.06 4294.78 6301.81 6307.15
734.12 1497.25 2684.18 2696.90 4288.46 6290.70 6295.56
734.96 1497.19 2596.20
813.42 1576.21 2417.53
5224.32 7659.20
5777.45 7005.13
t ~231.
’
4.8
938
M. R. KHALILetal.
60
-1,; I“!‘--/;/0
5
IO
15
20 Time
25 ,t (msec
30
35
40
45
1
Fig. 17. Average pressure for DRES panel (HOB 314).
YOUNG’S
MODULUS,
POISSON’S
RATIO,
YIELD STRAIN
STRESS PER
TIME
STEP,
=
UNIT
X IO’psi
0.3
IN SIMPLE
HARDENING
MASS
E * 0.3 v
with six strips between stiffeners amounting to a total of 21 strips. In all three models, a one mode approximation was used for each of v and w displacements and a two mode approximation was used for the u displacement. Also a scheme of 5-2-2 Gauss points was used in all three models. Figure 19 shows a comparison of the panel deflection responses for the three models. The convergence with strip refinement is quite apparent. The results that follow were calculated with the 21 strip model. Figure 20 shows the deflection at the centre of the central panel (D12) and at the midpoint of the
TENSION,
MODULUS, VOLUME
ET=
, p =0.733
0.544X 0.178 x IO-’
IO’Psi
./e
X106PSi lb.sec%n4
At - 5 microseconds
Fig. 16. DRES stiffened panel geometry and data
only two strips were used for modelling the portion of the plate between adjacent stiffeners. The eigenvalue calculation results for this model are shown in Table 5 along with some computed by the VAST finite element analysis package [24]. The first few frequencies compare reasonably well. The nonlinear elastic-plastic response of this model to the dynamic loading was significantly lower than expected. Hence, two more models were analysed, one with four strips between stiffeners amounting to a total of 16 strips for the model and the other
x,u
/’ Fig. 18. Eleven strip model of DRES stiffened panel.
Table 5. DRES stiffened panel frequencies (Hz) Mode
VAST7
F.S. (Consistent M)
F.S. (Lumped M)
40.9 41.21 41.86
40.33 41.7 42.54
40.97 42.12 42.81
Mainly panel modes
Mode (1) Mode (2) Mode (3)
Mainly stiffener modes
Mode (4) Mode (5)
-
105.13 109.42
117.64 119.78
Mixed modes
Mode (6) Mode (7)
-
143.92 156.77
155.7 158.54
t [24].
939
Nonlinear dynamic analysis of stiffened plates
-2.0
2.0
6.0
14.0
10.0
12.0
22.0
-2.0
2.0
6.0
TIME(MS) (a)
Panel
10.0
Ii.0
Ii.0
25.0
TIME(MS)
- Dl2 ..-.-. ---
(a) I I strip model 16 strip model 21 stripmodel
Pane
I - Dl2
89 _( -2s I -.
-0 u t-i41 $9 z2 Ptf -2.0
2.0
8.0
10.0
14.0
18.0
!.O
TIME(MS) (b) -2.0
i.0
60
lb.0
Ii.0
lb.0
2h.o
TIME(MS) (b)
Stiffener
-
- D13
Fig. 20. DRES stiffened panel response to HOB 314: predicted and experimental displacements.
Dl3
Fig. 19. DRES stiffened panel response to HOB 314: effect of strip refinement.
adjacent stiffener (D13) along with the experimental results for comparison. The predicted and experimental results agree well up to 2 msec but beyond this
time the finite strip model appears to quickly become stiffer than the experimental panel. It should be mentioned here that the experimental deflections were not measured directly, but were calculated from double integration of the accelerations. There are some factors that may explain the difference in the results. The first factor is that only a few modes were used to approximate the displacements. The second factor is the boundary condition of the real plate, which may not be exactly clamped or fully restrained in-plane. Some slippage at the boundary bolts is also a possibility. 3.4.2. Minor scale. A stiffened panel identical to the one described above was mounted flush with the ground surface and subjected to an air blast load from a ground explosion. The panel was embedded in a reinforced concrete foundation to simulate clamped boundary conditions. The pressure-time history shown in Fig. 21 was obtained by averaging the measurements of several pressure transducers after synchronizing the rise times. The following
Stiffener
analysis was simplified
was found in the Trial No. 314 calculations that the improvement in response by the use of the 21 strip model instead of the 16 strip model was not significant. Therefore, only the 11 strip model and the 16 strip model were analysed in this example. The central displacements at D12 and D13 locations of both the 11 strip model and the 16 strip model are shown in Fig. 22. The figure shows some higher mode contributions in the response of the 16 strip model. Also at this load level the plate panel has deflected more than the stiffener as was the case in Trial No. 314. The results from the 16 strip mode1 are compared to the experimentally measured displacements in Figs 23 and 24. It is seen that the predicted and measured deflections are quite close for the first 14msec but thereafter the experimental results are higher than the finite strip ones, especially for the stiffener. Some of
__ 40 h
by ignor-
ing the transit time (about 3 msec) of the shock front over the panel (parallel to the stiffeners), and assuming the pressure loading to be spatially uniform over the panel surface. A nonlinear elastic-plastic dynamic analysis was performed using the above pressure-time history. It
IO
I J
0
IO
20
30 Time.t (msec)
40
50
Fig. 21. Average pressure for minor scale test
60
M. R. KHALILet
al.
and for the I-beam example. These results are rather preliminary and more detailed comparisons still need to be carried out. The 16 strip model used about 14 min of CPU time on an Amdahl 5850 for the 50msec calculation shown here. 4. CONCLUSIONS
-2
6
14
22 Time,
30
38
46
t ( msec)
Fig. 22. Response of DRES stiffened panel (minor scale) using 11 and 16 strips.
this discrepancy may be due to the approximations made in the pressure representation. However the fact that the numerical results are on the stiff side is
consistent with the previous results for the 314 test
The finite strip method has been applied to the large deflection, elastic-plastic dynamic response of plate-type structures to air-blast loads. The loads have been modelled as uniformly distributed time dependent pressures. The method is ideally suited to this type of load because the distribution can be well represented by a single mode. The method has been applied to a wide variety of configurations including rectangular beams, symmetric and asymmetric Ibeams, isotropic plates and stiffened plates. The results reported herein are very encouraging. It appears that reasonable engineering accuracy can
4CENTRE
PANEL
-56# -6
2
IO
18 26 Tlme,t(msec)
34
42
50
Fig. 23. Response of DRES stiffened panel (minor scale) using 16 strips-D12.
0L ;-24 s
-
-6
2
IO
I8
26
34
42
50
Tlme,t(msec)
Fig. 24. Response of DRES stiffened panel (minor scale) using 16 strips-Dl3.
Nonlinear dynamic analysis of stiffened plates be achieved by modelling involving only a few modes and a few finite strips. The number of variables and the computer times required are typically an order of magnitude less than those required by corresponding finite element codes. Thus considerable savings in cost are possible for certain problems with only a minimal loss in accuracy. The computer code is also much smaller and can be put on a micro-computer. The method seems ideally suited for preliminary design work. So far the main limitation seems to be in the application to the case of clamped boundary conditions for a structure which is highly efficient in bending such as a symmetric I-beam. The plastic deformation of such a structure is not modelled well by the clamped beam vibration mode. This is in contrast to the simply supported case which apparently is represented quite well by a half sine wave. The former problem needs further work, perhaps in finding a better mode to represent it. On the other hand, the stiffened panel with clamped boundaries seems to be adequately represented by the aforementioned beam modes. In particular, the displacement predictions for the DRES stiffened panel (Minor Scale test) as reported herein compare quite well with the test results. The combination of a diagonalized mass matrix and the finite difference time stepping method seems ideally suited to the present applications. The requirement of small time steps for stability is very compatible with the need for updating strains and stresses to follow yielding and plastic deformation. Further it is also needed for those cases with rapidly changing pressure-time history such as in the above Minor Scale test.
941
On the behaviour and design of stiffened plates in ultimate limit state. J Ship Res. 22, 238-244 (1978). 7. T. R. Stagliano and L. J. Mente, Large deflection, elastic-plastic dynamic structural response of beams and stiffened or unstiffened panels-a comparison of finite element, finite difference and modal solutions. In Nonlinear Finite Element Analysis and ADINA, Proc. ADINA Conference (Edited by K. J. Bathe et al.) (1979). 8. R. W. H. Wu and E. A. Witmer, Nonlinear transient responses of structures by the spatial finite element method. AIAA J. 11 (1973). 9. D. S. Mofflin, M. D. Olson and D. L. Anderson, Finite strip analysis of blast loaded plates. In Finite Element Methods for Nonlinear Problems (Edited by P. G. Bergan et al.), pp. 539-554. Springer, Berlin (1986). 10. R. Khalil, M. D. Olson and D. L. Anderson, Large deflection, elastic-plastic dynamic response of air-blast loaded plate structures by the finite strip method. Strut. Series Rept. No. 33, Dept. of Civil Engineering, UBC (1987). 11. S. Abayakoon, Large deflection, elastic-plastic analysis of plate structures by the finite strip method. Ph.D. thesis, Dept. of Civil Engineering, UBC, Vancouver (1987). 12. S. B. S. Abayakoon, M. D. Olson and D. L. Anderson, Large deflection elastic-plastic analysis of plate structures by the finite strip method. Int. J. Numer. Meth. Engng (in press). 13. Y. K. Cheung, Finite Strip Method in Structural Analysis. Pergamon Press, New York (1976). 14. S. R. Bonder and P. S. Symonds, Experimental and theoretical investigation of the plastic deformation of cantilever beams subjected to impulsive loading. 6. T. H. Soreide, T. Moan and N. T. Nordsve.
ASME,
J. appl. Mech. 29, 719-728 (1962).
15. A. Mendelson, Plasticity Theory and Application. MacMillan New York (1968). 16. L. M. Kachanov, Fundamentals of the Theory of Plasticity. [Translated From the Russian by M. Konyaeva, Mir Publishers. Moscow (1974).1 17. R. D. Krieg and S. W. Key, Transient shell response by numerical time integration. Int. J. Numer. Meth. Engng 7, 273-286 (1973).
Acknowledgement-This
research was supported by the Canadian Department of National Defence though a contract from the Defence Research Establishment Suffield.
REFERENCES
1. M. A. Chrisfield, The automatic nonlinear analysis of stiffened plates and shallow shells using finite elements. Proc. Inst. Cin. Engng, 69 (1980). 2. E. Ramm., Geometrisch-nichtlineare elastostatik und finite elemente. University of Stuttgart, Report 76-2 (1976). 3. P. A. Frieze, P. J. Dowling and R. E. Hobbs, Ultimate load behaviour of plates in compression. Proc. of Conf. on Steel Plated Structures, London. Crosby Lockwood Staples, London (1976). 4. J. E. Harding, P. E. Hobbs and B. G. Neal, The elasto-plastic analysis of imperfect square plates under in-plane loading. Proc. Inst. Cio. Engng, 63 (1977). 5. S. E. Webb and P. J. Dowling, Large deflection elastoplastic behaviour of discretely stiffened plates. Proc. Inst. Cia. Engng. 69 (1980).
18. R. D. Cook, Concepts and Application of Finite Element Analysis, 2nd edn. John Wiley, New York (1981). 19. J. W. Leech, P. Hsu and E. W. Mack, Stability of a finite difference method for solving matrix equations. AIAA Jnl3 (ll), 2172-2173 (1965). 20. B. R. Folz, Numerical simulation of the nonlinear transient response of slender beams. M.A.Sc. thesis, U.B.C., Vancouver, British Columbia (1986). 21. R. Houlston, D. Ritzel, J. E. Slater, G. Rude and R. T. Schmitke, Air blast experiments on square plates. DRES Suffield Memorandum No. 114153, Project No. 27C77, Defence Research Establishment Suffield, Ralston, Alberta (1986). 22. G. R. Cowper, E. Kosko, G. M. Lindberg and M. D. Olson, A high precision triangular plate-bending element. National Research Council of Canada, Aeronautical Report LR-5 14, Ottawa (1968). 23. L. De Vito, G. Fichera, A. Fusciardi and M. Schaerf, Sul Calcolo degli Autovalori della Piastra Quadrate Incastrata Lungo il Bordo. Rend. Accad. Nazion. Lincei, 8 Ser. 40, 725-733 (1966).
24. R. Houlston, J. E. Slater, N. Pegg and C. Desrochers, Structural response of ship panels subjected to air blast loading. 5th ADINA Conference, June 1985.
0045.7949/a 13.00 + 0.00 Pergamon Press plc
1988
NONLINEAR DYNAMIC ANALYSIS OF STIFFENED PLATES M. R. KHALIL, M. D. OLSON and D. L. ANDERSON Department of Civil Engineering, University of British Columbia, Vancouver, B.C., Canada V6T lW5 (Receic& 10 September 1987) Abstract-A new finite strip formulation for the nonlinear analysis of stiffened plate structures subjected to transient pressure loadings is presented. The effects of large deflections, and strain rate sensitive yielding material properties are included. An explicit central difference/diagonal mass matrix time stepping method is adopted. Example results are presented for an I-beam, an isotropic plate and a five-bay stiffened panel and compared with other predictions and/or experimental results. It is observed that design level accuracy can be obtained for practical structures for a fraction of the cost of full finite element analyses.
shown in Fig. l(b), are treated much as large finite elements with degrees of freedom located at certain Over the years, numerous studies have been directed nodal displacements. It is the choice of the displacetowards predicting the nonlinear dynamic response of ment functions (which will vary depending on plate-like structures subjected to intense pressure so as to enable one to write the displacements anyloading. The problem is time dependent and involves where in the middle plane of the strip in terms of the large deflections, high strain rates and material nodal displacements. It is in the choice of the disnonlinearity. Finite element [ 1,2] and finite differ- placement functions (which will vary depending on ence [3,4] programs have been developed for the boundary conditions at the ends of the strips) predicting the static, elasto-plastic, large deflection based on a priori knowledge of the displacement behaviour of beams, plates and stiffened panels [S, 61. shape, that allows the use of so few degrees of There are also finite element [7] and finite difference [8] freedom for such a large element. programs for predicting the dynamic response of members subjected to surface pressure loadings of 2.1. Displacement functions varying time histories. However, these finite element Finite strip displacement functions vary across the and difference methods were found to be very time strip in the same manner as finite elements, but the consuming and expensive to run, and require a large amount of input data unless mesh generators are Y available. The finite strip method can be used to approximate the response of plates and plate assemblages under blast type loads at a cost and effort significantly less than that required for finite element or finite difference analyses [9]. An added advantage is that such problems can be run on smaller computers. This paper describes the development of such a finite strip approach and the associated computer program. Many test examples have been run and some of the results are presented herein. More details (al Typical Strip Division for a Plate are available [lo]. The present work relies heavily on the static investigations by Abayakoon [l 11,which will be reported in a companion paper [ 121. I. INTRODUCTION
,r
2. THE FINITE STRIP
APPROACH
The finite strip approach [13] is generally used for members that have constant cross-sectional properties along one axis. Although it is not restricted to such that restriction is made in the following analyses. The member is divided into a number of parallel strips as shown in Fig. l(a). The individual strips,
(b)
Strip
Parameters
Fig. 1. Finite strip modelling of a plate 929
M. R. KHALIL et al.
930
variation along the length is given by continuous functions selected to satisfy the boundary conditions at the ends of the strips as well as provide a reasonably accurate description of the displacement. In addition they should satisfy compatibility requirements between adjacent strips. Letting the degrees of freedom be the U, v, w displacements and the rotation (0) about the x axis (i.e. 0 = c?w/dy) at each of the two nodal positions, the displacements for a strip are given by (see Fig.
2. Clamped ends, including axial restraint: ah(5 ) g:(r) =
m=l,3,5
at i + added mode sin 4~5
gX[)=sinnnt;
n=l,3,5,...
(3)
g;(5) = 4,,(r) = aP (sinh a,,5 - sin a,,lj)
1(b))
+ cash a,,5 - cos ap(;
u = c [(I - rl)uhl+ su*mlg:(5) m
p=l,3,5..., lJ =
c w - rl)U,”+ W,“k~(i;)
(1)
n
It’ = c [(1 - 3$ + 2tj3)w,, P
where c(,,= (cos ap - cash a,)/(sinh ap - sin a,,). The a,, are solutions of the transcendental equation, cash ap = set ap, giving a, = 4.7300407448, a3 = 10.995607838.
+ (q - 2q2 + ‘lf)bB,, + (3q2 - 2?/3)w,,
(ii) Members plates).
+ +I3 - v2)b~,,k;(r)?
with stiffeners (I-beams
1. Simply supported axial restraint:
where
Typically a,,,, represents the contribution to the U, degree of freedom of the mth mode. In many cases only one mode is sufficient, but for some displacements in some cases more than one mode is needed. In these situations
4 =
c UlW m
variable
m = 1,3,5
+ added mode sin 27~5
g:(t) = sin nnc;
n = 1, 3, 5
g;(s) = sin&;
p = 1,3,5.
(4)
2. Clamped ends, including axial restraint:
r
w,(r)
(W
or stiffened
ends for bending,
cosma{;
g:(5) =
.
-;
g:(5) =
m=l,3,5...
at
[ + added mode sin 471{
The variation of the displacements across the strip is chosen to satisfy compatibility, i.e. the variation is linear for u and v and cubic for w. Along the length the displacements vary according to the gL,gL, gi functions, which for different boundary conditions and for the assumption of loads symmetrical about the middle of strip are as follows. (i) Members without stiffeners (beams of rectangular cross-section or uniform plates). 1. Simply supported restrained:
g:(t)
= sin mn<;
ends for bending,
axially
m = 2,4,6 . .
,&(5)=sinnn~;
n = 1,3,5...
g;(r)=sinPt<;
p=l,3,5...
gXO=9n(5);
n=l,3,5...
g;(5)=+,(5);
p=l,3,5...,
(5)
where the 4, functions are the same as in the clamped case above, and are the linear vibration mode shapes for a uniform clamped beam. Stiffeners, such as shown in Fig. 2, can be modelled by considering the strip turned on its edge. In the same way I beams can be built up using a series of strips as shown later. For these cases the v displacement in the web stiffener must match the w displacement in connecting strips, and thus they must have the same functional form and magnitude. In general the displacements can be expressed as
(2) =N6 IWJ
(6)
931
Nonlinear dynamic analysis of stiffened plates
Fig. 3. Linear strain hardening. 2.3. Constitutive relations
Fig. 2.
Plate-stiffener assemblage.
by combining eqns (2)-(S), as appropriate, with eqn (1). 6 is the strip nodal displacement vector made up ofalltheu,,,v,,,w,,,e,, ,..., U*,,,,... terms,andas a minimum 6 will contain eight terms. The N matrix is made up of the appropriate g:(t), g:(t) and g:(c) functions in conjunction with the functions in q. 2.2. Strain-displacement relations For strips used to model rectangular beams or uniform plates the w displacement will be much larger than the u or u displacements. However, for finite strips used as web members in beams or stiffened plates the v displacement will be larger than the u or w displacements. The minimum number of nonlinear terms that adequately account for in-plane stretching are thus
The material behaviour is considered to be isotropic, elastic-plastic, and strain rate sensitive. For simplicity a bilinear model has been assumed with linear strain hardening as shown in Fig. 3. The effect of strain rate is included by using the Cowper-Symonds relation [ 141,
to calculate the effective yield stress I?,,,where o,, is the static yield stress. Typical parameters for mild steel are p = 5 and D = 40 sec. The von Mises yield criteria and associated flow rule have been adopted, along with an isotropic hardening model. Figures 4 and 5 show the
Ex=--zg+;[(q+(g)‘]
E.v_;_z$+#2J+(gl (a 1 Perfectly
azw awaw au au axay +aX& +Z$
ya = @ + fi - 22 -
ay
In those cases corresponding Substituting into the above
ax
(7)
where the u or w terms are small, the nonlinear terms may be dropped. the displacement functions eqn (1) gives, in symbolic form
=B6 +C,@),
plastic
surface (b)
Isotropic
strain
hardening
( C)
Kinematic
strain
hardening
(8)
where B represents the linear strain-displacement terms, while C,, which is a function of the S parameters, includes the contribution of the nonlinear terms in the strain-displacement relations eqns (7).
Fig. 4. Hardening models.
M. R. KHALIL et al.
932
increased. Upon unloading and stress reversal the yield again depends upon the strain rate but also upon the isotropic increase in the size due to the hardening effects. Following the above assumptions it is possible to determine the stress increment for a strain increment. This is written as da = D, dc,
(10)
where D, is the elastic-plastic constitutive matrix. It depends on the stress state, which determines the direction of the possible strain increments, the previous strain history which determines the amount of strain hardening, and the current strain rate. ABCDE ABCJK ABCHI ABCFG
-
Kinematic hardening Actual case (experimets) Independent yield Isotropic hordenlng
Fig. 5. Bauschinger effect representations difference in the isotropic and kinematic hardening models, with Fig. 4 depicting the yield stress in the usual principal stress space representation, while Fig. 5 shows the stress-strain relation for either a uniaxial test case or for the equivalent stress and strain definitions of more complex stress states. It is seen that isotropic hardening assumes that the yield surface remains centered at the origin but gets larger as the material hardens, i.e. is strained beyond the elastic limit. Independent hardening on the other hand assumes that the yield surface on reversal remains the same size as the original. For monotonic loading both models give the same result, while on reversal of loading a difference occurs. Figure 5 also shows an idealized experimental result and it is seen that the yield stress on reversal usually falls between the yield stress predicted by the kinematic and independent models (Bauschinger effect). However, the isotropic model is used here for simplicity. The von Mises yield criteria and associated flow law are well documented for strain rate insensitive or for quasi-static loading [15, 161. Figure 6(a) shows the equivalent stress-strain relation for a strain rate insensitive material. If at some time the strain and stress are given by point A, and during a time increment the strain is incremented by dc, then the stress must be increased by da so that the stress strain point falls on point B. For a strain rate sensitive material under high strain rates the strain rate model schematically shown in Fig. 6(b) has been assumed. Assume at some time the strain and strain rate are such that the stress is u2, (Fig. 6(b)), i.e. the stress strain point falls on point A. If during an increment of time the strain increases by the amount dc and the strain rate increases from iZ to i,, then the stress is taken to be bj corresponding to the stress point B. On the other hand if during the increment the strain rate decreases, say to i,, then the stress would be assumed to decrease to 0, (point C) even though the strain has
2.4. Virtual work The dynamic equilibrium equations are obtained via the principle of virtual work. Defining 8 as a virtual displacement of the nodal degrees of freedom vector, the resulting virtual strains P are given by 3 = [B + C(S)]&
(11)
where
is a linear function
of 6.
(a) Strain
rate
insensitive
Stroin
rate
sensitive
(b)
material
material
Fig, 6. Idealization of strain-rate behaviour for elasticplastic strain hardening material.
933
Nonlinear dynamic analysis of stiffened plates Using d’Alembert’s principle to include the inertia forces, the virtual work equation is
where p = mass per unit area V = strip volume ii = acceleration vector p = applied load vector A = strip area = uxb. Substituting
2.5. Integration points
eqns (6), (11) and (12) gives
pNrNdV%+ sY
NrpdA, IA
[B+C]‘adl’= sY
(13)
where 8 is the second time derivative of the nodal degrees of freedom. Defining M =
Equation (16) as derived forms the equation of motion for one strip. The strip equations can then be assembled into a set of equations for the entire structure, that would have the same form as eqn (16) except that now 6 and P would represent all the generalized degrees of freedoms and forces, M would be the assembled mass matrix and V would represent the entire volume of the structure. In the following, eqn (16) will be taken to represent the assembled structure equations.
pNrN dV = consistent mass matrix jY
(14)
The volume integrals in eqns (14) and (16) are evaluated numerically using a Gauss integration scheme. The number of integration points depends on the mode shape and whether the material behaviour is plastic or remains elastic. For the simplest elastic problem with a single half sine wave displacement, it was found that a scheme with 5-2-2 integration points in the x-y-z directions, respectively, was adequate. Figure 7 shows the location and the weights of the S-2-2 integration points. For other problems the number of integration points was increased appropriately. 2.6. Loading characteristics
Nrp dA = consistent
P=
load vector
s”
(15)
gives M8+
[B+Cj’udI’=P.
(16)
s”
Equation (16) represents the equation of motion in terms of the stresses in the strip. Since the stress-strain relations are nonlinear, these equations are solved using an incremental time step procedure. The C matrix, which is a function of the displacements, and the stresses u, must be updated at each time step.
For the analysis herein, the loading is considered to be a pressure uniformly applied to each strip, i.e. spatial variations are not considered, nor are forces in the x and y directions. The program is designed to read a pressure-time file so that any arbitrary time history can be considered. In addition, many blast pressure histories can be modelled fairly accurately by the following expression, which can be called by an option in the program. p(t)
=
p,,,(l
= 0
-
t/T)[“‘B’-‘l
O,
t >
T,
0.
0.04
0.19390 1 Y-
0.2113b
0.2113b
JJyiTq-
Fig. 7. Gauss numerical
integration
points, S-2-2 (two points weights: 0.4786, 0.5689, 0.4786, 0.2369).
1.0, 1.0; five points weights: 0.2369.
(17)
934
M. R. KHALIL et al.
where
been shown [ 191 that for stability P,,, = peak pressure
At < 2/w,,,,
duration of load p = decay parameter.
(21)
7 =
2.7. Time stepping procedure Equation (16) are solved using a numerical time stepping procedure. The method selected here is an explicit central difference scheme, which has been shown [17, 181 to be a very efficient method for computing transient response if the mass matrix is diagonalized. In the central difference scheme, the displacement at time t + At is given by
where w,,. is the largest frequency. For most of the problems considered here the dynamic motion of interest is in the w direction; the inertia terms associated with the u and u displacements are small and can be ignored. Since it also turns out that the maximum frequencies are associated with the u and u displacements, it has been found possible to assign large masses to the u and u degrees of freedom which increases the smallest periods and allows for a larger time step. 3. NUMERICAL APPLICATIONS
6 , + A,
=
26,- 4 -A,
+
(At)*&
(18)
where At = time step, and the subscripts denote time. Knowing S, and d,, conditions at time t can be calculated giving, from eqn (16) [B+C]ra,dI’ ”
1
(19)
which then gives 6, + &,from eqn (18). If the material is strain rate sensitive the velocity 8, +A, must also be determined. This is given by 6I+ Ari ,+A!=-.
4
3.1. Symmetrical
At
In the program the consistent mass matrix is diagonalized by ignoring the off-diagonal terms while multiplying the diagonal terms by a factor. This factor was determined by trial and error such as to give the correct frequency for the fundamental mode of vibration. To be computationally stable the time step must be a fraction of the smallest period of vibration. It has
495lb
A new finite strip program based on the foregoing theoretical formulation has recently been developed and tested on many examples. Some of these examples and the results obtained are described in the following. The first two applications are to relatively simple beam and plate problems for which other analytical or numerical results were available, and are intended to illustrate the basic capabilities of the method. The last two applications are to test panels for which experiments were performed by the Defence Research Establishment in Suffield, Alberta (DRES). I-beam
Five strips were used to model a symmetrical I-beam with clamped ends as shown in Fig. 8. All the displacements were free at the nodal line mid-points. The calculated fundamental frequency values are shown in Table 1, where the lumped mass matrix mentioned was taken as equal to 1.3 times the diagonal terms of the consistent mass matrix. The exact frequency was calculated from simple beam theory and while the finite strip predictions are quite
/in.
R
E = 30 x 106psi cO= 36 000
psi
Y = 0.0 p = 0.733x
10-31b.sec2/in4
At = 20 microseconds Fig. 8. Symmetric
I-beam
with clamped ends
configuration
and data.
935
Nonlinear dynamic analysis of stiffened plates Table 1. Fundamental frequency of clamped I-beam (rad/sec)
F.S.
F.S.
Exact
(Consistent M)
(Lumped M)
122.73
120.9
120.55
close, they are slightly on the low side. All calculations use the same exact beam vibration mode shape, so the differences must be due to shear deformation and rotary inertia effects which are included in the finite strip modelling and also, perhaps to a lesser extent, to numerical integration. This has now been confirmed by a beam finite element calculation including shear deformation and rotary inertia which gave a fundamental frequency of 120.0 rad/sec. The dynamic response of this I-beam was obtained for a step load in time assuming elastic-perfectly plastic material behaviour and nonlinear geometry. The load intensity and yield stress used are shown in Fig. 8. The central displacement response history obtained is shown in Fig. 9 along with that from the FENTAB program [20]. The second w mode referred to in the figure was sin 47rl as given in Sec. 2.1. The lumped mass matrix was used, and the time step for the finite strip calculation was 0.02msec. The two results compare quite well for most of the response, but the finite strip is stiff compared to the FENTAB model as the deflections are less and the peak displacement is reached sooner. However, the displacement levels predicted here are very extreme being about four times the beam depth. At these levels, the beam is acting essentially as a plastic string and it is known that the clamped beam mode used in the finite strip model is not optimum for this case [12]. Hence considering all this, an error of only 15% in the maximum displacement is really exceptionally good. 3.2. Simply supported square plate Four strips were used to model the simply supported square plate [called plate (I)] shown in Fig. 10. The calculated fundamental frequencies are shown in Table 2 along with the exact result. The lumped mass terms were taken as pab/4 and pab2/96 for w and .0, respectively. The predicted frequencies are very close to the exact one with the consistent mass one higher as expected.
E = p
205000
= 7.9x
10’
FENTAB\,/--, I
36
/
‘;; 32 2 .: 28 & 24 E % 20 z ,; 16
TWO U-MODES 10-2-2 GAUSS
4
0
0
12
24 Time,t
36
48
(msec
60
72
1
Fig. 9. Symmetric I-beam response to step load. For the transient calculations, the mass terms for the u and u degrees of freedom were multiplied by 100 to increase the minimum period and allow the use of a relatively large time step. The linear response to a step load of 2310 N/m2 was calculated with a time step of 0.025 msec and the result is shown in Fig. 11. The result is almost the same as the exact solution (cf. the peak displacement and its location shown Fig. 11). Figure 12 shows a comparison of the finite strip results with other numerical results [7] for a nonlinear, elastic-plastic response of a simply supported plate with edges constrained against in-plane motion (plate (2)) under a transient uniform pressure. This example is one of very large deflections with a central deflection of approximately 30 thicknesses or 0.4 times the span. Increasing the number of strips, and using two modes instead of one for each of the displacement functions reduced the stiffness of the finite strip model and made it closer to the other predicted results shown in Fig. 12.
MPo
kg/m3
v = 0.3 Fig. 10. Square plate configuration and data
936
KHALIL er al
M. R.
Peak-Elosttc,Smoll Osflcciion Theory
/2 Pre6ent
6”“‘$/ 0
/
,
0.5
Analysis
/ /,
1.0 Time,t
1.5 (msec
)
2.0
1
Fig. 11.Linear response of simply supported square plate to a step load. 0
’
’
0.4
’
’
0.8
Time,t(msec
3.3. DRES square plate
Fig. 12. Response
A square steel plate of constant thickness was subjected to an air blast loading in the DRES blast chamber facility [21]. The panel was instrumented with pressure, acceleration and strain transducers. The pressure measurements showed that the assumption of a spatially uniform pressure distribution over the plate was reasonable. Figure 13 shows the average pressure-time history and the plate particulars. The plate was modelled by four finite strips using a one mode approximation for each of the u, v and w displacements. The in-plane boundary conditions
F.S. (Consistent M)
F.S. (Lumped M)
3042.9
3079.9
3039.5
’
1
of simply supported load.
plate to intense blast
190
150
0
I
2
3 Time,
Fig. 13. Average
pressure
1
I.6
were fully constrained, while two different bending conditions were considered, namely, simply supported and clamped. The plate vibration frequencies obtained from the eigenvalue problems are shown in Tables 3 and 4, for the two different bending boundary conditions, respectively. The numbers in the mode column refer to the number of half waves across the plate in the direction perpendicular to the strips. The finite strip model only gives results for modes with one half wave in the strip direction because of the mode shape assumptions. Overall the predicted frequencies compare well with the theoretical ones. It was found from these eigenvalue problem results that with the use of large masses for the u and v degrees of freedom, a time step of 0.2 msec was adequate for the dynamic analysis for either boundary condition.
Table 2. Fundamental frequency for S.S. square plate (Plate (1)) Exact
’
12
for
4 5 t (msec 1
6
7
DRES plate (Shot No. 7).
8
937
Nonlinear dynamic analysis of stiffened plates Table 3. Frequencies of DRES simply supported square plate (rad/sec) Mode
F.S. (Consistent M)
F.S. (Lumped M)
Exact?
(l,l) (1.2) (1.3) (134) (IS)
402.66 1007.00 2016.80 4118.26 5753.46
397.41 952.36 1757.16 3362.60 3998.79
402.69 1006.73 2013.46 3422.89 5235.00
t v-21.
-III,
0
The nonlinear elastic-plastic response of the above model to the pressure load shown in Fig. 13 was calculated for both boundary condition cases. Hopefully, the experimental boundary conditions should lie somewhere between these two. Figure 14 shows a comparison of the experimental central displacement response with the two calculated results. The results for the simply supported case are reasonably close to the experiment. However, it appears that higher modes are contributing to the experimental results. To check this, the simply supported case was analysed again using a two mode approximation for each of the U, v and w displacements and using eight finite strips instead of four. The numerical integration was increased to 10-2-2 for the second model, while the normal 5-2-2 scheme was used for the first. Also the time step had to be reduced to 0.05 msec to avoid instability. The results of the two mode analysis are shown and compared with the experimental ones in Fig. 15. A modest improvement in the correlation of the results with experiment was achieved by the use of the second mode and more strips. No yielding was observed in any of the above calculations so the numerical results shown here are all strictly large deflection elastic.
r
I
08
1.6
2.4
3.2
I
4.0
T1me.t
48 (msec
56
6.4
7 2
I
8.0
)
Fig. 14. Response to DRES square plate (four strips, one mode). 3.4.1. Trial No. 314. The stiffened panel was suspended over a 96” x 180” opening, with the edges bolted to the foundations such that a relatively rigid connection was made between the plate and the foundation, simulating a clamped boundary. The panel was subjected to a blast resulting from the detonation of TNT explosives above the panel at the Height-of-Burst (HOB) site located at DRES [21]. The pressure transducers were used to establish the average pressure-time history of the air blast which is shown in Fig. 17. The square symbols on the plot show the data points used for the piecewise linear representation used in the program. At first 11 finite strips were used to model half of the stiffened panel, with symmetry boundary conditions used at the centre line. As shown in Fig. 18
13r
3.4. DRES stiffened panel Several 8’ x 15’ T-beam stiffened steel panels with l/4” plating and 3” x 6” T-beam stiffeners were tested at DRES. Figure 16 shows the details for these panels. Several different load cases were considered at DRES but only two are presented in the following along with their corresponding finite strip analyses. Note that the same time step of 5 psec was used for all the runs described in the following.
-III0
’
08
’
1.6
’
24
’
3.2
’
4.0
’
5.6
’
64
’
7.2
’
8.0
Tlme,t (msec 1 Fig. 15. Response of DRES square plate with extra modes and strips (simply supported in bending).
Table 4. Frequencies of DRES clamped square plate (rad/sec) Mode
Upper bound?
Lower bound?
F.S. (Consistent M)
F.S. (Lumped M)
(171) (192) (1,3)--(3,l) (1.3) + (3.1) (1.4) (175)- (591) (ITS)+ (531)
734.13 1497.28 2684.33 2697.06 4294.78 6301.81 6307.15
734.12 1497.25 2684.18 2696.90 4288.46 6290.70 6295.56
734.96 1497.19 2596.20
813.42 1576.21 2417.53
5224.32 7659.20
5777.45 7005.13
t ~231.
’
4.8
938
M. R. KHALILetal.
60
-1,; I“!‘--/;/0
5
IO
15
20 Time
25 ,t (msec
30
35
40
45
1
Fig. 17. Average pressure for DRES panel (HOB 314).
YOUNG’S
MODULUS,
POISSON’S
RATIO,
YIELD STRAIN
STRESS PER
TIME
STEP,
=
UNIT
X IO’psi
0.3
IN SIMPLE
HARDENING
MASS
E * 0.3 v
with six strips between stiffeners amounting to a total of 21 strips. In all three models, a one mode approximation was used for each of v and w displacements and a two mode approximation was used for the u displacement. Also a scheme of 5-2-2 Gauss points was used in all three models. Figure 19 shows a comparison of the panel deflection responses for the three models. The convergence with strip refinement is quite apparent. The results that follow were calculated with the 21 strip model. Figure 20 shows the deflection at the centre of the central panel (D12) and at the midpoint of the
TENSION,
MODULUS, VOLUME
ET=
, p =0.733
0.544X 0.178 x IO-’
IO’Psi
./e
X106PSi lb.sec%n4
At - 5 microseconds
Fig. 16. DRES stiffened panel geometry and data
only two strips were used for modelling the portion of the plate between adjacent stiffeners. The eigenvalue calculation results for this model are shown in Table 5 along with some computed by the VAST finite element analysis package [24]. The first few frequencies compare reasonably well. The nonlinear elastic-plastic response of this model to the dynamic loading was significantly lower than expected. Hence, two more models were analysed, one with four strips between stiffeners amounting to a total of 16 strips for the model and the other
x,u
/’ Fig. 18. Eleven strip model of DRES stiffened panel.
Table 5. DRES stiffened panel frequencies (Hz) Mode
VAST7
F.S. (Consistent M)
F.S. (Lumped M)
40.9 41.21 41.86
40.33 41.7 42.54
40.97 42.12 42.81
Mainly panel modes
Mode (1) Mode (2) Mode (3)
Mainly stiffener modes
Mode (4) Mode (5)
-
105.13 109.42
117.64 119.78
Mixed modes
Mode (6) Mode (7)
-
143.92 156.77
155.7 158.54
t [24].
939
Nonlinear dynamic analysis of stiffened plates
-2.0
2.0
6.0
14.0
10.0
12.0
22.0
-2.0
2.0
6.0
TIME(MS) (a)
Panel
10.0
Ii.0
Ii.0
25.0
TIME(MS)
- Dl2 ..-.-. ---
(a) I I strip model 16 strip model 21 stripmodel
Pane
I - Dl2
89 _( -2s I -.
-0 u t-i41 $9 z2 Ptf -2.0
2.0
8.0
10.0
14.0
18.0
!.O
TIME(MS) (b) -2.0
i.0
60
lb.0
Ii.0
lb.0
2h.o
TIME(MS) (b)
Stiffener
-
- D13
Fig. 20. DRES stiffened panel response to HOB 314: predicted and experimental displacements.
Dl3
Fig. 19. DRES stiffened panel response to HOB 314: effect of strip refinement.
adjacent stiffener (D13) along with the experimental results for comparison. The predicted and experimental results agree well up to 2 msec but beyond this
time the finite strip model appears to quickly become stiffer than the experimental panel. It should be mentioned here that the experimental deflections were not measured directly, but were calculated from double integration of the accelerations. There are some factors that may explain the difference in the results. The first factor is that only a few modes were used to approximate the displacements. The second factor is the boundary condition of the real plate, which may not be exactly clamped or fully restrained in-plane. Some slippage at the boundary bolts is also a possibility. 3.4.2. Minor scale. A stiffened panel identical to the one described above was mounted flush with the ground surface and subjected to an air blast load from a ground explosion. The panel was embedded in a reinforced concrete foundation to simulate clamped boundary conditions. The pressure-time history shown in Fig. 21 was obtained by averaging the measurements of several pressure transducers after synchronizing the rise times. The following
Stiffener
analysis was simplified
was found in the Trial No. 314 calculations that the improvement in response by the use of the 21 strip model instead of the 16 strip model was not significant. Therefore, only the 11 strip model and the 16 strip model were analysed in this example. The central displacements at D12 and D13 locations of both the 11 strip model and the 16 strip model are shown in Fig. 22. The figure shows some higher mode contributions in the response of the 16 strip model. Also at this load level the plate panel has deflected more than the stiffener as was the case in Trial No. 314. The results from the 16 strip mode1 are compared to the experimentally measured displacements in Figs 23 and 24. It is seen that the predicted and measured deflections are quite close for the first 14msec but thereafter the experimental results are higher than the finite strip ones, especially for the stiffener. Some of
__ 40 h
by ignor-
ing the transit time (about 3 msec) of the shock front over the panel (parallel to the stiffeners), and assuming the pressure loading to be spatially uniform over the panel surface. A nonlinear elastic-plastic dynamic analysis was performed using the above pressure-time history. It
IO
I J
0
IO
20
30 Time.t (msec)
40
50
Fig. 21. Average pressure for minor scale test
60
M. R. KHALILet
al.
and for the I-beam example. These results are rather preliminary and more detailed comparisons still need to be carried out. The 16 strip model used about 14 min of CPU time on an Amdahl 5850 for the 50msec calculation shown here. 4. CONCLUSIONS
-2
6
14
22 Time,
30
38
46
t ( msec)
Fig. 22. Response of DRES stiffened panel (minor scale) using 11 and 16 strips.
this discrepancy may be due to the approximations made in the pressure representation. However the fact that the numerical results are on the stiff side is
consistent with the previous results for the 314 test
The finite strip method has been applied to the large deflection, elastic-plastic dynamic response of plate-type structures to air-blast loads. The loads have been modelled as uniformly distributed time dependent pressures. The method is ideally suited to this type of load because the distribution can be well represented by a single mode. The method has been applied to a wide variety of configurations including rectangular beams, symmetric and asymmetric Ibeams, isotropic plates and stiffened plates. The results reported herein are very encouraging. It appears that reasonable engineering accuracy can
4CENTRE
PANEL
-56# -6
2
IO
18 26 Tlme,t(msec)
34
42
50
Fig. 23. Response of DRES stiffened panel (minor scale) using 16 strips-D12.
0L ;-24 s
-
-6
2
IO
I8
26
34
42
50
Tlme,t(msec)
Fig. 24. Response of DRES stiffened panel (minor scale) using 16 strips-Dl3.
Nonlinear dynamic analysis of stiffened plates be achieved by modelling involving only a few modes and a few finite strips. The number of variables and the computer times required are typically an order of magnitude less than those required by corresponding finite element codes. Thus considerable savings in cost are possible for certain problems with only a minimal loss in accuracy. The computer code is also much smaller and can be put on a micro-computer. The method seems ideally suited for preliminary design work. So far the main limitation seems to be in the application to the case of clamped boundary conditions for a structure which is highly efficient in bending such as a symmetric I-beam. The plastic deformation of such a structure is not modelled well by the clamped beam vibration mode. This is in contrast to the simply supported case which apparently is represented quite well by a half sine wave. The former problem needs further work, perhaps in finding a better mode to represent it. On the other hand, the stiffened panel with clamped boundaries seems to be adequately represented by the aforementioned beam modes. In particular, the displacement predictions for the DRES stiffened panel (Minor Scale test) as reported herein compare quite well with the test results. The combination of a diagonalized mass matrix and the finite difference time stepping method seems ideally suited to the present applications. The requirement of small time steps for stability is very compatible with the need for updating strains and stresses to follow yielding and plastic deformation. Further it is also needed for those cases with rapidly changing pressure-time history such as in the above Minor Scale test.
941
On the behaviour and design of stiffened plates in ultimate limit state. J Ship Res. 22, 238-244 (1978). 7. T. R. Stagliano and L. J. Mente, Large deflection, elastic-plastic dynamic structural response of beams and stiffened or unstiffened panels-a comparison of finite element, finite difference and modal solutions. In Nonlinear Finite Element Analysis and ADINA, Proc. ADINA Conference (Edited by K. J. Bathe et al.) (1979). 8. R. W. H. Wu and E. A. Witmer, Nonlinear transient responses of structures by the spatial finite element method. AIAA J. 11 (1973). 9. D. S. Mofflin, M. D. Olson and D. L. Anderson, Finite strip analysis of blast loaded plates. In Finite Element Methods for Nonlinear Problems (Edited by P. G. Bergan et al.), pp. 539-554. Springer, Berlin (1986). 10. R. Khalil, M. D. Olson and D. L. Anderson, Large deflection, elastic-plastic dynamic response of air-blast loaded plate structures by the finite strip method. Strut. Series Rept. No. 33, Dept. of Civil Engineering, UBC (1987). 11. S. Abayakoon, Large deflection, elastic-plastic analysis of plate structures by the finite strip method. Ph.D. thesis, Dept. of Civil Engineering, UBC, Vancouver (1987). 12. S. B. S. Abayakoon, M. D. Olson and D. L. Anderson, Large deflection elastic-plastic analysis of plate structures by the finite strip method. Int. J. Numer. Meth. Engng (in press). 13. Y. K. Cheung, Finite Strip Method in Structural Analysis. Pergamon Press, New York (1976). 14. S. R. Bonder and P. S. Symonds, Experimental and theoretical investigation of the plastic deformation of cantilever beams subjected to impulsive loading. 6. T. H. Soreide, T. Moan and N. T. Nordsve.
ASME,
J. appl. Mech. 29, 719-728 (1962).
15. A. Mendelson, Plasticity Theory and Application. MacMillan New York (1968). 16. L. M. Kachanov, Fundamentals of the Theory of Plasticity. [Translated From the Russian by M. Konyaeva, Mir Publishers. Moscow (1974).1 17. R. D. Krieg and S. W. Key, Transient shell response by numerical time integration. Int. J. Numer. Meth. Engng 7, 273-286 (1973).
Acknowledgement-This
research was supported by the Canadian Department of National Defence though a contract from the Defence Research Establishment Suffield.
REFERENCES
1. M. A. Chrisfield, The automatic nonlinear analysis of stiffened plates and shallow shells using finite elements. Proc. Inst. Cin. Engng, 69 (1980). 2. E. Ramm., Geometrisch-nichtlineare elastostatik und finite elemente. University of Stuttgart, Report 76-2 (1976). 3. P. A. Frieze, P. J. Dowling and R. E. Hobbs, Ultimate load behaviour of plates in compression. Proc. of Conf. on Steel Plated Structures, London. Crosby Lockwood Staples, London (1976). 4. J. E. Harding, P. E. Hobbs and B. G. Neal, The elasto-plastic analysis of imperfect square plates under in-plane loading. Proc. Inst. Cio. Engng, 63 (1977). 5. S. E. Webb and P. J. Dowling, Large deflection elastoplastic behaviour of discretely stiffened plates. Proc. Inst. Cia. Engng. 69 (1980).
18. R. D. Cook, Concepts and Application of Finite Element Analysis, 2nd edn. John Wiley, New York (1981). 19. J. W. Leech, P. Hsu and E. W. Mack, Stability of a finite difference method for solving matrix equations. AIAA Jnl3 (ll), 2172-2173 (1965). 20. B. R. Folz, Numerical simulation of the nonlinear transient response of slender beams. M.A.Sc. thesis, U.B.C., Vancouver, British Columbia (1986). 21. R. Houlston, D. Ritzel, J. E. Slater, G. Rude and R. T. Schmitke, Air blast experiments on square plates. DRES Suffield Memorandum No. 114153, Project No. 27C77, Defence Research Establishment Suffield, Ralston, Alberta (1986). 22. G. R. Cowper, E. Kosko, G. M. Lindberg and M. D. Olson, A high precision triangular plate-bending element. National Research Council of Canada, Aeronautical Report LR-5 14, Ottawa (1968). 23. L. De Vito, G. Fichera, A. Fusciardi and M. Schaerf, Sul Calcolo degli Autovalori della Piastra Quadrate Incastrata Lungo il Bordo. Rend. Accad. Nazion. Lincei, 8 Ser. 40, 725-733 (1966).
24. R. Houlston, J. E. Slater, N. Pegg and C. Desrochers, Structural response of ship panels subjected to air blast loading. 5th ADINA Conference, June 1985.