- Email: [email protected]
Dynamic response of thick plates and shells
Computers d Smrrura Vol. 33. No. I. PP. 63-72. Prinlcd in Great Britain.
1989
DYNAMIC
0
RESPONSE OF THICK AND SHELLS
0045.7949/89 $3.00 + 0.00 1989 Pergaman Press plc
PLATES
DEBAL BAGCHI,? NABIL F. GRACES and JOHN 8. KENNED@ TAdvanced Technology Department, Giffels Assoc., Southtieid, MI 48086-5025,
U.S.A.
@epartment of Civil Engineering, Lawrence Institute of Technology, Southheld, MI 4075, U.S.A. pbpartmcnt of Civil Engineering, University of Windsor, Windsor, Ontario, Canada N9B 3P4 (Received I4 July 1988) Abstract-A method for solving the nonlinear dynamic problem of moderately thick plates and shells is presented using an incremental finite element method. In the solution the nonlinear shear factor is derived and plasticfracture stress-strain relationships of progressively fracturing solids are used. Furthermore, an incremental force-displacement relationship is developed using the principle of virtual work. The minimum number of degrees of freedom for an element is employed in order to reduce the computer solution time. Numerical studies are presented to demonstrate the accuracy and efficiency of the method of solution.
1. INTRODUCTION
Solutions for the elastic bending of plates and shells accounting for transverse shear deformation have been derived by several inv~tigations [1,2]. In some cases the equilibrium equations were formulated with the transverse shear force included to arrive at a solution; in other cases an approximate displacement function, including shear rotation, was assumed. Thus the effect of transverse shear OR deflection, stresses and actual frequencies was determined. Considerable work has also been carried out on the inelastic response of plates and shells [l, 3-61. The objective of these works was to determine the effect of nonlinear material behavior on the deformation and stress distribution in a structure. Some of the most successful solutions were developed by means of the finite element method. Previous works using the finite element method had the following deficiencies: (1) constant strain within the element was assumed; (2) elements incorporating transverse shear formed from displacements do not satisfy equilibrium; (3) elements required postulated crack direction; (4) elements do not account for the variation of shear through the thickness of the plate or shell; (5) approximate form of nonlinear moment/ curvature relation ignoring transverse shear was used; (6) the isoparametric plate and shell elements used led to rather extremely large matrices which are unsuitable to analyze large structures. In this paper the above deficiencies are avoided by using an incremental finite element matrix formulation [7]. The element stiffness is developed considering higher-order transverse shear strain and nonlinear plastic-fracture stress-strain relations expressed in incremental form. The rotational degree of CAS 3311-E
freedom of the element includes total rotation due to pure bending and transverse shear; this facilitates the incorporation of the effect of rotating inertia. Subsidiary conditions are applied to reduce the number of degrees of freedom of the element as well as to improve the accuracy of the solution. A solution procedure is developed to examine the inelastic behavior of thick plates and shells of reinforced concrete or of laminated structures subjected to dynamic impact and seismic loadings encountered in nuclear and other industries.
2. THEORRTKAL FORMULATION
Consider a structure of volume V and density p in a state of dynamic eq~lib~~ with a known loading history. When such a structure is subjected to a small increment in an external force, d& (i = 1,2,3), on its surface St, changes in the internal stresses, deti, and in the displacements, dp,, will be induced. If the force increments are kept suffi~ently small, then the govequations relating strain-displacement, erning stress-strain, etc. can be linearized [8]. Furthermore, it is assumed that displacements are small and the resulting strains are of medium level so that the original geometry of the structure remains sensibly unaffected. Using the incremental relations given in [S], and based on the above assumptions, the principle of virtual work may be written for the present problem in a time domain from t to t + dr as
-
s
d&6(dp,) dS $1
1
dt ~0,
(1)
64
D. BAOCHIet al.
where dc, is the incremental strain tensor, d& is the incremental velocity; and dt is the time increment. The above principle of virtual work in the discrete finite element method is applied to each element of the discretized structure with the element nodal degrees of freedom treated as generalized coordinates to obtain the stiffness matrix [A?‘],mass matrix [M’] and nodal force vector (dF) of the element. This principle is then applied to the whole structure in which the elements are connected at the nodes. To include the energy dissipation due to damping [9] used in the equation of motion in structural dynamics, it is assumed that the dissipation force is proportional to the instantaneous mass matrix by a multiplier, 07, and to the instantaneous stiffness matrix by a multiplier, 8. Thus, applying eqn (1) to the whole structure and integrating by parts yields
+ [Kl {dp ) - (dF))S (dp ) dr = 0,
can be expanded in a Taylor series in terms of xj [8] as
x3 + . **(higher terms),
(4)
where dz, is the incremental midplane displa~ment. Neglecting the higher terms yields dpi=dE, +x,dp;T
(4a)
where dp; is the derivative of dpl with respect to x3; it should be noted that dp; and dz, are functions of xI , x2 only, and $, = 0 for small displacement theory. Now confining our attention to bending only, we can write strains for a thick plate element as
dtz =
(2)
ad.4 x3 dx; 2
where [Ml, [Kl, (dF) are the assembled mass, stiffness and nodal force matrices of the structure, respectively, and {dp} is the generalized displacement matrix of the nodes. Since eqn (2) is true for any arbitrary value of variation of displacement b(dp}, the equation of motion of the structure in incremental form is obtained by equating the expression in the first bracket to zero. Or, for any time step r,
ad&+adpj
d&12”
de
x3
(
--ax,
ax2
a +3
--+d/.i;; ax,
r3-
and
adlt3
(5)
dez3= -+d&. a x2
2.1. Derivation of mass, stiffness, and load matrices for a thick plate element at any time step Consider the thick plate element in Fig. 1, with coordinates x, , x2 and x,. The resulting incremental displacement dp,, being a function of xi (i = I, 2,3),
Consider a rectangular plate element of four nodes with three degrees of freedom at each node. Its bending displa~ments can be assumed in ~lynomi~ series as d@, = -xj (aI + olzxl + @,x2) ab -X3-p
dlcz=
-x3&
+
9x1
+
w2)
a dii, --x3-Q dp3 = a1 + agxl + a9x2 + ~~~~~x2 + a,,x: f ai2x: + at3x: + a,,xix2 + aIsxIxi + a16xi 2
2
+ a,,xix2 + a18xIxi + alpxIx2,
Fig. I. Deformed shape of a thick plate element.
(6)
where a, to a,9 are constants and h is the plate element thickness. The fourth term of the power series, in-
Dynamic response of thick plates and shells
eluded in dpr and dgr, reflects the influence of the parabolic shear distribution. At any stage of loading, the 19 constants given in eqn (6) are evaluated as follows:
obtained can be expressed in matrix form as [H”](a)
(3) Combining (1) Substituting the 12 unknown nodal degrees of freedom in eqn (6) results in the matrix equation {@f =
WI(~),
65
= [E’] dq = (E).
(10)
eqns (7) and (10) yields
(a} = [If-’
‘*
11E
;
(11)
(7) where
where (#‘I is the column matrix of nodal degrees of freedom; {a} is the column matrix of unknown constants; and [H’] = 12 x 19 rectangular matrix. (2) Considering the incremental stress-strain relations, incremental moments and shear forces can be written as h/2 dMii =
- x3 cjfil dck,dx,; I -h/2
$j=
k,l=
1,2
~d~~=[~l~~*~+[~~ldq,
1,3
where [A] relates the displacements to the nodal degrees of freedom and [AA] relates the displacements to the distributed loading. Substituting eqn (12) into (5), the strain~ispla~ent relationship can now be expressed as
h/Z
dQ, =
The 19 constants, a),. . . , a,9, can now be evaluated from the 19 equations. After evaluating the 19 constants, eqn (6) can be written in a matrix form in terms of the 12 nodal degrees of freedom (S *) and the distributed loading dq, as
I
Cm ki h ; -h/Z
02)
i=l,2 (de} = k, I= 1,3,
where I?,, is the incremental plastic-fracture material moduli [IO] defined later on. The incremental equations of equilibrium for the plate are given by
a d&f,* a d&it2 --ax+dQrr=O; ax, I 8 dMr, ax2 a
a dM,, --dQ,,=O; + ax,
dQx + a dQ23 ax+dq ax, z
PI 16”)f 101dq,
(13)
(8)
=0,
(9)
where {de>‘= {de,, d&22de,r d+ dezJj, [B] is a matrix relating the strains to the nodal degree of freedom and [B] is a matrix relating the strains to the distributed loading. The second terms in eqns (12) and (13) are used to calculate the additional effect on deflections, strains and stresses within the element, but will not be used in the derivation of the mass, load and stiffness matrices for the plate element. The constitutive relations [lO-141 for the plate material at any time step in a linear incremental form can be represented by (14)
where dq is an incremental distributed load. It should be noted that eqn (8) shows the coupling between the direct stress and shear strain as a consequence of the assumed inelastic stress-strain relation as well as the shift of the neutral axis. Substituting eqn (8) into (9) and evaluating strains from eqns (5) and (6) yields six equations relating the constants a,, . . . , a,9, from the first two equations of (9) after equating the coefficients of x1 and x2 to zero. The third equation of (9) will relate the constants with the distributed load. It should be noted that the effect of variation in the moduli through the thickness of plate is ignored by considering average moduli. Furthe~ore, the inertia terms ignored in eqn (9) are applied as a correction term based on the solution of the previous time step. Thus, the seven equations
where C,,, is the incremental plastic-fracture modulus [lo]. While eqn (14) is available for plain concrete, it requires modifi~tion for reinforced concrete structures. For the latter it is assumed that the composite incremental modulus can be written as Pa)
where C&, is the materiai modulus of reinforcement and V, is the fraction of concrete in unit volume. The material modulus due to reinforcement is evaluated by considering material modulus of the ~info~ment about its axis times its volume fraction and then transformed to the concrete axis. It is assumed that the presence of reinforcement will not
D. BAOCHIef al.
66
affect the plastic-fracture modulus of plain concrete. Equation (14) can be written in a matrix form as {da} = [C]{de}.
tuting these strains into eqn (16) yields
(14b)
Applying eqn (1) to the plate element, using eqns (12x14) the mass, stiffness and load matrices for a moderately thick plate element become and pATA dV, Y
[M’] =
[K’] =
s”
{dF’} =
B=CB dV,
dqA ‘dS.
(15)
To derive the nonlinear shear factor, approximate forms of the following equilibrium equations are used: a dM,, --de,,=0 axi
2.2. Derivation of mass, stiffness, and load matrices for a thick shell at any time step
a dMzz --dQ,,=O,
ax2
where h/2 dMii 1
- x, Ciiiide, dx,
;
s -h/2 and
s hi2
dQo =
where a, 6 are the length and width of the plate element; the factors 0, and G2 are the nonlinear shear factors in the directions of the curvatures rc, and K~. It can be observed that these factors can also be defined as the ratio of bending rigidity to the shear rigidity based on tangential material moduli.
c,,, dc, dx, , i = 1,2;
-h/2
and strains dc,, dc, are calculated using eqns (5) and (6) in terms of the unknown constants, {a}. Substi-
Consider the shell element in Fig. 2(a) with curvilinear axes x,, x2 and x3. Figure 2(b) shows an axisymmetric shell in which do is the angle in the meridional plane between the apex and the reference point in the shell at which the radii of curvatures R, and ROin the meridional direction (x, ) and circumferential direction (x2) are considered. Expressions for the displacement, dpi, for the thick shell are given by eqn (4) where d:, and d,ui are functions of x, , x2 only and d@; = 0 for small displacement theory. The straindisplacement relations for a thick shell element including the effect of transverse shear deformation can be readily found elsewhere [2,8]. The incremental displacement functions for a general axisymmetric shell with nonaxisymmetric load-
b) Fig, 2. (a) Deformed shape of a thick shell element; (b) an axisymmetric shell element with four d.o.f. at each node.
Dynamic response of thick plates and shells
ing can be represented as dp, = dp,,, cos nfl, +dM,,cos4
-rdQ,,=O
dp2 = dpg sin n0, dp, = dp,, cos no,
(1%
$-trdQ13)+n
dQ23
I
where
-
G,n = $,. -x3
(af,+aix,)
where hi* C,,kld6k,dX3,i,i=l,2
dNii =
k,l=l,3.
(21a)
s -h/2
G2n= 4,
-xl
(ai+a:x,)
4x:
( )
x 1--
-5
3h2
x(-n
r
dp3,+d0,,“sin4);
dp3,=a~+a~x,+a~x: + ai xi + a:‘~:;
d:,” = a1:+ aA0x,; d:,” = af’ + ai* x,;
dpi, is the incremental displacement for Fourier component, n; 4, is the incremental midplane displacement for Fourier component, n; ai are the constants for Fourier component, n and i = 1,. . . ,13; h is the thickness of the shell; and R, and r are defined in Fig. 2(b). At any stage of loading, the 13 constants given by eqn (19) for each Fourier component, n, are evaluated as follows. (1) Substituting the eight nodal degrees of freedom, defined in Fig. 2(b), in eqns (19), will result in
Evaluating the shell strain-displacement relations [2,8] in terms of variable x, and a and substituting them into eqn (21), yields four equations from the fist two equations of (21) by equating the coefficients of the variable x, for each equation to zero. The third equation of (21) will relate the constants a to the distributed load. For the purpose of calculating the constants, moments dM2*, dM,, are assumed to vary linearly in the x, direction. As in thick plate, the effect of variation in the moduli through the thickness of the shell is ignored by considering average moduli. The resulting five equations can be expressed in matrix form similar to eqn (10). (3) Combining eqns (20) and the above five equations results in an equation similar to (11). Solving this will yield the 13 constants in terms of the eight nodal degrees of freedom and the distributed load. After evaluating these constants and expressing the distributed load in terms of a cosine and sine series, eqn (19) can be written for each Fourier term in matrix form as IdKj =
[Al{@ij+[AAl{dqnIt
where (Si} are the nodal displacements and {dq,,} is the distributed load for the nth Fourier component. Using eqn (22), the strain-displacement relations for each Fourier component take the form:
W,} = PI{%} + Pl{dqnIv which is similar in form to eqn (7) except that [H’] is an 8 x 13 rectangular matrix. (2) The incremental stress-strain relation and moment and shear relations given in eqns (8) and (14) are now substituted in the following equations of equilibrium for a given Fourier component n for axisymmetric shell [2]:
&trdMll)+n
dM,2
I
-dM,,cos4
-rdQ,,=O
(22)
(23)
where {de,}r= {dCun, dC22nrdc,,, dt,3n, dcr,,}. The remaining procedure is similar to that for the thick plate element. It should be mentioned that under a nonaxisymmetric loading, the plastic moduli [C] changes with the angle 0 and as such, decoupling of eqn (3) for each harmonic cannot be done readily. Since integration of the element from 0 to 2x at each step is extremely time-consuming, three approaches are suggested, namely: the average modulus approach; the average effective stress approach; and the approach in which the stress-strain relations are expanded in a Fourier series. However, due to decoupling difficulties, the first approach was applied
68
D. BAGCHIet al.
herein. Thus an average modulus independent angle, 8, can be calculated as
of the
c =; i c,, ’
k=l
(24)
where C, is the actual modulus at an angle, 8,; and t is the number of points along the circumference where the actual modulus is calculated. As in the case of the thick plate element the incremental stiffness, mass and load matrices for a moderately thick shell element can now be obtained from eqn (15) for each Fourier component n. In deriving the nonlinear shear factors for the shell element, the following approximate forms of the shell equilibrium equations for moments are used: &dY,,)-d&cos+ I
-rdQ,,=O n dMz2 - r de,, = 0,
(25)
where dA4,i, dQ13are defined before in eqn (8). Evaluating dM,, dQa from the stress-strain relations and strain-displacement relations [2,8] in terms of variable x, and a and substituting in eqns (25), yields
where @,, and 9,, are the nonlinear the x, and x2 directions. 2.3. Solution procedure
The incremental global (structure) mass, stiffness and force matrices, [Ml, [K] and [F] are determined by assembly of the element mass, stiffness and force matrices [eqn (15)] after a suitable coordinate transformation. The dynamic equation of motion at a time step t [eqn (3)] is solved numerically by means of a backward finite difference scheme and a Gaussian integration procedure [ 1, 151. For example, the numerical solution of eqn (3) at a particular time step t, is found based on solutions at previous time steps (t - 2) and (t - 1). Thus, for a = 0, one can write for a small time interval At:
&JW’sJ =
{AF,
+ 1 +
-&It-2 dM,, cos 4
and
where oi, = ai +2azl+ 3ai1* + 4ai313; and I is the length of the element in the x, direction. For the case of a cylindrical shell, R, = CC and rj~= 90”, eqn (25) becomes h/2 6c x: Cm, dx3 ,
shear factors in
(& +1)Wi4
&fl{b-J (At)2
&-,I.
(28)
Depending on the nature of loading, it may be necessary to modify eqn (28) for the first two time steps. At the end of every time step solution, it is necessary to modify the mass and stiffness matrices to reflect changes in the stress-strain relation and in the shear rigidity. For a particular time t, the total strains and stresses are determined from the summation of incremental strains and stresses calculated previously. In general, the solution has been found to be stable (for example, up to the time before yielding of the steel in a reinforced concrete structure, say). However, experience has shown that instability in the solution can be circumvented by considering smaller time steps. A computer flow chart for the incremental solution is shown in Fig. 3.
3. DISCUSSION
OF RESULTS
The convergency of the solution was examined. Figure 4 shows the rapid rate of convergence of the result for the first natural frequency with increase in the number of elements for a cylindrical shell simply supported at both ends. Another comparison is made in Fig. 5 for the radial deflection response for a simply supported axisymmetric shell under suddenly applied ring load; good correspondence is shown with
Dynamic
READ:
response of thick plates and shells
STRUCTURAL GEOMETRY, MATERIAL PROPERTIES, BOUNDARY CONDITIONS 6 LOAD FUNCTION
I
4
FORM MASS h
ELASTIC STIFFNESS MATRICES OF EACH ELEMENT @ t=o
I FORM
GLOBAL MASS [Ml AND STIFFNESS MATRIX [Kl
<
1
FORM LOAD MATRIX
(Fj
SOLVE EQUATION (28)
FROM INPUT
LOAD FUNCTION1
TO OBTAIN DISPLACEMENT, Apt
CALCULATE FOR EACH ELEMENT: (i) STRAINSAND STRESSES; (ii) NONLINEAR SHEAR FACTOR USING CURRENTMATERIAL
MODULI
. REVISE MATERIALMODULI BASED ON REVISED TOTAL STRAINS AND STRESSESIN THE ELENBNT
ALL TIME STEPS
FORM MASS 6 STIFFNESS - MATRICES OF EACH ELEMENT BASED ON REVISED MODULI
PRINT RESULTS
Fig. 3. Flow chart for the incremental
results from an exact solution [la], and those from a finite element solution [ 171. 3.1. Results from elastic solution Deflection results at the center of a clamped square thin plate (h/a = 0.01) were compared with results from the classical solution (3) which neglects shear effects; the comparison shows that the results from the present solution are about 10% higher than those given by [ 181;this indicates that shear has an influence on the deflections even in thin plates. Deflection results for a simply supported relatively thick plate (h/a = 0.2) are only 2% higher than those derived using an eight-node thick plate element [l]. Further-
solution.
more, results for a clamped square plate show that transverse shear and rotatory inertia have no significant effect on the translational mass (M). However, they do have a significant effect on the natural frequency as the (h/a) ratio increases. [A decrease of 7% is noted for h/a = 0.1 (shear factor Q = 0.014) and a decrease of 21% for h/a = 0.3 (shear factor @ = 0.128).] The effect of the shear factor 4~ on the radial deflection of a short (I = 12.7’) and thick (h = 5’) cylindrical shell (radius r = a = 12’) was examined. The shell was fixed at the base and subjected to an axisymmetric load and then to a nonaxisymmetric edge load P at the upper free edge of the shell. For
D. BAGCHIet al.
IQ
I
x 14oDD-
r
5 12000-
3 10000z
SOOO-
t, f
6000-
(Ref. 17)
2 4000I Ia- 2DoD04 NO. OF ELEMENTS
USED
Fig. 4. Convergence of solution with increase in number of elements (mode, m = 0).
the axisymmetric load the present solution shows an increase of 10% in the radial deflection over the values derived from the classical solution [19] which neglects the shear effect; the increase in deflection is of course due to the shear effect; however, this percentage does not vary significantly with changes in @. For the nonsymmetric load (P cos 0) the shear factor @was found to be very significant, as expected; a 40% increase in deflection was observed when the shear factor increased from @ = 0.00084 to 0 = 0.337.
compressive strength of concrete = 4000 psi with longitudinal, circumferential and diagonal steel reinforcements with ratios of 1.5, 1.5 and l.O%, respectively. It is assumed that the shell is fixed at its two ends but axially free to move at its top end, with a length of 60 ft, wall thickness of 5 ft, and a radius of 60 ft. It is observed from Fig. 6 that when the longitudinal and circumferential steels yield and the concrete cracks near the mid-height of the shell (point U on curve OW, Fig. 6) the resulting loss of stiffness effects a rapid increase in the deflection (point V on curve OW, Fig. 6); continued loading triggers complete yielding of the longitudinal steel. This causes the transfer of load to the diagonal reinforcing steel which can resist the load in both the circumferential and longitudinal directions; also, the rate of increase in the radial deflections with time decreases as shown by curve VW in Fig. 6. It is evident that by increasing
ii SUDDENLY APPLIED RING(LINE)LOAl
3.2. Results from nonlinear-plastic solution Results for the nonlinear dynamic behavior of a reinforced concrete thick shell subjected to an axisymmetric impulsive load are shown in Fig. 6. Such loading may be due to an impact of seismic event. The translational mass matrix has been magnified by a factor of 10 to account for the effect of the surrounding water mass. The nonlinear radial deflection-time response is based on thick-she11 theory and on the plastic-fracture behavior of concrete as described in[lO]. It is assumed that the
TIME
t(* 10%~~)
Fig. 5. Radial deflection-time response of shell under suddenly applied ring load.
Diagonal Steal at bottom edge Started Yielding
Yielding of Circumf. 6 La&: Steel in middle of Shell, and Crac mg of Concrete
Fig. 6. Inelastic radial deflection-time
response at location A of reinforced concrete shell under axisymmetric ramp load.
Dynamic response of thick plates and shells
TIME
71
increase in the strain cl2 with the shear stress u12 remaining sensibly constant. Subsequent severe cracking of the concrete leads to a transfer of load to the diagonal reinforcing steel and its eventual yielding. After this, the increase of stress with strain is essentially due to material hardening effect. The same cylindrical shell was analyzed for a nonaxisymmetric load (P cos 0) as shown in the inset of Fig. 7; the shell was assumed to be fixed at one end and free at the other end. The results in Fig. 7 are based on Bazant’s plastic-fracture constitutive law [lo] and use the average modulus concept explained earlier. The radial deflection results for point B near the free edge show significant increase after the initiation of yielding in the diagonal reinforcing, accompanied by extensive cracking in the concrete. The shear stress-strain (u~~/c,~) relation for point A on the shell (see inset in Fig. 7) in Fig. 8 reveals an initial drop in stiffness at the initiation of cracking in the concrete followed by a subsequent
in seconds
Fig. 7. Inelastic radial ddlection-time response at location B of reinforced concrete shell under nonaxisymmetric first cosine ramp load.
the percentage of diagonal reinforcing to its limit value, the radial deflections at the center of the shell will be significantly reduced. Results for the stress-strain relations for the concrete were also determined but are not shown herein for brevity. The corresponding circumferential stress-strain (a2r/czz) relation was significantly nonlinear, exhibiting degrading stiffness until the final load was reached at a stress a,, = 5640 psi for a strain cu = 0.0075. The axial stress-strain (ur, /L,,) relation was essentially linear, perhaps due to the triaxial nature of compression in concrete. The results also show that transverse shear stress-strain (uu/+) relation is linear; however, upon yielding of the circumferential reinforcing and initial cracking, continued loading will incur plastic flow in the material, i.e. significant
Cracking
0
.OOl
I
222
444 LOAD
666
,
008
I
1110
(P) in psi
Fig. 9. Variation of shear factor with increase in load.
of Concrete
.002 SHEAR
Fig. 8. Shear stress-strain
0
.003 STRAIN
.004
C,,
relation at location A at the inner surface of shell in Fig. 6 at angle B= 90”.
D. BAGCHI et al.
72
increase in stiffness due to the presence of the longitudinal and circumferential steels. A subsequent drop in stiffness is observed after yielding of the latter steels; again, due to the presence of the diagonal reinforcing steel there is a subsequent increase in
(4) Composite material such as reinforced concrete exhibits severe nonlinearity in its shear factor @ as a result of tension cracks; this state of affairs would lead to structural failure at a much lower load than that predicted by elastic analysis.
stiffness followed by a levelling off in stiffness when ail steels had yielded. When compared with the response of the structure in Fig. 6, it is observed from
REFERENCES
the results in Figs 7 and 8 that the load-carrying capacity of the shell shown in Fig. 7 is much less. The variation of the nonlinear shear factor @, defined earlier by eqn (27), with load is presented in Fig. 9 for the shells and loadings in Figs 6 and 7; for the latter shell, two values of 6 are considered, namely 6 = 0” and 90”. It is observed that the shear factor remains essentially unaffected by increase in the load; however, when the concrete cracks and the diagonal reinforcing steel had yielded subsequently, the shear rigidity decreased rapidly, causing a sharp increase in the nonlinear shear factor @. Furthermore, the difference in the values of Q, for location C in Fig. 7 for inclination 6 = 0 and 8 = 90” does not appear to be significant.
0. C. Zienkiewicz, The Finite Element Methods in Engineering Science. McGraw-Hill, London (1977). H. Kraus, Thin Elastic Shells. John Wiley, New York
(1967). D. Bushnell, Buckling of elasto-plastic shells of revolu-
tion with discrete elasto-plastic-dng stiffeners. Znt. J. Sot%& Sirwt.
ticity.
4. CONCLUSIONS An incremental for~isplacement fo~ula~on for the solution of the nonlinear dynamic problem of thick plates and shells under arbitrary loading has been introduced using a discretized finite element energy approach. It is shown that: (1) The plate and shell elements developed can predict fairly accurately the structural response of relatively thick plate and shell structures; the deformation of these efficient elements is represented by a minimum number of degrees of freedom. (2) Effects of rotatory inertia and transverse shear strain are accurately represented; furthermore, the shear factor is estimated from the tangent moduli at every increment of loading. (3) The three-dimensional plastic-fracture characteristics of heterogeneous and anisotropic materials, such as reinforced concrete and laminated composites, can be readily treated, including cycles of loading-unloading of the structure due to cracks. Therefore, the method does not rely on a postulated crack direction nor on the use of an approximate moment-curvature relationship.
12, 51-66 11976).
4. M. P. Bieneck add J. R. l%naio, Elastic-plastic behavior of plates and shells. Tech. Report, Contract No. DNAOOl-76-C-0125, Defense Nuclear Agency, Washington, DC (1976). 5. R. F. Jones and J. E. Roderick, Incremental solution procedure for nonlinear structural dynamic problem. J. Pressure Vessel Technol. 97. 95-100 (1975). 6. S. Klein, The nonlin~r dynamic analysis’of shells of revolution with asymmetric properties by the finite element method. Trans. ASME 169, 163-170 (1975). 7. D. K. Bagchi, Inelastic bending of beam including transverse shear. Int. J. Nurner. Meth. Engng 14, 1323-1333 (1979). 8. K. Washizu, Variation Methods in Elasticity and Plas-
Pergamon Press, Oxford (1968).
9. H. F. Rubenstein and W. C. Hurty, Dynamics of Structures. Prentice-Hall, Englewood Cliffs, NJ (1964). 10.Z. P. Bazant and Sung-Sin Kim, Plastic-fmctu~ng theory for concrete. Engng Mech. Div., ASCE 105, (EM3), 407-428 (1979). 11. R. Hill, The Mathematical Theory of Plasticity. Clarendon Press, Oxford (1950). 12. B. Budiansky and R. J. O’Connell, Elastic moduli of a cracked solid. Int. J. Solids Struct. 12, 81-97 (1976). 13. J. W. Dougill and M. A. M. Rid, Further consideration of progressively fracturing solids. Errgag Mech. Div., ASCE 106(EM5), 1021-1038 (1980). 14. Zvi Hashin, D. K. Bagchi and B. Rosen, Nonlinear behaviour of fiber composite laminates. NASA Contract Report prepared under Contract NASI-11284 (1973). 15. R. S. Burington, Handbook ofMathematical Tables and Formulas. 5th Edn. McGraw-Hill. New York (1973). 16. H. Reismann and J. Padlog, Forced, axi-symmetric motions of cylindrical shells. J. Franklin Inst. 284(S), 308-319 (1967).
17. S. Ghosh and E. Wilson, Dynamic stress analysis of axisymmetric structures under arbitrary loading. Rep. No. EERC 69-10, College of Engineering, University of California, Berkeley (1969). 18. S. Timoshenko and S. Woinowesky-Krieger, Theory of Plates and Shells, 2nd Edn. McGraw-Hill, New York (1961). 19. R. J. Roark, Formulae for Stress and Strain. McGraw-Hill, New York (1965).
1989
DYNAMIC
0
RESPONSE OF THICK AND SHELLS
0045.7949/89 $3.00 + 0.00 1989 Pergaman Press plc
PLATES
DEBAL BAGCHI,? NABIL F. GRACES and JOHN 8. KENNED@ TAdvanced Technology Department, Giffels Assoc., Southtieid, MI 48086-5025,
U.S.A.
@epartment of Civil Engineering, Lawrence Institute of Technology, Southheld, MI 4075, U.S.A. pbpartmcnt of Civil Engineering, University of Windsor, Windsor, Ontario, Canada N9B 3P4 (Received I4 July 1988) Abstract-A method for solving the nonlinear dynamic problem of moderately thick plates and shells is presented using an incremental finite element method. In the solution the nonlinear shear factor is derived and plasticfracture stress-strain relationships of progressively fracturing solids are used. Furthermore, an incremental force-displacement relationship is developed using the principle of virtual work. The minimum number of degrees of freedom for an element is employed in order to reduce the computer solution time. Numerical studies are presented to demonstrate the accuracy and efficiency of the method of solution.
1. INTRODUCTION
Solutions for the elastic bending of plates and shells accounting for transverse shear deformation have been derived by several inv~tigations [1,2]. In some cases the equilibrium equations were formulated with the transverse shear force included to arrive at a solution; in other cases an approximate displacement function, including shear rotation, was assumed. Thus the effect of transverse shear OR deflection, stresses and actual frequencies was determined. Considerable work has also been carried out on the inelastic response of plates and shells [l, 3-61. The objective of these works was to determine the effect of nonlinear material behavior on the deformation and stress distribution in a structure. Some of the most successful solutions were developed by means of the finite element method. Previous works using the finite element method had the following deficiencies: (1) constant strain within the element was assumed; (2) elements incorporating transverse shear formed from displacements do not satisfy equilibrium; (3) elements required postulated crack direction; (4) elements do not account for the variation of shear through the thickness of the plate or shell; (5) approximate form of nonlinear moment/ curvature relation ignoring transverse shear was used; (6) the isoparametric plate and shell elements used led to rather extremely large matrices which are unsuitable to analyze large structures. In this paper the above deficiencies are avoided by using an incremental finite element matrix formulation [7]. The element stiffness is developed considering higher-order transverse shear strain and nonlinear plastic-fracture stress-strain relations expressed in incremental form. The rotational degree of CAS 3311-E
freedom of the element includes total rotation due to pure bending and transverse shear; this facilitates the incorporation of the effect of rotating inertia. Subsidiary conditions are applied to reduce the number of degrees of freedom of the element as well as to improve the accuracy of the solution. A solution procedure is developed to examine the inelastic behavior of thick plates and shells of reinforced concrete or of laminated structures subjected to dynamic impact and seismic loadings encountered in nuclear and other industries.
2. THEORRTKAL FORMULATION
Consider a structure of volume V and density p in a state of dynamic eq~lib~~ with a known loading history. When such a structure is subjected to a small increment in an external force, d& (i = 1,2,3), on its surface St, changes in the internal stresses, deti, and in the displacements, dp,, will be induced. If the force increments are kept suffi~ently small, then the govequations relating strain-displacement, erning stress-strain, etc. can be linearized [8]. Furthermore, it is assumed that displacements are small and the resulting strains are of medium level so that the original geometry of the structure remains sensibly unaffected. Using the incremental relations given in [S], and based on the above assumptions, the principle of virtual work may be written for the present problem in a time domain from t to t + dr as
-
s
d&6(dp,) dS $1
1
dt ~0,
(1)
64
D. BAOCHIet al.
where dc, is the incremental strain tensor, d& is the incremental velocity; and dt is the time increment. The above principle of virtual work in the discrete finite element method is applied to each element of the discretized structure with the element nodal degrees of freedom treated as generalized coordinates to obtain the stiffness matrix [A?‘],mass matrix [M’] and nodal force vector (dF) of the element. This principle is then applied to the whole structure in which the elements are connected at the nodes. To include the energy dissipation due to damping [9] used in the equation of motion in structural dynamics, it is assumed that the dissipation force is proportional to the instantaneous mass matrix by a multiplier, 07, and to the instantaneous stiffness matrix by a multiplier, 8. Thus, applying eqn (1) to the whole structure and integrating by parts yields
+ [Kl {dp ) - (dF))S (dp ) dr = 0,
can be expanded in a Taylor series in terms of xj [8] as
x3 + . **(higher terms),
(4)
where dz, is the incremental midplane displa~ment. Neglecting the higher terms yields dpi=dE, +x,dp;T
(4a)
where dp; is the derivative of dpl with respect to x3; it should be noted that dp; and dz, are functions of xI , x2 only, and $, = 0 for small displacement theory. Now confining our attention to bending only, we can write strains for a thick plate element as
dtz =
(2)
ad.4 x3 dx; 2
where [Ml, [Kl, (dF) are the assembled mass, stiffness and nodal force matrices of the structure, respectively, and {dp} is the generalized displacement matrix of the nodes. Since eqn (2) is true for any arbitrary value of variation of displacement b(dp}, the equation of motion of the structure in incremental form is obtained by equating the expression in the first bracket to zero. Or, for any time step r,
ad&+adpj
d&12”
de
x3
(
--ax,
ax2
a +3
--+d/.i;; ax,
r3-
and
adlt3
(5)
dez3= -+d&. a x2
2.1. Derivation of mass, stiffness, and load matrices for a thick plate element at any time step Consider the thick plate element in Fig. 1, with coordinates x, , x2 and x,. The resulting incremental displacement dp,, being a function of xi (i = I, 2,3),
Consider a rectangular plate element of four nodes with three degrees of freedom at each node. Its bending displa~ments can be assumed in ~lynomi~ series as d@, = -xj (aI + olzxl + @,x2) ab -X3-p
dlcz=
-x3&
+
9x1
+
w2)
a dii, --x3-Q dp3 = a1 + agxl + a9x2 + ~~~~~x2 + a,,x: f ai2x: + at3x: + a,,xix2 + aIsxIxi + a16xi 2
2
+ a,,xix2 + a18xIxi + alpxIx2,
Fig. I. Deformed shape of a thick plate element.
(6)
where a, to a,9 are constants and h is the plate element thickness. The fourth term of the power series, in-
Dynamic response of thick plates and shells
eluded in dpr and dgr, reflects the influence of the parabolic shear distribution. At any stage of loading, the 19 constants given in eqn (6) are evaluated as follows:
obtained can be expressed in matrix form as [H”](a)
(3) Combining (1) Substituting the 12 unknown nodal degrees of freedom in eqn (6) results in the matrix equation {@f =
WI(~),
65
= [E’] dq = (E).
(10)
eqns (7) and (10) yields
(a} = [If-’
‘*
11E
;
(11)
(7) where
where (#‘I is the column matrix of nodal degrees of freedom; {a} is the column matrix of unknown constants; and [H’] = 12 x 19 rectangular matrix. (2) Considering the incremental stress-strain relations, incremental moments and shear forces can be written as h/2 dMii =
- x3 cjfil dck,dx,; I -h/2
$j=
k,l=
1,2
~d~~=[~l~~*~+[~~ldq,
1,3
where [A] relates the displacements to the nodal degrees of freedom and [AA] relates the displacements to the distributed loading. Substituting eqn (12) into (5), the strain~ispla~ent relationship can now be expressed as
h/Z
dQ, =
The 19 constants, a),. . . , a,9, can now be evaluated from the 19 equations. After evaluating the 19 constants, eqn (6) can be written in a matrix form in terms of the 12 nodal degrees of freedom (S *) and the distributed loading dq, as
I
Cm ki h ; -h/Z
02)
i=l,2 (de} = k, I= 1,3,
where I?,, is the incremental plastic-fracture material moduli [IO] defined later on. The incremental equations of equilibrium for the plate are given by
a d&f,* a d&it2 --ax+dQrr=O; ax, I 8 dMr, ax2 a
a dM,, --dQ,,=O; + ax,
dQx + a dQ23 ax+dq ax, z
PI 16”)f 101dq,
(13)
(8)
=0,
(9)
where {de>‘= {de,, d&22de,r d+ dezJj, [B] is a matrix relating the strains to the nodal degree of freedom and [B] is a matrix relating the strains to the distributed loading. The second terms in eqns (12) and (13) are used to calculate the additional effect on deflections, strains and stresses within the element, but will not be used in the derivation of the mass, load and stiffness matrices for the plate element. The constitutive relations [lO-141 for the plate material at any time step in a linear incremental form can be represented by (14)
where dq is an incremental distributed load. It should be noted that eqn (8) shows the coupling between the direct stress and shear strain as a consequence of the assumed inelastic stress-strain relation as well as the shift of the neutral axis. Substituting eqn (8) into (9) and evaluating strains from eqns (5) and (6) yields six equations relating the constants a,, . . . , a,9, from the first two equations of (9) after equating the coefficients of x1 and x2 to zero. The third equation of (9) will relate the constants with the distributed load. It should be noted that the effect of variation in the moduli through the thickness of plate is ignored by considering average moduli. Furthe~ore, the inertia terms ignored in eqn (9) are applied as a correction term based on the solution of the previous time step. Thus, the seven equations
where C,,, is the incremental plastic-fracture modulus [lo]. While eqn (14) is available for plain concrete, it requires modifi~tion for reinforced concrete structures. For the latter it is assumed that the composite incremental modulus can be written as Pa)
where C&, is the materiai modulus of reinforcement and V, is the fraction of concrete in unit volume. The material modulus due to reinforcement is evaluated by considering material modulus of the ~info~ment about its axis times its volume fraction and then transformed to the concrete axis. It is assumed that the presence of reinforcement will not
D. BAOCHIef al.
66
affect the plastic-fracture modulus of plain concrete. Equation (14) can be written in a matrix form as {da} = [C]{de}.
tuting these strains into eqn (16) yields
(14b)
Applying eqn (1) to the plate element, using eqns (12x14) the mass, stiffness and load matrices for a moderately thick plate element become and pATA dV, Y
[M’] =
[K’] =
s”
{dF’} =
B=CB dV,
dqA ‘dS.
(15)
To derive the nonlinear shear factor, approximate forms of the following equilibrium equations are used: a dM,, --de,,=0 axi
2.2. Derivation of mass, stiffness, and load matrices for a thick shell at any time step
a dMzz --dQ,,=O,
ax2
where h/2 dMii 1
- x, Ciiiide, dx,
;
s -h/2 and
s hi2
dQo =
where a, 6 are the length and width of the plate element; the factors 0, and G2 are the nonlinear shear factors in the directions of the curvatures rc, and K~. It can be observed that these factors can also be defined as the ratio of bending rigidity to the shear rigidity based on tangential material moduli.
c,,, dc, dx, , i = 1,2;
-h/2
and strains dc,, dc, are calculated using eqns (5) and (6) in terms of the unknown constants, {a}. Substi-
Consider the shell element in Fig. 2(a) with curvilinear axes x,, x2 and x3. Figure 2(b) shows an axisymmetric shell in which do is the angle in the meridional plane between the apex and the reference point in the shell at which the radii of curvatures R, and ROin the meridional direction (x, ) and circumferential direction (x2) are considered. Expressions for the displacement, dpi, for the thick shell are given by eqn (4) where d:, and d,ui are functions of x, , x2 only and d@; = 0 for small displacement theory. The straindisplacement relations for a thick shell element including the effect of transverse shear deformation can be readily found elsewhere [2,8]. The incremental displacement functions for a general axisymmetric shell with nonaxisymmetric load-
b) Fig, 2. (a) Deformed shape of a thick shell element; (b) an axisymmetric shell element with four d.o.f. at each node.
Dynamic response of thick plates and shells
ing can be represented as dp, = dp,,, cos nfl, +dM,,cos4
-rdQ,,=O
dp2 = dpg sin n0, dp, = dp,, cos no,
(1%
$-trdQ13)+n
dQ23
I
where
-
G,n = $,. -x3
(af,+aix,)
where hi* C,,kld6k,dX3,i,i=l,2
dNii =
k,l=l,3.
(21a)
s -h/2
G2n= 4,
-xl
(ai+a:x,)
4x:
( )
x 1--
-5
3h2
x(-n
r
dp3,+d0,,“sin4);
dp3,=a~+a~x,+a~x: + ai xi + a:‘~:;
d:,” = a1:+ aA0x,; d:,” = af’ + ai* x,;
dpi, is the incremental displacement for Fourier component, n; 4, is the incremental midplane displacement for Fourier component, n; ai are the constants for Fourier component, n and i = 1,. . . ,13; h is the thickness of the shell; and R, and r are defined in Fig. 2(b). At any stage of loading, the 13 constants given by eqn (19) for each Fourier component, n, are evaluated as follows. (1) Substituting the eight nodal degrees of freedom, defined in Fig. 2(b), in eqns (19), will result in
Evaluating the shell strain-displacement relations [2,8] in terms of variable x, and a and substituting them into eqn (21), yields four equations from the fist two equations of (21) by equating the coefficients of the variable x, for each equation to zero. The third equation of (21) will relate the constants a to the distributed load. For the purpose of calculating the constants, moments dM2*, dM,, are assumed to vary linearly in the x, direction. As in thick plate, the effect of variation in the moduli through the thickness of the shell is ignored by considering average moduli. The resulting five equations can be expressed in matrix form similar to eqn (10). (3) Combining eqns (20) and the above five equations results in an equation similar to (11). Solving this will yield the 13 constants in terms of the eight nodal degrees of freedom and the distributed load. After evaluating these constants and expressing the distributed load in terms of a cosine and sine series, eqn (19) can be written for each Fourier term in matrix form as IdKj =
[Al{@ij+[AAl{dqnIt
where (Si} are the nodal displacements and {dq,,} is the distributed load for the nth Fourier component. Using eqn (22), the strain-displacement relations for each Fourier component take the form:
W,} = PI{%} + Pl{dqnIv which is similar in form to eqn (7) except that [H’] is an 8 x 13 rectangular matrix. (2) The incremental stress-strain relation and moment and shear relations given in eqns (8) and (14) are now substituted in the following equations of equilibrium for a given Fourier component n for axisymmetric shell [2]:
&trdMll)+n
dM,2
I
-dM,,cos4
-rdQ,,=O
(22)
(23)
where {de,}r= {dCun, dC22nrdc,,, dt,3n, dcr,,}. The remaining procedure is similar to that for the thick plate element. It should be mentioned that under a nonaxisymmetric loading, the plastic moduli [C] changes with the angle 0 and as such, decoupling of eqn (3) for each harmonic cannot be done readily. Since integration of the element from 0 to 2x at each step is extremely time-consuming, three approaches are suggested, namely: the average modulus approach; the average effective stress approach; and the approach in which the stress-strain relations are expanded in a Fourier series. However, due to decoupling difficulties, the first approach was applied
68
D. BAGCHIet al.
herein. Thus an average modulus independent angle, 8, can be calculated as
of the
c =; i c,, ’
k=l
(24)
where C, is the actual modulus at an angle, 8,; and t is the number of points along the circumference where the actual modulus is calculated. As in the case of the thick plate element the incremental stiffness, mass and load matrices for a moderately thick shell element can now be obtained from eqn (15) for each Fourier component n. In deriving the nonlinear shear factors for the shell element, the following approximate forms of the shell equilibrium equations for moments are used: &dY,,)-d&cos+ I
-rdQ,,=O n dMz2 - r de,, = 0,
(25)
where dA4,i, dQ13are defined before in eqn (8). Evaluating dM,, dQa from the stress-strain relations and strain-displacement relations [2,8] in terms of variable x, and a and substituting in eqns (25), yields
where @,, and 9,, are the nonlinear the x, and x2 directions. 2.3. Solution procedure
The incremental global (structure) mass, stiffness and force matrices, [Ml, [K] and [F] are determined by assembly of the element mass, stiffness and force matrices [eqn (15)] after a suitable coordinate transformation. The dynamic equation of motion at a time step t [eqn (3)] is solved numerically by means of a backward finite difference scheme and a Gaussian integration procedure [ 1, 151. For example, the numerical solution of eqn (3) at a particular time step t, is found based on solutions at previous time steps (t - 2) and (t - 1). Thus, for a = 0, one can write for a small time interval At:
&JW’sJ =
{AF,
+ 1 +
-&It-2 dM,, cos 4
and
where oi, = ai +2azl+ 3ai1* + 4ai313; and I is the length of the element in the x, direction. For the case of a cylindrical shell, R, = CC and rj~= 90”, eqn (25) becomes h/2 6c x: Cm, dx3 ,
shear factors in
(& +1)Wi4
&fl{b-J (At)2
&-,I.
(28)
Depending on the nature of loading, it may be necessary to modify eqn (28) for the first two time steps. At the end of every time step solution, it is necessary to modify the mass and stiffness matrices to reflect changes in the stress-strain relation and in the shear rigidity. For a particular time t, the total strains and stresses are determined from the summation of incremental strains and stresses calculated previously. In general, the solution has been found to be stable (for example, up to the time before yielding of the steel in a reinforced concrete structure, say). However, experience has shown that instability in the solution can be circumvented by considering smaller time steps. A computer flow chart for the incremental solution is shown in Fig. 3.
3. DISCUSSION
OF RESULTS
The convergency of the solution was examined. Figure 4 shows the rapid rate of convergence of the result for the first natural frequency with increase in the number of elements for a cylindrical shell simply supported at both ends. Another comparison is made in Fig. 5 for the radial deflection response for a simply supported axisymmetric shell under suddenly applied ring load; good correspondence is shown with
Dynamic
READ:
response of thick plates and shells
STRUCTURAL GEOMETRY, MATERIAL PROPERTIES, BOUNDARY CONDITIONS 6 LOAD FUNCTION
I
4
FORM MASS h
ELASTIC STIFFNESS MATRICES OF EACH ELEMENT @ t=o
I FORM
GLOBAL MASS [Ml AND STIFFNESS MATRIX [Kl
<
1
FORM LOAD MATRIX
(Fj
SOLVE EQUATION (28)
FROM INPUT
LOAD FUNCTION1
TO OBTAIN DISPLACEMENT, Apt
CALCULATE FOR EACH ELEMENT: (i) STRAINSAND STRESSES; (ii) NONLINEAR SHEAR FACTOR USING CURRENTMATERIAL
MODULI
. REVISE MATERIALMODULI BASED ON REVISED TOTAL STRAINS AND STRESSESIN THE ELENBNT
ALL TIME STEPS
FORM MASS 6 STIFFNESS - MATRICES OF EACH ELEMENT BASED ON REVISED MODULI
PRINT RESULTS
Fig. 3. Flow chart for the incremental
results from an exact solution [la], and those from a finite element solution [ 171. 3.1. Results from elastic solution Deflection results at the center of a clamped square thin plate (h/a = 0.01) were compared with results from the classical solution (3) which neglects shear effects; the comparison shows that the results from the present solution are about 10% higher than those given by [ 181;this indicates that shear has an influence on the deflections even in thin plates. Deflection results for a simply supported relatively thick plate (h/a = 0.2) are only 2% higher than those derived using an eight-node thick plate element [l]. Further-
solution.
more, results for a clamped square plate show that transverse shear and rotatory inertia have no significant effect on the translational mass (M). However, they do have a significant effect on the natural frequency as the (h/a) ratio increases. [A decrease of 7% is noted for h/a = 0.1 (shear factor Q = 0.014) and a decrease of 21% for h/a = 0.3 (shear factor @ = 0.128).] The effect of the shear factor 4~ on the radial deflection of a short (I = 12.7’) and thick (h = 5’) cylindrical shell (radius r = a = 12’) was examined. The shell was fixed at the base and subjected to an axisymmetric load and then to a nonaxisymmetric edge load P at the upper free edge of the shell. For
D. BAGCHIet al.
IQ
I
x 14oDD-
r
5 12000-
3 10000z
SOOO-
t, f
6000-
(Ref. 17)
2 4000I Ia- 2DoD04 NO. OF ELEMENTS
USED
Fig. 4. Convergence of solution with increase in number of elements (mode, m = 0).
the axisymmetric load the present solution shows an increase of 10% in the radial deflection over the values derived from the classical solution [19] which neglects the shear effect; the increase in deflection is of course due to the shear effect; however, this percentage does not vary significantly with changes in @. For the nonsymmetric load (P cos 0) the shear factor @was found to be very significant, as expected; a 40% increase in deflection was observed when the shear factor increased from @ = 0.00084 to 0 = 0.337.
compressive strength of concrete = 4000 psi with longitudinal, circumferential and diagonal steel reinforcements with ratios of 1.5, 1.5 and l.O%, respectively. It is assumed that the shell is fixed at its two ends but axially free to move at its top end, with a length of 60 ft, wall thickness of 5 ft, and a radius of 60 ft. It is observed from Fig. 6 that when the longitudinal and circumferential steels yield and the concrete cracks near the mid-height of the shell (point U on curve OW, Fig. 6) the resulting loss of stiffness effects a rapid increase in the deflection (point V on curve OW, Fig. 6); continued loading triggers complete yielding of the longitudinal steel. This causes the transfer of load to the diagonal reinforcing steel which can resist the load in both the circumferential and longitudinal directions; also, the rate of increase in the radial deflections with time decreases as shown by curve VW in Fig. 6. It is evident that by increasing
ii SUDDENLY APPLIED RING(LINE)LOAl
3.2. Results from nonlinear-plastic solution Results for the nonlinear dynamic behavior of a reinforced concrete thick shell subjected to an axisymmetric impulsive load are shown in Fig. 6. Such loading may be due to an impact of seismic event. The translational mass matrix has been magnified by a factor of 10 to account for the effect of the surrounding water mass. The nonlinear radial deflection-time response is based on thick-she11 theory and on the plastic-fracture behavior of concrete as described in[lO]. It is assumed that the
TIME
t(* 10%~~)
Fig. 5. Radial deflection-time response of shell under suddenly applied ring load.
Diagonal Steal at bottom edge Started Yielding
Yielding of Circumf. 6 La&: Steel in middle of Shell, and Crac mg of Concrete
Fig. 6. Inelastic radial deflection-time
response at location A of reinforced concrete shell under axisymmetric ramp load.
Dynamic response of thick plates and shells
TIME
71
increase in the strain cl2 with the shear stress u12 remaining sensibly constant. Subsequent severe cracking of the concrete leads to a transfer of load to the diagonal reinforcing steel and its eventual yielding. After this, the increase of stress with strain is essentially due to material hardening effect. The same cylindrical shell was analyzed for a nonaxisymmetric load (P cos 0) as shown in the inset of Fig. 7; the shell was assumed to be fixed at one end and free at the other end. The results in Fig. 7 are based on Bazant’s plastic-fracture constitutive law [lo] and use the average modulus concept explained earlier. The radial deflection results for point B near the free edge show significant increase after the initiation of yielding in the diagonal reinforcing, accompanied by extensive cracking in the concrete. The shear stress-strain (u~~/c,~) relation for point A on the shell (see inset in Fig. 7) in Fig. 8 reveals an initial drop in stiffness at the initiation of cracking in the concrete followed by a subsequent
in seconds
Fig. 7. Inelastic radial ddlection-time response at location B of reinforced concrete shell under nonaxisymmetric first cosine ramp load.
the percentage of diagonal reinforcing to its limit value, the radial deflections at the center of the shell will be significantly reduced. Results for the stress-strain relations for the concrete were also determined but are not shown herein for brevity. The corresponding circumferential stress-strain (a2r/czz) relation was significantly nonlinear, exhibiting degrading stiffness until the final load was reached at a stress a,, = 5640 psi for a strain cu = 0.0075. The axial stress-strain (ur, /L,,) relation was essentially linear, perhaps due to the triaxial nature of compression in concrete. The results also show that transverse shear stress-strain (uu/+) relation is linear; however, upon yielding of the circumferential reinforcing and initial cracking, continued loading will incur plastic flow in the material, i.e. significant
Cracking
0
.OOl
I
222
444 LOAD
666
,
008
I
1110
(P) in psi
Fig. 9. Variation of shear factor with increase in load.
of Concrete
.002 SHEAR
Fig. 8. Shear stress-strain
0
.003 STRAIN
.004
C,,
relation at location A at the inner surface of shell in Fig. 6 at angle B= 90”.
D. BAGCHI et al.
72
increase in stiffness due to the presence of the longitudinal and circumferential steels. A subsequent drop in stiffness is observed after yielding of the latter steels; again, due to the presence of the diagonal reinforcing steel there is a subsequent increase in
(4) Composite material such as reinforced concrete exhibits severe nonlinearity in its shear factor @ as a result of tension cracks; this state of affairs would lead to structural failure at a much lower load than that predicted by elastic analysis.
stiffness followed by a levelling off in stiffness when ail steels had yielded. When compared with the response of the structure in Fig. 6, it is observed from
REFERENCES
the results in Figs 7 and 8 that the load-carrying capacity of the shell shown in Fig. 7 is much less. The variation of the nonlinear shear factor @, defined earlier by eqn (27), with load is presented in Fig. 9 for the shells and loadings in Figs 6 and 7; for the latter shell, two values of 6 are considered, namely 6 = 0” and 90”. It is observed that the shear factor remains essentially unaffected by increase in the load; however, when the concrete cracks and the diagonal reinforcing steel had yielded subsequently, the shear rigidity decreased rapidly, causing a sharp increase in the nonlinear shear factor @. Furthermore, the difference in the values of Q, for location C in Fig. 7 for inclination 6 = 0 and 8 = 90” does not appear to be significant.
0. C. Zienkiewicz, The Finite Element Methods in Engineering Science. McGraw-Hill, London (1977). H. Kraus, Thin Elastic Shells. John Wiley, New York
(1967). D. Bushnell, Buckling of elasto-plastic shells of revolu-
tion with discrete elasto-plastic-dng stiffeners. Znt. J. Sot%& Sirwt.
ticity.
4. CONCLUSIONS An incremental for~isplacement fo~ula~on for the solution of the nonlinear dynamic problem of thick plates and shells under arbitrary loading has been introduced using a discretized finite element energy approach. It is shown that: (1) The plate and shell elements developed can predict fairly accurately the structural response of relatively thick plate and shell structures; the deformation of these efficient elements is represented by a minimum number of degrees of freedom. (2) Effects of rotatory inertia and transverse shear strain are accurately represented; furthermore, the shear factor is estimated from the tangent moduli at every increment of loading. (3) The three-dimensional plastic-fracture characteristics of heterogeneous and anisotropic materials, such as reinforced concrete and laminated composites, can be readily treated, including cycles of loading-unloading of the structure due to cracks. Therefore, the method does not rely on a postulated crack direction nor on the use of an approximate moment-curvature relationship.
12, 51-66 11976).
4. M. P. Bieneck add J. R. l%naio, Elastic-plastic behavior of plates and shells. Tech. Report, Contract No. DNAOOl-76-C-0125, Defense Nuclear Agency, Washington, DC (1976). 5. R. F. Jones and J. E. Roderick, Incremental solution procedure for nonlinear structural dynamic problem. J. Pressure Vessel Technol. 97. 95-100 (1975). 6. S. Klein, The nonlin~r dynamic analysis’of shells of revolution with asymmetric properties by the finite element method. Trans. ASME 169, 163-170 (1975). 7. D. K. Bagchi, Inelastic bending of beam including transverse shear. Int. J. Nurner. Meth. Engng 14, 1323-1333 (1979). 8. K. Washizu, Variation Methods in Elasticity and Plas-
Pergamon Press, Oxford (1968).
9. H. F. Rubenstein and W. C. Hurty, Dynamics of Structures. Prentice-Hall, Englewood Cliffs, NJ (1964). 10.Z. P. Bazant and Sung-Sin Kim, Plastic-fmctu~ng theory for concrete. Engng Mech. Div., ASCE 105, (EM3), 407-428 (1979). 11. R. Hill, The Mathematical Theory of Plasticity. Clarendon Press, Oxford (1950). 12. B. Budiansky and R. J. O’Connell, Elastic moduli of a cracked solid. Int. J. Solids Struct. 12, 81-97 (1976). 13. J. W. Dougill and M. A. M. Rid, Further consideration of progressively fracturing solids. Errgag Mech. Div., ASCE 106(EM5), 1021-1038 (1980). 14. Zvi Hashin, D. K. Bagchi and B. Rosen, Nonlinear behaviour of fiber composite laminates. NASA Contract Report prepared under Contract NASI-11284 (1973). 15. R. S. Burington, Handbook ofMathematical Tables and Formulas. 5th Edn. McGraw-Hill. New York (1973). 16. H. Reismann and J. Padlog, Forced, axi-symmetric motions of cylindrical shells. J. Franklin Inst. 284(S), 308-319 (1967).
17. S. Ghosh and E. Wilson, Dynamic stress analysis of axisymmetric structures under arbitrary loading. Rep. No. EERC 69-10, College of Engineering, University of California, Berkeley (1969). 18. S. Timoshenko and S. Woinowesky-Krieger, Theory of Plates and Shells, 2nd Edn. McGraw-Hill, New York (1961). 19. R. J. Roark, Formulae for Stress and Strain. McGraw-Hill, New York (1965).