Longitudinal vibration analysis of partially-filled ellipsoidal tanks

Longitudinal vibration analysis of partially-filled ellipsoidal tanks

Compurers& Structures, Voi. 3, pp. 205-215. PcrgamcmPress 1973.printedin Great Britain LONGITUDINAL VIBRATION ANALYSIS OF PARTIALLY-FILLED ELLIPSOIDA...

697KB Sizes 4 Downloads 76 Views

Compurers& Structures, Voi. 3, pp. 205-215. PcrgamcmPress 1973.printedin Great Britain



Martin Marietta Corporation,

Baltimore, Maryland, U.S.A.

and T. JAMES RUDD§ Bellcomm, Inc., Washington,

D.C., U.S.A.

Abstract-To evaluate the POGO Stability of large liquid fueled rockets, such as the SATURN V, the Iongitudinal vibration characteristics must be accurately known. In particular this moans the inertial properties of the liquid propellant must be properly represented. Rather than just the analysis of liquid motion in an elastic container this requires derivation of a liquid mass matrix which is compatible with the mass matrix used in an overall structural analysis of the vehicle. In this paper a new method of finding this liquid mass matrix based on a finite-difference approximation to Laplace’s equation is described. The method is used to find liquid mass matrices for the liquid hydrogen and liquid oxygen tanks of the SATURN S-IV B stage at various flight times. The matrices are incorporated into a structural model of the S-IV- B and a vibration analysis carried out. The results of the vibration analysis are compared to the spectrogram traces from the Apollo 12 mission.

IN THE past, two approaches have been used to derive liquid mass matrices. The first approach [l, 21, which is restricted to axisymmetric propellant-tank systems, involves discretizing the fluid region into a series of cylindrical or toroidal rings and using an assumed velocity field to derive the liquid mass matrix. In general this velocity field will not be irrotational. Since the motion of an inviscid propellant which starts from a rest, must be irrotational these methods [1, 2] do not yield very reliable results. The second approach [3, 43, the so called hydroelastic method, is based on the assumption that the velocity field is irrotational and derivable from a potential function. The kinetic energy of the liquid is then expressable as an integral over the boundary of the liquid whose integrand involves the velocity potential and its normal derivative. After evaluating t:?is integral using a finite number of points and then utilizing a relationship between the velocity potential and its normal derivative, it is a simple matter to come up with the liquid mass matrix. The central task, then, in the second approach [3,4] is to derive the relationship between the velocity potential and its normai derivative at a given set of points on the liquid boundary. In 131,Guyan et al. use a simple source distribution on the boundary to represent the liquid behaviour. The velocity potential and normal velocity are then solved in terms of the source density and hence in terms of each other. Khabbaz [4] uses a similar procedure

tpresented at the National Symposium on Computerized Structural Analysis and Design at the School of Engineering and Applied Sciences, George Washington University, Washington D. C., 27-29. March 1971. &Senior Research Scientist, Research Institute for Advanced Studies. $Presently Staff Scientist, Ensco, Inc., Springfield, Virginia. 205



to Guyan et al. [3] but whereas [3] is restricted to axisymmetric systems, [4] can be used for the general three dimensional problem. In the work presented here the artifice of a source distribution is not required. Indeed since the liquid is assumed incompressible it follows from the continuity equation that the potential function satisfies Laplace’s equation. Using a finite-difference solution to this equation (in which the grid lines may be chosen to intersect the liquid boundary at preassigned points) it is a straight forward task to find the desired relationship’between the velocity potential and the normal velocity on the liquid boundary. The observation that for hydroelastic vibrations the potential function may be taken as identically zero on the free surface, further simplifies the solution. The method may be applied to any liquid tank system, axisymmetric or not, provided the liquid is assumed incompressible and its motion irrotational. The tank geometry can also have discontinuities in its slope (such as in the S-IV LH, tank). THEORETICAL


The problem presented here is to derive a mass matrix by finite differences for a liquid contained in a flexible axisymmetric tank. The mass matrix in this instance, is to be compatible with the mass and stiffness matrices generally used in vibration analysis. Although the solution obtained is applicable to an arbitrary axisymmetric tank, we will restrict our discussion to longitudinal and radial motion of an elliptical tank whose cross-section is as shown in Fig. 1. This tank, which has the same geometry as the SATURN S-IVB liquid FREE SURFAC






Fro. 1. Elliptical tank system showing boundary points and degrees of freedom. oxygen tank, consists of two ellipsoidal bulkheads joined along their maximum openings. The origin of the ellipsoidal bulkhead is at 01 while for the upper ellipsoidal bulkhead the origin is at 02. We assume a cylindrical co-ordinate system (r, 8, Z) with origin at 01. Since the problem is axisymmetric (the Z axis being the axis of symmetry) none of the variables have any 0 dependence. The liquid is assumed to be incompressible and inviscid, these being valid assumptions for the liquid propellants considered. The free surface of the liquid can be anywhere between the lowest and highest point of the tank. The space between the free surface and the top of the tank is occupied by ullage gas whose dynamics in this derivation have been ignored.


Vibration Analysis of Partially-Filled

Ellipsoidal Tanks


The size of the mass matrix for the liquid depends on the number of degrees-of-freedom chosen to represent its motion. As will be shown the motion of the fluid is completely determined once the motion of its boundary is known. Hence the number of degrees-offreedom for the liquid (with the exception of the liquid free surface) can be chosen to be the same as those used to represent the motion of the tank structure. These degrees-offreedom are simply the vertical displacements, Z,, and radial displacements, R,, at the selected tank boundary points, Pi, shown in Fig. 1. While these boundary points will usually be chosen to make the stiffness matrix calculation for the tank as simple as possible, the finite-difference solution requires that corresponding points on the lower and upper bulkhead lie on the same vertical lines (see Fig. 1). Derivation of liquid mass matrix

Because of the connection between the mass and kinetic energy of a liquid, once a discrete form of the kinetic energy is known, a mass matrix can be formed. The present procedure for determining the mass matrix is specifically aimed at evaluating the kinetic energy of the liquid solely in terms of the tank wall motion. Basically the kinetic energy of a contained liquid is given by the expression:

where p is the liquid mass density v=v(r,Z,t) is the axisymmetric fluid velocity vector field, and the integration is carried out over the volume V occupied by the liquid. If we consider only small perturbations of the fluid, which we assume is inviscid and initially at rest, then the fluid motion remains ii-rotational and a velocity potential Q(r,z’,t) can be defined such that its gradient yields the fluid velocity vector, v



from the continuity equation for an incompressible fluid, div v=O


and equation (2) one obtains for the fluid region the linear Laplace relationship

[email protected]=0.


From Gauss’s theorem we can rewrite equation (1) solely in terms of the total bounding surface S enclosing the fluid so that the kinetic energy



@!%js an


where &D/an is the velocity of the fluid in the direction of the outward normal to the liquid surface and ds is the incremental surface area on the boundary S. The integration is extended over the total bounding surface, S, of the liquid which includes both the free surface, S,, and the wetted area along the tank wall, S,.




At this stage we make an important simplification in (5) through the following consideration of the free surface boundary condition. First of all on the free surface the pressure of the liquid equals the ullage pressure. During the vibration of the tank, ullage pressure perturbations are usually so small that the ullage pressure can be assumed to be constant. The linearized Bernoulli equation for the free surface then has the form

[email protected] at+a’q=O


where a, is the magnitude of the longitudinal inertial acceleration vector of the rocket relative to the local gravitational acceleration and r(r,t) the instantaneous wave height of the free surface. As Goldman [5] and others have noted the second term in (6) contributes solely to the potential energy of the free surface. Its influence on vibration modes that are essentially associated with tank wall motions is small and as a result may be omitted. In the present derivation equation (6) is further reduced to

[email protected], --



on the free surface. For oscillating motion this is equivalent to saying that on the free surface @=O. (8) Since Q is zero on the free surface, we need only evaluate equation (5) along the tank wall. Thus the kinetic energy is given by the simpler expression



2s ss an


where now the integration is taken solely over the wetted area, S,, of the actual tank wall. The integral in (9) can be approximated by using a banded symmetric integration matrix [G] to tie together the values of the integrand at any N discrete points on the wall. Thus the fluid kinetic energy in matrix notation becomes


(L =l, N)


where DL and

am 0 an, are now the specific values of the velocity potential and outward normal velocity of the liquid at the L discrete boundary points. The banded symmetric matrix [G] can be found from equation (9) by introducing a set of interpolation functions such that the velocity potential and outward normal velocity along the tank wall vary linearly with arc length between any two consecutive boundary points P, and PL+l (see Fig. 1). (For more details of the calculation of the [G] matrix see Ref. [5]). Since the liquid is inviscid, the compatibility condition between the liquid and tank is


Vibration Analysis of Partially-Filled


0;; -


Ellipsoidal Tanks





where uL is the normal outward velocity of the tank structure at boundary point L. From a finite difference solution of (4), subject to the boundary conditions (8) and (1 l), we can find (as will be shown below) a matrix [Q] whose elements result in the interconnection relationship (12)


Substituting the results of equations (11) and (12) into equation (10) transposed gives the discrete form for the kinetic energy solely in terms of the elements of [c?] and [Q] and the normal outward velocity of the tank wall.

T=${uL}'[G][Q]{u&=1, N


Finally from Lagrange’s equation (the governing equations fro the liquid motion) we identify the inertial properties of the field as being contained in the Lagrangian expression (14) where [M] is now the liquid mass matrix of the fluid. Substituting the kinetic energy expression for T from equation (13) into (14), we obtain directly the desired fluid mass matrix



The central problem in the derivation of the mass matrix, equation (15), by finitedifferences is the formulation of the [Q] matrix. As indicated, this is found by constructing a solution for Laplace’s equation (4) subject to the hydroelastic boundary conditions (8) and (11). Since the motion is assumed axisymmetric we need only consider the solution of @ in the half section shown in Fig. 1. From the assigned bounday points the grid is constructed as shown in Fig. 1. In cylindrical coordinates, Laplace’s equation (4) has the form

For convenience we introduce the non-dimensional variables I=r/R and z=ZIR, is the maximum radius of the tank. In non-dimensional form (16) becomes






Approximating the terms in (17) by central differences we find for the general interior point (i, j), Fig. 2, the equation

where Qpi,j is the approximation at (i- 1, j), and so on.

to @ at the point (i, j), @i- r,j is the approximation

to @




Bi &j-l)





(i, j+l)




FICA 2. General interior point.

For points on the axis we note that Laplace’s equation has the special form [email protected]~~--CDZ,=o.


The corresponding finite difference approximation for an axis point (j= 1) is now -1 -1 :tQi, ~~-~(~~-~+~i)~i~l’l+pI(~~-~+~~)~i+”l l=O* C20) Ul 2



To complete the finite difference solution we must add equations for the boundary points. Since the potential at the free surface is zero we need only write equations for the boundary points in contact with the tank wall. At these points (taking account of the normalization)

aQr 0 an.


1, IV)


where as previously stated II~ is the normal velocity of the tank at the boundary point L. The next step is to suitably express

[email protected] z-l 0


in terms of the potential at the boundary and interior points. Here a first order approximation to the derivative

[email protected] an. 0

Longitud~al vibration Analysis of P~ti~ly-FilIed

is sufficient since the approximations the expression for

Etlipsoidal Tanks


to the interior points are only to this order [6]. Since

adt, an, 0 varies with the geometry of the tank, the equations for the boundary points on the upper bulkhead are considered separately from those on the lower bulkhead. Boundary points on upper and lower bulkheads Referring to Fig. 3 we see that on the upper bulkhead the normal at the point PL makes an angle of QL with the horizontal. The mesh points adjacent to P, in the horizontal and vertical directions are respectively. Pn andPL-r.

FIG. 3.

Boundary point on upper bulkhead.

Now resolving along the normal gives

Using first order approximations

(23) Substituting equations (23) into (22) we find for the upper boundary point that sin CILQ CosszLmL+ aL




SimiIarly for the boundary point on the lower bulkhead, shown in Fig. 4, we find that



FIG. 4. Boundary

point on lower bulkhead.

solution .for the Q metric From the finite-difference equations for the interior points (equation (IS)), axis points (equation (20)) and boundary points (equations (24 and 25)) we derive the matrix relationship

CA1 iar*,j>= RIvLj>


where [A] is the coefficient matrix obtained from the finite difference equations, {Bi, jj is a column matrix of velocity potentials at every grid point (interior, axis and boundary) and {vi, j} is a column vector of normal velocities such that vi, j=O when point i, j is not on the boundary vi, j=uL when point i, j coincides with boundary point L.


From equation (26) we find that

By judicious partitioning of R[A]-’ we can pick out those values of @;,j and vi, j that lie on the boundary and establish the relationship

The resulting [QJ matrix is exactly the interconnection relationship specified previously in equation (12) and leads directly to the desired mass matrix formulation in equation (15). Computer requirements Since in all hydroelastic models a liquid mass matrix has to be obtained using a digital computer the computer storage requirements and run times are an important consideration in the development of such a model. One of the attractive features of the model presented here is that the computer requirements are relatively low. Indeed for the special but not trivial case of the elliptical tank shown in Fig. I, the program run time is only fractionally more than that required to invert an Nx N matrix, where N is the number of designated mass or boundary points. Storage requirements are equally low being in this case equal to 3N2 locations.


Vibration Analysis of Partially-Filled

Ellipsoidal Tanks


An important factor in reducing the computer requirements stems from the fact that no degrees of freedom are required for the free surface. The only effect of this approximation is that the model does not predict the pure longitudinal sloshing modes, however, the frequencies of those modes are well below those of interest for POGO analysis. For a tank with more complicated geometry, such as the liquid hydrogen tank, shown in Fig. 6, the computer requirements are more than for the elliptical tank of Fig. 1. The reason for this is that in general to find the [Q] matrix one has to invert the [A] matrix (see equation 28) whose dimensions are NTx NT, where NT is the total number of grid points, this number is of course greater than N the number of boundary points. It turns out in the case of the elliptical tank that the topology of the grid net allows one to find the [Q] matrix without in fact inverting the [A] matrix (see Vandergraft [6]). APPLICATION The method described here for generating liquid mass matrices has been used extensively in the POGO related structural dynamics work done at Bellcomm on the SATURN V. Results similar to those reported by Ujihara and Guyan [7] were previously obtained for the S-II LOX tank. In this section we describe the more recent results of a finite-difference hydroelastic analysis of the SATURN S-IV B stage. A schematic diagram of the lumped parameter model of the S-IVB is shown in Fig. 5. The model has 51 masspoints with a total of 88 degrees-of-freedom. (All points off the THRUST STRUCTURE


FIG. 5. Longitudinal


model for Saturn S-IV B stage showing mass points.

axis have both radial and vertical degrees-of-freedom.) Considerable detail was included for those parts comprising a large proportion of the total mass (e.g. the liquid propellants, the engine and the thurst structure) while those parts less massive and more removed from the engine were more simply modeled. The details of the LOX tank have been given already in Fig. 1 and the derivation of the liquid mass matrix discussed. The LH2 tank is shown in Fig. 6 with the corresponding boundary points (mass points) and finite-difference grid. The derivation of the liquid mass matrix is completely analogous to that of the LOX, although the boundary is slightly more complicated to deal with. The total mass matrix for the S-IV B stage, at any flight time (note that the mass matrix is time variable due to the consumption of propellant) is found by adding the contribution from the LOX, LH2 and the structural mass of the stage. The structural mass is found by a simple lumping technique.














FIG. 6. Axisymmetric cylindrical tank system showing finite difference grid, boundary points and associated degrees of freedom.

The stiffness properties of the stage were found by use of the axisymmetric shell program of Kaufman [S]. For the part forward of the LH, tank a simple beam-rod analysis was used to estimate the stiffness. The natural frequencies of the stage were found from the equation C&M] - [K]=O


where [M] is the total mass matrix at any time, t [lu] is the total stiffness matrix at time, t w is the circular frequency. The modal frequencies are shown as a fraction of time in Fig. 7. The modes that are primarily connected with the liquid propellant are those for which N=3,4,5 and6. The first mode (N=l) is a low gain mode connected to the matrix motion of the lunar module while the second mode N=2 is the first longitudinal compression mode. As a check on the analysis, flight data frequencies obtained from Bellcomm spectrograms of Apollo 12 have been superimposed on the graph. It can be seen that the theoretical model closely predicts all the major modes that show up in the flight data.


Vibration Analysis of Partially-Filled


FIG. 7. S-IVB/Apollo

Ellipsoidal Tanks



12 hydroelastic model frequencies compared with flight data.

REFERENCES PI J. S. ARCHERand D. P. RUBIN, Improved analytical longitudinal response analysis for axisymmetric launch vehicles. I, NASA CR-345 (Dec. 1965).

PI P. J. GRIMES,L. D. MCTIGUE, G. F. RILEY and D. 1. TILDEN, Advancements

in Structural Dynamic Technology Resulting from Saturn V Programs, II, NASA CR-1540 (June 1970). [31 R. J. GUYAN, B. H. UJIHARAand P. W. WELCH, Hydroelastic analysis of axisymmetric systems by finite element method, Proc. 2nd Conf: Matrix Methods in Structural Mechanics, AFFDL-TR-68-150 (Dec. 1969). [41 G. R. KHABBAZ,Dynamic Behavior of Liquids in Elastic Tanks, LMSC 60-80-70-23, (Aug., 1970), Lockheed Missiles and Space Co., Palo Alto, California, U.S.A. PI R. L. GOLDMAN,Longitudinal Vibration Analysis of Partially-Filled Ellipsoidal Tanks by Finite Differences, TR 70-6C (Aug., 1970), Research Institute for Advanced Studies, Baltimore, Md., U.S.A. A Finite Difference Method for Solving Laplace’s Equation, Using a Non-Uniform [cl J. S. VANDERGRAFT, Mesh, TM-71-2031-1 (March, 1971), Bellcomm, Inc., Washington, D.C., U.S.A. Vl B. H. UJIHARAand R. J. GUYAN, Hydroelastic Properties of a Full Scale S-II LOX Tank, AZAA Paper Paper No. 72-173, AIAA 10th Aerospace Sciences Meeting(Jan. 1972). PI S. KAUFMAN,Axisymmetric Shells, TM-69-2031-3, (July 1969)., Bellcomm, Inc., Washington, D.C., U.S.A. (Received 22 February 1972)