- Email: [email protected]
Dynamic analysis of a pressure vessel
0045.7949/a 53.00 + 0.00 01~mPQtt~ple
Conlpcms & Stmctwex Vol. 29, No. I, pp. 161470, 1988 Flinted in Qnnt Britain.
DYNAMIC
ANALYSIS
OF A PRESSURE VESSEL
V. RAMAMURTI and P. SESHU Departmentof Applied Mechanics, Indian Institute of Technology, Madras 600 036, India (Received 13 Jdy 1987) AMract-The finite dilTerence method together with the Novoxhilov shell theory has been used to formulate the seismic response of a pmssure vessel. In the circmnferential direction, a Fourier series variation is used for loads and deformations. The static problem (i.e. internal pmssm’e problem) is solved using Potters’ scheme. Lumped mass method is adopted for mass matrix determination. Simultaneous iteration method is used to solve the unsymmetric, non-standard eigenvahte problem. Medalsuperposition technique is used to determine the time-domain response of the structure when subjected to an earthquake excitation.
NOTATION
reported that the modal superposition method is very efficient for earthquake response problems. In this method, the lowest few modes of free vibration of the structure are considered. The ~nt~bution of each of these modes to the structural response is determined. A simplified method of calculation of these bending moment per unit length Fourier coe&icnt of Me mode participation factors has been reported by &cmnferential wave number Tedesco et af. [2]. Using finite element analysis, Altes Potter blocks er al. [3] analysed the seismic behaviour of a reactor. shell loads per unit area In the present paper, treating the pressure vessel as a Fourier cwe&ients of P,. P, and P, Fourier coef&ient of K&hhoffs radial shear shell of revolution, the method of iinite differences is thickness of the shell applied to develop a totally compute&d model of reference thickness determining the seismic response in the time-domain. Fourier coefficient of membrane force per unit An efficient way of solving the resulting unsymmetric fen%h eigenvalue problem is described. It is well known that Fourier coethcient of Kirchhog’s transverse Shy for structures such as these under a horizontal earthmeridional, circumferential and normal dis- quake excitation, only the first one or two modes of placements on the middle surface the first circumferential harmonic (n = 1) contribute Fourier coefficients of U, V, W significantly to the response. In the present analysis, eigen vector-straight problem the first four modes of the structure corresponding to ground displacement vector relative displacement vector n = 1 are taken for calculation of the response. eigen vector-transpose problem The detailed dimensions of the pressure vessel and non~rn~~o~ me~dio~l coordinate discretisation adopted, are shown in Fig. 1. It is a inclination change at discontinuity Fourier coefbcient of rotation in the meridional multi-layered pressure vessel with a flat-head, long direction cylindrical portion and a spherical bottom, with a mass density cylindrical skirt branching off from the cylindrical reference stress level body. The static stresses when the vessel is subjected hysteritic damping factor to internal pressure are determined. The dynamic natural frequency (rad/seo) response of the structure to an earthquake excitation non~~nsio~ curvature appropriate diagonai matrices at the junction is obtained (based on Koyna earthquake 1967, node India). non-dimensional interval length Novoxhilov’s shell theory is used to determine the governing differential equations of equilibrium. For the purpose of analysis, the method of finite differences is used with Fourier series ~rn~~tion, 1. ~ODU~ON for all the pertinent variables, in the circumferential D&rrnination of the seismic response of a pressure direction. vessel is of real practical importance. In a numerical approach, the modal solution techniques are gener2. MATDEMATICALFDDMULATION ally adopted for determming the structural response The whole structure is treated as a shell of revolugiven the excitation input, Based on a critical evaluation of the modal methods, Anagnostopoulus [l] has tion. An element of the shell considered for analysis refmnce length iXdii&nts of the Novozhilovtheery Young’smodulus stiffnessmatrix massmatrix
161
V. RAMAMURTI and P. SESHU
162
Fig. 2. Typical shell element.
where d (I’=@ The non-dimensional Fourier coefficients pt , ps and the intensities of pressure in the three respective directions for static loading. They are replaced by inertia forces per unit area in the three respective directions for dynamic problems. Also, p, represent
iv=
Wl{z’~+ LwZ~~
(2)
where
Fig. 1. Structure discretisation (dimensions in mm).
is shown in Fig. 2. The governing motion when the shell is subjected loading are given by [4]
equations of to a general
a,un + a& + a,24 + qu’ + a,o + UhW’+ u,w +u,m;
(1,011’+ u,,u + a,+” + u,+’ + +u,,mc=
+u,~w’+u~~w
u,,u’+u,u
+a,,w
u,u’+u~,u
The coelkients a, and 6, are listed in [5]. The boundary conditions in their most general form can be stipulated as
+ u9mt = -pc
+u*,o”+u,,u’+a,,u
+
,, q4u
+
u,,w
-pO
rnllY> + [nIbI = VJv +u,w”
u2,w + a,,m; + uam; + u,,m, = -p, +u,,u
+u,w’+u,,w
where [a] and [/i] are appropriate diagonal matrices and {L} is a column matrix specifying the stress resultants and displacements on the boundary. 2.1 Continuity conditions
+u33wn +u,,mc=O,
(3)
(1)
Whenever the shell geometry changes, e.g. from flat head to cylinder or from cylinder to sphere, the
Dynamic analysis of 8 prca8luc vessel
163
Fig. 3. Continuity conditions. following relations can be easily real&d with reference to Fig. 3.
{Y)‘=[WY~@I’ =
(4)
Fw~-
(5)
main cylindrical portion and skirt. Two nodes on the structure are taken to constitute a partition for Potters’ scheme solution [6] with a uniform bandwidth of 16. The formation of Potter blocks is explained below. For any two successive nodes i and i + 1, eqn (7) gives
where
rcQscp0
-sincp
01
or Using eqn (2), eqns (4) and (5) can be combined to give
[HI+@‘I+ + wl’[‘pl- [~lvl-l~z~-V’lWl- @I- = (0). 64 2.2 Dlyerence equations lP‘48x (I
Using finite difference scheme, the governing equations given in eqn (1) can be recast [4], for a continuous geometry, as ordinary simultaneous equations in four unknowns. For a general node i
[Atl{z,+I} + IBil{zt)+ LcilIzi-II = bI*
C7)
For the first and last nodes, the forward and backward differences, respectively, are used. At a discontinuity node i, with
{z;}+=((zi+,)-{z,}+)/A+
(8)
IzIl- =(bJ- -h-,))/A-,
(9)
i.e.
+[pAl,{Z’+,} = {FORCE}j, (13) where [PA], [PSI, [PC] are the blocks used in Potters’ scheme and j the partition number. In the skirt portion, the alternating node numbers will affect the locations of [A,], [B,], [C,] in the blocks of eqn (13).
and Branch
1
Branch
eqn (6) can be fitted into the master form eqn (7). Using eqn (5), [C,] of eqn (7) at a node following the discontinuity is replaced by [C,][‘P] in order to retain only (z}-. 3. RANDWIDTH CONSIDEtRATION-CROUPING OF NODES
Branch
111
An
alternating scheme of node numbering is used in the skirt portion of the structure as shown in Fig. 1. The nodes are numbered alternately between the
Fig. 4. Branches of shell.
11
V. RAwn
164
and P. !bnu
0.5
-0.5 0
2
L
10
8
6 Time in sac
Fig. 5. Accelerogram of Koyna earthquake of 11Dec. 1967 (horizontal component along axis of Koyna Dam).
4.1 First approach
4. JUNCIION NODE
At node 249 (Fig. 1) of the structure, three branches of the shell meet. With reference to Fig. 4 the compatibility of displacements and forces yields [4]
This is a complete adaptation of the method suggested in [41. From Potters’ scheme, {Ii{:}=
-[PI-‘I{:“‘}+{,-‘},
(18)
r.Bl{~“‘~ + [‘IHY”‘)= uN’p”1~z”~+ rvl{y”1 = [B][Y’]{z’) + [q]{ Y’}
(14)
wh;s;njg;$‘$n
number.
wl{Y”‘}+ [tll(z”‘)= [Blw”~~y”~ + hl{z”~
Z’ Z”
11
+[B1Wl{y’~ + rttl{z’). (19 Using eqn (2), the difference equations of eqns (14) and (15) can be found to be
is eliminated from eqns (16) and (17) to result in the following equation at the junction
[Ql{zi+~)+[RI{z”‘)
w1l4.*
=[S”]{Z*‘}+ [P]{zfl_,}
where [DD l] and (002)
=[S’]{z’}+ [T’l{zf_,} [M]{z,+,}+[N]{z”‘}
(16) the elimination process. Considering
-[X”l{z”}+[XX”l{zf’_,}
+[X’l{z’} +[XX’l{zf_,}, where i = 249, the junction
11” *x, = P2h.‘, 1 1
are obtained as products of
node i + 1, eqn (7) gives
[Ai+‘l{zi+*j+[Bi+‘l{zi+‘}
(17)
+[ci+‘l{zl”) = {gi+l},
node.
0.6 ; ._ Z 3 8 8
0
-0.5 I
0
2
I ,
I
I
I
I
6
6
10
Time
(19)
in sac
Fig. 6. Accelerogram-used
for analysis.
C20)
Dynamic analysis of a pressure vessel ---She\\
-.-FFEM
theory - First
-Sheti
theory - Second
approach
165
4.2.1 Nodes (i - 1)’ and (i - 1)“. From eqn (7),
approach
t~:-,l(z’}+lBf-*l(zf-,} +Kz,lH-,~
t-
= Id-ll
tw
[A1’_,]{z”}+[Bf~,l(rft,} +[c,‘i.,t(z~‘-tf=(g:f_,). (23) On grouping,
/-86
Nfmm’
12xI
s. M. wuc 83NlmmZ
(24) ~03Nlmm2
-.-
FEM -
) ---_ W/
_ c--
/
---
Shell theory-first Shell theory-Second
approach approach
>
Fig. 7. Meridional stress (static problem). or, on grouping i and i + 1,
S. M.Value 176N/mm*
4.2 Second approach In order to avoid the tedious mathematical manip ulations involved in eliminating
t77NImmZ r
21
21
(1 ii zII
’
Z’I
z
i
III
was identified as the state variable at the junction node i. This greatly simpMkd the mathematical formulation and is readily applicable even when there are many branches meeting at the junction.
Fig. 8. Circumf’krentialstress (static problem).
V. RAMAMuRn and P. SPSHU
166
/II\
/
\
I
t-
/
I
i
I---_
1
\ ,,
\ \ \ \
I / / ‘/
/
1’
/ I
(d)
1
Fig. 9.
Mode shapes (n = 1). (a) w = 22.37 rad/sec. (b) o = 158.9rad/sec. (c) w = 386.2rad/sec. (d) w = 556.3 rad/sec.
4.2.2 Node i, the junction. Writing eqns (16) and (17) in matrix form,
[ i]tz,+
1)4X1+ [ 1:
-[T’]
Ii]{
I}
5. MASS MATRIX
wll2 x 12
PA112 x4
+
12;
one of programming the computer for variable block size. Also, the problem of many branches meeting at a junction can be easily tackled by the present approach whereas the first approach becomes more labofious.
0
-[T”] 0 -[XX”] -[XX’] m-l12x8
(25)
Lumped mass method is used. For dimensionalisation chosen,
(27)
PC=
4.2.3 Node i + 1. From the field eqn (7),
the non-
Inertia force per unit area is given by Pt=ptid
therefore
(1 u0=
E
u,
(28)
= (gi+I 14x I * (26) (29)
12x I
It can be seen that the problem area has been shifted from tedious mathematical manipulation to
where F = t/to.
167
Dynamic analysis of a pmssurc vessel ,
_/
I
I----
\
>
.’
cl \
>
,’
d
k
(a)
r
t
(b)
I
/
r
/
td)
Fig. 10. Mode shapes (n = 2). (a) o = 333.9 rad/sec. (b) o = 364.5rad/sec. (c) o = 465.1 rad/sec. (d) o = 618.1rad/sec. The quantity within the brackets of eqn (29) represents a typical mass matrix element. From the formulation of the problem itself, every fourth element of the mass matrix is zero. Also, at the discontinuity nodes and at the junction node the mass matrix elements are taken to be zero in accordance with the governing equations at these nodes. With this semi-definite, diagonal mass matrix the unsymmetric eigen value problem is solved in its non-standard form by simultaneous iteration method [7,8].
Let
{WWl{W = PI.
(32)
If m eigen vectors are taken, [D] is an unsymmetric matrix of size (m x m). Decomposing [D],
where
6. ITERATION SCHEME
Two problems are identified:
rfmw = ~2wl{~l [KIT{ V} = WQf]‘{ V).
(30) (3 1)
The bi-orthogonality condition is used to orthonormalise the eigen vectors with respect to the mass matrix as is shown below. C.A.S. 29/l-L
L 0 Let
{Xl) and vectors and let
(X2)
{Xl] =
kJ be the ortho-normalised
{Whl
IX21 = 1 V
hl;
(34 (35)
V. RAMAWJR~~ and P. SESIU
168
1--
(a)
(b)
(d)
Fig. 11. Mode shapes (n = 3). (a) w = 926.5 rad/sec. (b) w = 953.9 rad/sec. (c) o = 979 rad/sec.
(d) o = 994.9 rad/sec.
7. SEISMIC ANALYSIS
therefore
The governing equation when the structure subjected to an earthquake is
WIW
the identity matrix. Using eqn (33),
h1~‘[~,1L1%1= [II,
(36)
from which ]s,l=
IJTT
[+I= L-r.
(37) (38)
Using eqns (34), (35), (37) and (38) the orthonormalised vectors can be found and used in the next iteration. As the stiffness matrix is already available in the form of blocks, the out of core solution offered by Potters’ scheme is used in solving eqns (30) and (31).
+[clIlj,I +[m~,~ = 4~1~lf,).
is
(39)
The details pertaining to the Koyna earthquake of 1967 [9] are used. The structure is assumed to be uniformly subjected to the ground acceleration in the horizontal plane and the corresponding inertia force is impressed on the structure at all the nodes. In the actual accelerogram as charted in [9] and reproduced in Fig. 5, the peak points are identified and a linear variation of acceleration is used between two consecutive peaks. The resulting graph, which is used for seismic input, is shown in Fig. 6. Modal superposition technique [lo] is used to solve eqn (39). As was proved in [ll], for this type of structure, only n = 1 modes are excited by the horizontal earthquake disturbance. So, the first four modes pertaining to n = 1 are taken for solution. As the fundamental natural frequency was 22.37 rad/sec for n = 1, a time step of about l/50 set is taken for response determination.
Dynamic analysis of a pressure vessel
169
Fig. 12. Cantilever model (dimensions in mm).
For hysteritic damping in the material, an l/r value of 300 (for steel) [12] is taken. An arbitrary value of 0.5 for the damping ratio is set as the upper limit in the fifth mode and Rayleigh damping curve fit is used to interpolate at the other natural frequencies.
8. RESULTS
AND DISCUSSIONS
8.1 Static analysis For a uniform internal pressure of 220 atmospheres, the static stresses are determined. The stresses as obtained by the two different approaches at the junction node are compared with those obtained by the finite element method (using the shape function which is a linear polynomial in radial and axial directions) in Figs 7 and 8. Considering the pressure vessel as a closed ended thin cylinder subjected to internal pressure, the stress values predicted by elementary strength of materials are calculated and are shown in Figs 7 and 8. The ratio of circumferential to meridional stress in the main cylindrical body is obtained as 2.13 whereas strength of material predicts a value of 2.0.
-10. 0
8
10
Time in set
I 12
6 4 Time in set
2
la)
8.2 Dynamic analysis The eigen-pairs are determined for three values of n (1, 2 and 3). The first four mode shapes for each value of n are shown in Figs 9-l 1. A crude approximation of the pressure vessel as a cantilever beam is shown in Fig. 12. This is used to validate the results for n = 1. Since u displacement is very small compared to I) and w (ratio varying from & to h) and the latter two are of the same order (because n = I), only w displacements need be considered for potential and kinetic energy terms. Total mass of the structure, M = 360 x 10’ kg Mass of the lid, m = 28.5 x 10’ kg Length (span) of the cantilever, L = 23.0 m
I
-7.51 (bl
10
Fig.
13.
Relative
I 11
acceleration
excnation.
I
responx+horizontal
V. RAMAMUIUIand P. SESHU
-501
I
0
I
I
2
L
I
I
I
Time
I
I
I
6
6
10
I
I
12
in ret
Fig. 14. Dynamic stress plot-horizontal
Second moment of area, Z = 2.75 m4 Young’s modulus for steel = 2.058 x 10” N/m*
It is to be noted that the span of the equivalent cantilever beam (determined only very approximately) can strongly influence the value of the natural frequency as obtained above. In the mode shapes it is observed that a node for axial displacement (a) corresponds to an anti-node for radial displacement (w) and vice versa. Also for n >O,uisobservedtobeoftheorderof(-w/n).The participation of the lid is minimal for II = 3. 8.3 Seismic response The relative acceleration response of the structure is shown in Fig. 13. In Fig. 13(b), the time scale is chosen so as to indicate the time period of the response. When the ground acceleration is given in the horizontal plane i.e. D and w directions of the structure, a peak response of 8.8 g (i.e. a dynamic magnification factor of 14.7) is obtained at node 25 on the cylindrical portion, 800 mm from the flat head. The dynamic stresses in the material of the structure are also obtained. Node 232,500 mm away from the support on the skirt, is found to have the maximum stress. The dynamic stress plot for this node is shown in Fig. 14. In the main cylindrical body of the structure a peak stress of about f 12.20 kgf/mm* is obtained. 9. CONCLUSlONS
The dynamic response of a pressure vessel to a typical earthquake has been obtained. A modification
excitation.
to the existing method of treatment of a junction node is presented. It has been observed that the skirt portion of the structure is the maximum stress zone. The stress versus time plot over a period of 12.5 set for a node on the skirt is presented. REFERENCES
1. S. A. Anagnostopoulos, Wave and earthquake response of offshore structures; evaluation of modal solutions, AXE, J. Struct. Div. loS, 2175-2191 (1982). 2. J. W. Tedesco, C. N. Kostem and A. Kahtins, Rapid calculation of mode participation factors for circular cylindrical shells. Comput. Sfruct. Zo, 509-515 (1985). 3. J. Altes, H. Grafii and D. Koschmieder, Investigation of the earthquake behaviour of the research reactor FRJ-2 (Dido). Report No. JUEL-SPEZ-222 (1983). 4. B. Budiansky and P. P. Radkowski, Numerical analysis of unsymmetric bending of shells of revolution. AZAA Jnl 1, 1833-1842 (1963). 5. V. Ramamurti and R. S. Alwar, Stress analysis of tube mills. J. Strain Anal. 8, 200-208 (1973). 6. M. L. Potters, A matrix method for the solution of a second order difference equation in two variables. Mathematische Centrum, Amsterdam, Holland, Report No. MR 19 (1955). 7. Alan Jennings, Matrix Computation for Engineers and Scientists. John Wiley, London (1977). 8. V. Ramamurti and 0. Mahrenholta, The application of the simultaneous iteration method to flexural vibration problems. Int. J. Mech. Sci. 16, 269-283 (1974). 9. J. Krishna and A. R. Chandrasekharan, Elements of Earthquake Engineering. Sarita Prakashan, Meerut (1976). 10. K. J. Bathe and E. L. Wilson, Numerical Methods in Finite Element Analysis. Prentice-Hall, New Delhi (1978). 11. E. L. Wilson, Earthquake analysis of reactor structures. In Symposium on Seismic Analysti of Pressure Vessel and Piping Components, pp. 53-73. ASME, New York (1971). 12. W. K. Wilson, Practical Solutions of Torsional Vibration Problems, Vol. 2. Chapman & Hall, London (1963).
Conlpcms & Stmctwex Vol. 29, No. I, pp. 161470, 1988 Flinted in Qnnt Britain.
DYNAMIC
ANALYSIS
OF A PRESSURE VESSEL
V. RAMAMURTI and P. SESHU Departmentof Applied Mechanics, Indian Institute of Technology, Madras 600 036, India (Received 13 Jdy 1987) AMract-The finite dilTerence method together with the Novoxhilov shell theory has been used to formulate the seismic response of a pmssure vessel. In the circmnferential direction, a Fourier series variation is used for loads and deformations. The static problem (i.e. internal pmssm’e problem) is solved using Potters’ scheme. Lumped mass method is adopted for mass matrix determination. Simultaneous iteration method is used to solve the unsymmetric, non-standard eigenvahte problem. Medalsuperposition technique is used to determine the time-domain response of the structure when subjected to an earthquake excitation.
NOTATION
reported that the modal superposition method is very efficient for earthquake response problems. In this method, the lowest few modes of free vibration of the structure are considered. The ~nt~bution of each of these modes to the structural response is determined. A simplified method of calculation of these bending moment per unit length Fourier coe&icnt of Me mode participation factors has been reported by &cmnferential wave number Tedesco et af. [2]. Using finite element analysis, Altes Potter blocks er al. [3] analysed the seismic behaviour of a reactor. shell loads per unit area In the present paper, treating the pressure vessel as a Fourier cwe&ients of P,. P, and P, Fourier coef&ient of K&hhoffs radial shear shell of revolution, the method of iinite differences is thickness of the shell applied to develop a totally compute&d model of reference thickness determining the seismic response in the time-domain. Fourier coefficient of membrane force per unit An efficient way of solving the resulting unsymmetric fen%h eigenvalue problem is described. It is well known that Fourier coethcient of Kirchhog’s transverse Shy for structures such as these under a horizontal earthmeridional, circumferential and normal dis- quake excitation, only the first one or two modes of placements on the middle surface the first circumferential harmonic (n = 1) contribute Fourier coefficients of U, V, W significantly to the response. In the present analysis, eigen vector-straight problem the first four modes of the structure corresponding to ground displacement vector relative displacement vector n = 1 are taken for calculation of the response. eigen vector-transpose problem The detailed dimensions of the pressure vessel and non~rn~~o~ me~dio~l coordinate discretisation adopted, are shown in Fig. 1. It is a inclination change at discontinuity Fourier coefbcient of rotation in the meridional multi-layered pressure vessel with a flat-head, long direction cylindrical portion and a spherical bottom, with a mass density cylindrical skirt branching off from the cylindrical reference stress level body. The static stresses when the vessel is subjected hysteritic damping factor to internal pressure are determined. The dynamic natural frequency (rad/seo) response of the structure to an earthquake excitation non~~nsio~ curvature appropriate diagonai matrices at the junction is obtained (based on Koyna earthquake 1967, node India). non-dimensional interval length Novoxhilov’s shell theory is used to determine the governing differential equations of equilibrium. For the purpose of analysis, the method of finite differences is used with Fourier series ~rn~~tion, 1. ~ODU~ON for all the pertinent variables, in the circumferential D&rrnination of the seismic response of a pressure direction. vessel is of real practical importance. In a numerical approach, the modal solution techniques are gener2. MATDEMATICALFDDMULATION ally adopted for determming the structural response The whole structure is treated as a shell of revolugiven the excitation input, Based on a critical evaluation of the modal methods, Anagnostopoulus [l] has tion. An element of the shell considered for analysis refmnce length iXdii&nts of the Novozhilovtheery Young’smodulus stiffnessmatrix massmatrix
161
V. RAMAMURTI and P. SESHU
162
Fig. 2. Typical shell element.
where d (I’=@ The non-dimensional Fourier coefficients pt , ps and the intensities of pressure in the three respective directions for static loading. They are replaced by inertia forces per unit area in the three respective directions for dynamic problems. Also, p, represent
iv=
Wl{z’~+ LwZ~~
(2)
where
Fig. 1. Structure discretisation (dimensions in mm).
is shown in Fig. 2. The governing motion when the shell is subjected loading are given by [4]
equations of to a general
a,un + a& + a,24 + qu’ + a,o + UhW’+ u,w +u,m;
(1,011’+ u,,u + a,+” + u,+’ + +u,,mc=
+u,~w’+u~~w
u,,u’+u,u
+a,,w
u,u’+u~,u
The coelkients a, and 6, are listed in [5]. The boundary conditions in their most general form can be stipulated as
+ u9mt = -pc
+u*,o”+u,,u’+a,,u
+
,, q4u
+
u,,w
-pO
rnllY> + [nIbI = VJv +u,w”
u2,w + a,,m; + uam; + u,,m, = -p, +u,,u
+u,w’+u,,w
where [a] and [/i] are appropriate diagonal matrices and {L} is a column matrix specifying the stress resultants and displacements on the boundary. 2.1 Continuity conditions
+u33wn +u,,mc=O,
(3)
(1)
Whenever the shell geometry changes, e.g. from flat head to cylinder or from cylinder to sphere, the
Dynamic analysis of 8 prca8luc vessel
163
Fig. 3. Continuity conditions. following relations can be easily real&d with reference to Fig. 3.
{Y)‘=[WY~@I’ =
(4)
Fw~-
(5)
main cylindrical portion and skirt. Two nodes on the structure are taken to constitute a partition for Potters’ scheme solution [6] with a uniform bandwidth of 16. The formation of Potter blocks is explained below. For any two successive nodes i and i + 1, eqn (7) gives
where
rcQscp0
-sincp
01
or Using eqn (2), eqns (4) and (5) can be combined to give
[HI+@‘I+ + wl’[‘pl- [~lvl-l~z~-V’lWl- @I- = (0). 64 2.2 Dlyerence equations lP‘48x (I
Using finite difference scheme, the governing equations given in eqn (1) can be recast [4], for a continuous geometry, as ordinary simultaneous equations in four unknowns. For a general node i
[Atl{z,+I} + IBil{zt)+ LcilIzi-II = bI*
C7)
For the first and last nodes, the forward and backward differences, respectively, are used. At a discontinuity node i, with
{z;}+=((zi+,)-{z,}+)/A+
(8)
IzIl- =(bJ- -h-,))/A-,
(9)
i.e.
+[pAl,{Z’+,} = {FORCE}j, (13) where [PA], [PSI, [PC] are the blocks used in Potters’ scheme and j the partition number. In the skirt portion, the alternating node numbers will affect the locations of [A,], [B,], [C,] in the blocks of eqn (13).
and Branch
1
Branch
eqn (6) can be fitted into the master form eqn (7). Using eqn (5), [C,] of eqn (7) at a node following the discontinuity is replaced by [C,][‘P] in order to retain only (z}-. 3. RANDWIDTH CONSIDEtRATION-CROUPING OF NODES
Branch
111
An
alternating scheme of node numbering is used in the skirt portion of the structure as shown in Fig. 1. The nodes are numbered alternately between the
Fig. 4. Branches of shell.
11
V. RAwn
164
and P. !bnu
0.5
-0.5 0
2
L
10
8
6 Time in sac
Fig. 5. Accelerogram of Koyna earthquake of 11Dec. 1967 (horizontal component along axis of Koyna Dam).
4.1 First approach
4. JUNCIION NODE
At node 249 (Fig. 1) of the structure, three branches of the shell meet. With reference to Fig. 4 the compatibility of displacements and forces yields [4]
This is a complete adaptation of the method suggested in [41. From Potters’ scheme, {Ii{:}=
-[PI-‘I{:“‘}+{,-‘},
(18)
r.Bl{~“‘~ + [‘IHY”‘)= uN’p”1~z”~+ rvl{y”1 = [B][Y’]{z’) + [q]{ Y’}
(14)
wh;s;njg;$‘$n
number.
wl{Y”‘}+ [tll(z”‘)= [Blw”~~y”~ + hl{z”~
Z’ Z”
11
+[B1Wl{y’~ + rttl{z’). (19 Using eqn (2), the difference equations of eqns (14) and (15) can be found to be
is eliminated from eqns (16) and (17) to result in the following equation at the junction
[Ql{zi+~)+[RI{z”‘)
w1l4.*
=[S”]{Z*‘}+ [P]{zfl_,}
where [DD l] and (002)
=[S’]{z’}+ [T’l{zf_,} [M]{z,+,}+[N]{z”‘}
(16) the elimination process. Considering
-[X”l{z”}+[XX”l{zf’_,}
+[X’l{z’} +[XX’l{zf_,}, where i = 249, the junction
11” *x, = P2h.‘, 1 1
are obtained as products of
node i + 1, eqn (7) gives
[Ai+‘l{zi+*j+[Bi+‘l{zi+‘}
(17)
+[ci+‘l{zl”) = {gi+l},
node.
0.6 ; ._ Z 3 8 8
0
-0.5 I
0
2
I ,
I
I
I
I
6
6
10
Time
(19)
in sac
Fig. 6. Accelerogram-used
for analysis.
C20)
Dynamic analysis of a pressure vessel ---She\\
-.-FFEM
theory - First
-Sheti
theory - Second
approach
165
4.2.1 Nodes (i - 1)’ and (i - 1)“. From eqn (7),
approach
t~:-,l(z’}+lBf-*l(zf-,} +Kz,lH-,~
t-
= Id-ll
tw
[A1’_,]{z”}+[Bf~,l(rft,} +[c,‘i.,t(z~‘-tf=(g:f_,). (23) On grouping,
/-86
Nfmm’
12xI
s. M. wuc 83NlmmZ
(24) ~03Nlmm2
-.-
FEM -
) ---_ W/
_ c--
/
---
Shell theory-first Shell theory-Second
approach approach
>
Fig. 7. Meridional stress (static problem). or, on grouping i and i + 1,
S. M.Value 176N/mm*
4.2 Second approach In order to avoid the tedious mathematical manip ulations involved in eliminating
t77NImmZ r
21
21
(1 ii zII
’
Z’I
z
i
III
was identified as the state variable at the junction node i. This greatly simpMkd the mathematical formulation and is readily applicable even when there are many branches meeting at the junction.
Fig. 8. Circumf’krentialstress (static problem).
V. RAMAMuRn and P. SPSHU
166
/II\
/
\
I
t-
/
I
i
I---_
1
\ ,,
\ \ \ \
I / / ‘/
/
1’
/ I
(d)
1
Fig. 9.
Mode shapes (n = 1). (a) w = 22.37 rad/sec. (b) o = 158.9rad/sec. (c) w = 386.2rad/sec. (d) w = 556.3 rad/sec.
4.2.2 Node i, the junction. Writing eqns (16) and (17) in matrix form,
[ i]tz,+
1)4X1+ [ 1:
-[T’]
Ii]{
I}
5. MASS MATRIX
wll2 x 12
PA112 x4
+
12;
one of programming the computer for variable block size. Also, the problem of many branches meeting at a junction can be easily tackled by the present approach whereas the first approach becomes more labofious.
0
-[T”] 0 -[XX”] -[XX’] m-l12x8
(25)
Lumped mass method is used. For dimensionalisation chosen,
(27)
PC=
4.2.3 Node i + 1. From the field eqn (7),
the non-
Inertia force per unit area is given by Pt=ptid
therefore
(1 u0=
E
u,
(28)
= (gi+I 14x I * (26) (29)
12x I
It can be seen that the problem area has been shifted from tedious mathematical manipulation to
where F = t/to.
167
Dynamic analysis of a pmssurc vessel ,
_/
I
I----
\
>
.’
cl \
>
,’
d
k
(a)
r
t
(b)
I
/
r
/
td)
Fig. 10. Mode shapes (n = 2). (a) o = 333.9 rad/sec. (b) o = 364.5rad/sec. (c) o = 465.1 rad/sec. (d) o = 618.1rad/sec. The quantity within the brackets of eqn (29) represents a typical mass matrix element. From the formulation of the problem itself, every fourth element of the mass matrix is zero. Also, at the discontinuity nodes and at the junction node the mass matrix elements are taken to be zero in accordance with the governing equations at these nodes. With this semi-definite, diagonal mass matrix the unsymmetric eigen value problem is solved in its non-standard form by simultaneous iteration method [7,8].
Let
{WWl{W = PI.
(32)
If m eigen vectors are taken, [D] is an unsymmetric matrix of size (m x m). Decomposing [D],
where
6. ITERATION SCHEME
Two problems are identified:
rfmw = ~2wl{~l [KIT{ V} = WQf]‘{ V).
(30) (3 1)
The bi-orthogonality condition is used to orthonormalise the eigen vectors with respect to the mass matrix as is shown below. C.A.S. 29/l-L
L 0 Let
{Xl) and vectors and let
(X2)
{Xl] =
kJ be the ortho-normalised
{Whl
IX21 = 1 V
hl;
(34 (35)
V. RAMAWJR~~ and P. SESIU
168
1--
(a)
(b)
(d)
Fig. 11. Mode shapes (n = 3). (a) w = 926.5 rad/sec. (b) w = 953.9 rad/sec. (c) o = 979 rad/sec.
(d) o = 994.9 rad/sec.
7. SEISMIC ANALYSIS
therefore
The governing equation when the structure subjected to an earthquake is
WIW
the identity matrix. Using eqn (33),
h1~‘[~,1L1%1= [II,
(36)
from which ]s,l=
IJTT
[+I= L-r.
(37) (38)
Using eqns (34), (35), (37) and (38) the orthonormalised vectors can be found and used in the next iteration. As the stiffness matrix is already available in the form of blocks, the out of core solution offered by Potters’ scheme is used in solving eqns (30) and (31).
+[clIlj,I +[m~,~ = 4~1~lf,).
is
(39)
The details pertaining to the Koyna earthquake of 1967 [9] are used. The structure is assumed to be uniformly subjected to the ground acceleration in the horizontal plane and the corresponding inertia force is impressed on the structure at all the nodes. In the actual accelerogram as charted in [9] and reproduced in Fig. 5, the peak points are identified and a linear variation of acceleration is used between two consecutive peaks. The resulting graph, which is used for seismic input, is shown in Fig. 6. Modal superposition technique [lo] is used to solve eqn (39). As was proved in [ll], for this type of structure, only n = 1 modes are excited by the horizontal earthquake disturbance. So, the first four modes pertaining to n = 1 are taken for solution. As the fundamental natural frequency was 22.37 rad/sec for n = 1, a time step of about l/50 set is taken for response determination.
Dynamic analysis of a pressure vessel
169
Fig. 12. Cantilever model (dimensions in mm).
For hysteritic damping in the material, an l/r value of 300 (for steel) [12] is taken. An arbitrary value of 0.5 for the damping ratio is set as the upper limit in the fifth mode and Rayleigh damping curve fit is used to interpolate at the other natural frequencies.
8. RESULTS
AND DISCUSSIONS
8.1 Static analysis For a uniform internal pressure of 220 atmospheres, the static stresses are determined. The stresses as obtained by the two different approaches at the junction node are compared with those obtained by the finite element method (using the shape function which is a linear polynomial in radial and axial directions) in Figs 7 and 8. Considering the pressure vessel as a closed ended thin cylinder subjected to internal pressure, the stress values predicted by elementary strength of materials are calculated and are shown in Figs 7 and 8. The ratio of circumferential to meridional stress in the main cylindrical body is obtained as 2.13 whereas strength of material predicts a value of 2.0.
-10. 0
8
10
Time in set
I 12
6 4 Time in set
2
la)
8.2 Dynamic analysis The eigen-pairs are determined for three values of n (1, 2 and 3). The first four mode shapes for each value of n are shown in Figs 9-l 1. A crude approximation of the pressure vessel as a cantilever beam is shown in Fig. 12. This is used to validate the results for n = 1. Since u displacement is very small compared to I) and w (ratio varying from & to h) and the latter two are of the same order (because n = I), only w displacements need be considered for potential and kinetic energy terms. Total mass of the structure, M = 360 x 10’ kg Mass of the lid, m = 28.5 x 10’ kg Length (span) of the cantilever, L = 23.0 m
I
-7.51 (bl
10
Fig.
13.
Relative
I 11
acceleration
excnation.
I
responx+horizontal
V. RAMAMUIUIand P. SESHU
-501
I
0
I
I
2
L
I
I
I
Time
I
I
I
6
6
10
I
I
12
in ret
Fig. 14. Dynamic stress plot-horizontal
Second moment of area, Z = 2.75 m4 Young’s modulus for steel = 2.058 x 10” N/m*
It is to be noted that the span of the equivalent cantilever beam (determined only very approximately) can strongly influence the value of the natural frequency as obtained above. In the mode shapes it is observed that a node for axial displacement (a) corresponds to an anti-node for radial displacement (w) and vice versa. Also for n >O,uisobservedtobeoftheorderof(-w/n).The participation of the lid is minimal for II = 3. 8.3 Seismic response The relative acceleration response of the structure is shown in Fig. 13. In Fig. 13(b), the time scale is chosen so as to indicate the time period of the response. When the ground acceleration is given in the horizontal plane i.e. D and w directions of the structure, a peak response of 8.8 g (i.e. a dynamic magnification factor of 14.7) is obtained at node 25 on the cylindrical portion, 800 mm from the flat head. The dynamic stresses in the material of the structure are also obtained. Node 232,500 mm away from the support on the skirt, is found to have the maximum stress. The dynamic stress plot for this node is shown in Fig. 14. In the main cylindrical body of the structure a peak stress of about f 12.20 kgf/mm* is obtained. 9. CONCLUSlONS
The dynamic response of a pressure vessel to a typical earthquake has been obtained. A modification
excitation.
to the existing method of treatment of a junction node is presented. It has been observed that the skirt portion of the structure is the maximum stress zone. The stress versus time plot over a period of 12.5 set for a node on the skirt is presented. REFERENCES
1. S. A. Anagnostopoulos, Wave and earthquake response of offshore structures; evaluation of modal solutions, AXE, J. Struct. Div. loS, 2175-2191 (1982). 2. J. W. Tedesco, C. N. Kostem and A. Kahtins, Rapid calculation of mode participation factors for circular cylindrical shells. Comput. Sfruct. Zo, 509-515 (1985). 3. J. Altes, H. Grafii and D. Koschmieder, Investigation of the earthquake behaviour of the research reactor FRJ-2 (Dido). Report No. JUEL-SPEZ-222 (1983). 4. B. Budiansky and P. P. Radkowski, Numerical analysis of unsymmetric bending of shells of revolution. AZAA Jnl 1, 1833-1842 (1963). 5. V. Ramamurti and R. S. Alwar, Stress analysis of tube mills. J. Strain Anal. 8, 200-208 (1973). 6. M. L. Potters, A matrix method for the solution of a second order difference equation in two variables. Mathematische Centrum, Amsterdam, Holland, Report No. MR 19 (1955). 7. Alan Jennings, Matrix Computation for Engineers and Scientists. John Wiley, London (1977). 8. V. Ramamurti and 0. Mahrenholta, The application of the simultaneous iteration method to flexural vibration problems. Int. J. Mech. Sci. 16, 269-283 (1974). 9. J. Krishna and A. R. Chandrasekharan, Elements of Earthquake Engineering. Sarita Prakashan, Meerut (1976). 10. K. J. Bathe and E. L. Wilson, Numerical Methods in Finite Element Analysis. Prentice-Hall, New Delhi (1978). 11. E. L. Wilson, Earthquake analysis of reactor structures. In Symposium on Seismic Analysti of Pressure Vessel and Piping Components, pp. 53-73. ASME, New York (1971). 12. W. K. Wilson, Practical Solutions of Torsional Vibration Problems, Vol. 2. Chapman & Hall, London (1963).