Two-dimensional unsteady flow in propellant tank under low gravity by finite difference methods - reduction to boundary value proble

Young, C. T. K.

English

NASA Technical Reports Server

TECHNICAL NOTE R-164 
1, 
0 w c ANALYSIS OF TWO-DIMENSIONAL, 
UNSTEADY FLOW IN A PROPELLANT 
TANK UNDER LOW GRAVITY BY 
FINITE DIFFERENCE METHODS 
October 1965 
Hard copy (HC) 
Microfiche ( M  F) 
Y 6 5 3  Julb 65 
RESEARCH LABORATORIES 
BROWN ENGINEERING COMPANY,  INC. 
H U N T S  VIL L E,  ALABAMA 
https://ntrs.nasa.gov/search.jsp?R=19660004273 2020-03-16T21:29:12+00:00Z
+I. 
TECHNICAL NOTE R-164 
ANALYSIS OF TWO-DIMENSIONAL, UNSTEADY 
FLOW IN A PROPELLANT TANK UNDER LOW 
GRAVITY B Y  FINITE DIFFERENCE METHODS 
October 1965 
P repa red  For 
PROPULSION DIVISION 
P & V E  LABORATORY 
GEORGE C. MARSHALL SPACE FLIGHT CENTER 
RESEARCH LAB ORATORIES 
BROWN ENGINEERING COMPANY, INC. 
Contract No. NAS8-20073 
Prepared By 
C a r l  T .  K. Young 
.. 
ABSTRACT 
This report  descr ibes  a numerical procedure for solving a two- 
dimensional, unsteady flow problem. The fluid is in a tank, has a f r e e  
surface,  a periodic source (to be terminated) on the side of the tank, and 
is subject to a near ze ro  gravitational field. 
incompressible and i r rotat ional  so tha t  the problem reduces to a boundary 
The flow is assumed to  be 
value problem governed by the Laplace equation. 
Approved 
Mechanics and Thermodynamics 
Department 
Approved 
Director of Research 
ii 
TABLE OF CONTENTS 
INTRODUCTION 
NUMERICAL ANALYSIS 
Problem De sc riptions 
Finite - Difference Method 
Calculation of Flow Field 
Calculation Procedures  
c ONC LVSIONS 
REFERENCES 
iii 
Page 
1 
3 
3 
4 
10  
14 
15 
16 
b 
LIST OF SYMBOLS 
H 
h0 
T 
t 
U 
UO 
V 
W 
X 
Y 
A 
rl 
e 
Mesh size,  f t  
Depth of slot ,  f t  
Fraction of mesh  s ize  
Effective gravity, f t / s ec2  
Liquid height in static equilibrium, f t  
Average liquid height under static equilibrium condition, ft 
Surface tension, lbf/ft 
Time, sec  
x-component velocity, f t /  sec 
Source velocity, f t /  sec  
y-component velocity, f t /  sec 
Width of the tank, ft 
x-coordinate 
y-coordinate 
Diff e r enc e 
Perturbation of liquid height, ft 
Contact angle, degree 
Density, s1ug/ft3 
Velocity potential, f t 2 /  sec  
Re laxation factor 
i v  
LIST O F  SYMBOLS (Continued) 
Superscripts 
k Number of iterations 
S Point on the f r ee  surface 
Subscripts 
i Index for x-coordinate 
j Index for y-coordinate 
n Index for t-coordinate 
V 
. 
INTRODUCTION 
The study of dynamic behavior of liquids in greatly reduced 
gravitational fields has  become a subject of intense in te res t  in  recent 
years  a s  space technology has progressed. 
of the oscil latory motion of fuel in  a large rocket propellant tank can 
help the designer to prevent its being in resonance with the motion of 
the vehicle or  of its control system, thus avoiding such undesirable 
effects a s  dynamic instability. 
face is important to  the engine res ta r t  operation af ter  a prolonged 
coasting period, a s  the liquid and vapor tend to mix together when the 
gravity force is almost  absent. Other applications can be found in the 
life support system, fuel ce l l  system, e t c . ,  in a space vehicle. 
Knowledge of the frequency 
The behavior of the liquid-vapor inter-  
The present  study is concerned with the dynamic behavior of 
the f r ee  surface,  and the flow field in general ,  of liquid hydrogen in 
the propellant tank of the S-IVB rocket stage during flight after the 
main fuel supply has  been cut off. A water hammer  effect is created 
by such a sudden cutoff operation and becomes the periodic source of 
the tank. 
reduced to a magnitude of about l o m 5  g and surface tension can no 
longer be neglected when considering the behavior of the f r ee  surface.  
For  such a problem one is  required to solve the ful l  Navier-Stokes 
equations in two dimensions with a moving boundary, the so-called 
Stefan's problem. Aside from the additional complications ar is ing 
f r o m  the presence of a moving boundary, the analytical o r  numerical  
solutions of Navier-Stokes equations present  tremendous difficulties. 
Most of the attempts in the past  have failed in the obtaining of a numerical  
solution to the ful l  Navier-Stokes equations. 
with simple boundary conditions has limited success  been achieved. 
During this stage of flight the gravitational field has been 
Only in  a few problems 
Therefore ,  in order  to obtain a 
problem, some simplifications 
more pract ical  solution in the present  
must be made. 
1 
. 
be neg 
dition 
The f i r s t  simplification is to assume that the viscous effects can 
ected and that the flow i s  irrotational. 
s used with the continuity equation to  obtain Laplace 's  equation in 
Satterlee and Reynolds' have found that 
The i r rotat ional i ty  con- 
t e r m s  of the velocity potential. 
in the case  of a c i rcu lar  cylindrical  tank the inviscid analysis  gives a 
na tura l  frequency of oscil lation only a few percent  lower than actual. 
The assumption of potential flow i s  therefore  useful  and pract ical ,  and 
the re  is no reason  to  believe that it cannot be  applied to  the two-dimensional 
case .  
Next, the concept of s m a l l  perturbations is introduced, i. e . ,  i t  is 
assumed that the f r e e  su r face  can  be expressed  as the s u m  of the static 
equilibrium sur face  and a s m a l l  perturbation. 
simplifying the f r e e  surface conditions and fixing the mesh  points on the 
f r e e  sur face  when finite-difference techniques a r e  employed. Also, fluid 
proper t ies ,  including the density and surface tension, a r e  considered 
c ons tant . 
This has the aclvzctage c?f 
Even with such simplifications, an analytical  solution of the 
present  problem i s  not i n  sight and numerical  solutions must  be attempted. 
The following chapter outlines the procedure for  the solution of this prob- 
lem by finite - difference methods. 
2 
NUMERICAL ANALYSIS 
dH d2H 3 -  - a? - -  dx dx2 -
ax 
PROBLEM DESCRIPTIONS 
The assumption of an  incompressible,  potential flow leads to  the 
following boundary value problem. 
= o  (1 )  , 
Governing. Differential Eauation 
a+a = o  in the region R 
a x 2  a y 2  
Boundarv Conditions 
On Solid W a i i  
3 = o  n normal to the wall an 
At the Source 
3 = uo(t) ax 
Dynamic Surface Condition 
Kinematic Surface Condition 
3 
Contact Angle Condition (No Hysteresis)  
Initial Conditions 
utx, y, 0 )  = 0 
The equations for the surface conditions a r e  derived in  Refer- 
ence 2 .  
the source velocity uo(t) is a periodic function. 
The fluid is assumed to be ir, stztic equilibrium initially, and 
FINITE-DIFFERENCE METHOD 
The continuous derivatives 
in the differential equations above 
a r e  replaced by discrete  finite 
differences as follows: 
N W  - 
W -
s w  - 
't- 
hS 
1 
I 
4 
NE - 
E - 
SE - 
F o r  equal m e s h  s izes ,  the truncation e r r o r s  in  each of these 
formulas  would be O(h2) a s  h+O, as compared to  O(h) for  different mesh  
s i zes .  This shows the g rea t  advantage of using equally spaced nets.  
A network of m e s h  points is laid over the flow field as shown in  
F igure  1. 
e r r o r s ,  s torage capacity, running t ime of the computer ,  and programming 
ease .  The mesh  s ize  must  be chosen such that both the width of the tank 
and the depth of the slot  (h, - h l )  a r e  integral  multiples of the mesh  s ize .  
Ordinar i ly  this i s  impossible without making the mesh  s ize  intolerably 
sma l l ;  however, with slight adjustment of e i ther  o r  both of W and d, the 
m e s h  s ize  may be chosen with relative ease.  
should a l so  be adjusted, if necessary,  to meet  cer ta in  programming 
requirements .  
without appreciably affecting the flow charac te r i s t ics .  
then the same throughout most  of the flow field. 
su r f ace  mesh  points a r e  selected a t  the intersect ions of the static equi- 
l ib r ium surface and the ver t ica l  l ines,  therefore  the mesh  s izes  a r e  
different.  
Choice of mesh  s ize  a depends on considerations of truncation 
The position of the slot  
F o r  all pract ical  purposes these adjustments can be made 
The mesh  s ize  is 
Adjacent to  the f r ee  
5 
6 
At each of the regular  interior mesh  points the Laplace equation 
is replaced by  a finite-difference equation using the five point formula.  
1 ( 3 )  +i, j , n t  1 = 4 (+i- 1, j , nt  1 t- +i, j - 1, n+ 1 + +it 1, j nt 1 ' +i, j t  1, n+ 1 1 
F o r  mesh  points on the solid wall  the Neumann boundary condition 
a+/ an = 0 can be satisfied by locating dummy points outside the wall  at the 
images  of the points immediately inside the wall, so  that 
On the left wall, 
1 
9 1, j ,  n t  1 = 4 (2  $2, j n t  1 + + 1, j -  1, n t  1 + + 1, jt 1, n t  1) 
On the right wall, 
1 - -  
+I, j ,  n t  1 - 4 ( 2  +I- 1, j ,  n t  1 +I, j -  1, n t  1 + +I, j t  1, n t  1) 
On the bottom of the tank, 
1 
) +i, 1 ,  n t  1 - 4 ( 2  +i, 2 n t  1 + +i - I, 1 nt 1 + +it 1,1, n t  1 - -  
At the two co rne r s ,  
At the source ,  since ?k = u0(t) ,  ax 
(4)  
( 5 )  
7 
On the f r ee  sur face ,  
dH d2H 3 T -  -
dx dx2 - 
5/2 
L 
avs 
+ vi, j ,  n + (F)~, j ,  q i , n  
S S 
where velocities on the surface U i ,  j ,  
procedure to  be descr ibed later.  (Actually they a r e  velocities a t  points 
with coordinates corresponding to  the s ta t ic  equilibrium sur face .  ) For 
and vi, j, can  be obtained b y  a 
avs (x)~, j ,  a backward difference is used. 
S S 
- vi, j ,  n - vi, j ,  n-i - 
(%)i, j ,  n A t 1  
8 
a r e  obtained by modifying the appropriate 
d2H and -dH The t e r m s  - 
dx dx2 
equations of Reference 3. They a r e  
1 
2 
-
1 
- 1  dH - = f  
t 2 c o s 8  -  W 
dx H - H o  
W 
where 
Ho = H when x = W/2 
Bo = p g WZ/T = Bond number. 
is 
dH depends on the value of x. F o r  x > W / 2 ,  - dH The sign for  - 
dx dx 
positive; - is negative for x < W / 2 .  dH 
dx 
t 2 cos e -  
W W - -  - d 2H 
3 
t 2 c o s e - -  - 
2 
dx2 
L 
where 
W 
ho = $ lo H dx . 
Also,  
9 
S 
Thus, +i ,  j ,  n t  1 can be calculated by Equation 10. 
t e r m s  of a Taylor s e r i e s  expansion, 
Using the first two 
1 S ’i, J- 1, n 
t AY 2 
Vi ,  j ,  n t 2  
where 
(F) can be  calculated s imilar ly .  
CALCULATION OF FLOW FIELD 
In a typical ca se  there  a r e  at l eas t  severa l  hundred mesh  points 
in the flow field and there  a r e  an equal number of equations to  be solved 
simultaneously at any time level. Since the resulting matrix is spa r se ,  
the bes t  way to solve these equations is by an i terative method. 
i terat ive scheme to be used is the successive over-relaxation method 
o r  the extrapolated Liebmann method when applied to Laplace’s equation, 
which has the fo rm 
The 
k t l  - -  w k  k t l  t + i , j t i , n t l t ~ i , j - i , n t i ) - ( w - l ) + i , j , n t l  k k t l  k 
+i ,  j , nt 1 - 4 (+it 1, j ,  n t  1 t +i- 1, j ,  n t  1 
10 
The relaxation factor,  w,  lies somewhere between 1 and 2 .  The 
optimum value of w,  which gives the highest convergence rate ,  is given 
by Young's4 formula 
where A is the spec t ra l  radius of the Gauss-Seidel method. 
of Laplace's  equation and for  a rectangle with s ides  Ra and Sa, where a 
F o r  the case  
is the mesh  s ize  and R and S 
1 A = -  
4 
a r e  integers. X is given by 
(cos g Tr t cos q 
s -  
Since the configuration under consideration is not rectangular,  it 
is suggested that one use either a rectangle which contains the given 
region o r  one which has approximately the same a r e a  and proportions. 
These formulas  have been found to ag ree  ve ry  well with the resu l t s  of 
tes t  runs,  re fe r red  to la te r ,  for  which the liquid surface was assumed 
to  be flat. 
The successive over-relaxation method has been proven to be 
superior  over the Gauss-Siedel method in t e r m s  of convergence ra te .  
This is so when applied to the five point formula,  provided that the 
ordering of the points is properly chosen. Experience shows that by 
i terating column wise f rom left to right with a r rows  pointing upward 
in each column, the highest convergence ra te  can be obtained. 
It should be pointed out, however, that in obtaining mesh  points 
immediately below the f r ee  surface,  i t  is usually undesirable to apply 
the modified five-point formula since the mesh  s ize  ra t ios  in such cases  
m a y  be s o  large that the iterations may  converge very  slowly and some- 
t imes  may even diverge. 
formula is the so-called "interpolation of degree one": 
A different method is needed. The following 
11 
. 
.', o r  (xi ,  j - a ,  y i , j )  is a point outside the surface.  If one of them is such a 
point, the value of + at the surface point intersected by the horizontal 
line y = ( j - 1 )  a i s  obtained through linear interpolation of the two nearest  
surface points. 
of (3) ax i ,  j ,  n t i  
I 
I 
I 
This new surface point is then used for  the calculation 
. For  example, as shown below, it is intended to find uc. 
S 
+i, j t  1, 
+i, j ,  n t  
S 
+ i , j - i , n t i t  + i , j t i , n t I  
(12) 
- 
+i, j ,  n t l  - l t f  
+i, j -  1, n t  1 E 
The truncation e r r o r s  in such interpolations a r e  O(aL), compared 
with O(a)  in the modified five-point formula and O(a2) in  the five-point 
formula.  This shows the advantage in using the interpolation formula. 
Note a l so  that Equation 1 2  is of positive type with diagonal dominance, 
which is excellent for  numerical  solutions. 
After the velocity potentials, + I s ,  have been calculated for 
every  mesh point on the new time level, the velocities will be computed. 
The following procedure is to be  followed: 
vi, j ,  n t  1 = (F) Y i, j , n t l  
To find (*) , it is necessary to  tes t  whether ( X i ,  j t a, yi,  j )  
ax i ,  j , n + l  
12 
Point L, with coordinates (xc - a, yc) ,  is located outside the surface.  
Thus,  
L dl' J' 'S -/A. . t 1 , n t i  E 
I 
I 
Then 
Also , 
F o r  velocities on the f ree  surface,  apply a Taylor s e r i e s  
expansion to f i r s t  degree 
1 A Y l +  Ay2 A Yz2 ui, j - l , n t l  A y l  ( A y l  t ui, j-2, n t l  - S j ,  n i l  - 
Similarly, 
CALCULATION PROCEDURES 
Before proceeding to  calculations fo r  points on a new t ime level, 
the t ime increment must  be selected. 
on the stability and convergence requirements,  and is  mostly determined 
by  experience. Surface conditions a r e  then obtained by Equations 10 and 
11, and values of 4 in  the flow field a r e  approximated by extrapolations 
through the corresponding points on the two previous t ime levels. 
the i teration process  s ta r t s .  Firs t ,  values of 4 for points immediately 
below the f ree  surface a r e  obtained by Equation 12. Next, mesh  points 
a r e  revised column wise f rom left to right, with the direction in each 
column pointing upward, according to  Equations 3 through 9. 
procedures a r e  ca r r i ed  out alternatively until the maximum e r r o r  in 4 
at any point drops below the previously se t  tolerance,  usually around l o e 6 ,  
and i teration is terminated. Velocities a r e  then computed, thus com- 
pleting the calculations of the flow field on a new time level. 
increment is again selected, and the whole procedure is repeated. 
The value of the increment depends 
Then 
These two 
A new t ime 
14 
CONCLUSIONS 
Several  tes t  runs have been made to  determine the validity of 
Young's formula for predicting the value of the optimum relaxation factor.  
Assuming a flat surface with constant 4 a c r o s s  the width of the tank, the 
flow field is laid with 44 X 38 mesh points. The best  value for  o was found 
to be around 1.94 which agrees  well with the value of 1. 92 given by Young's 
formula.  With the value of 1 .94 f o r  w ,  it took approximately 200 i terations 
to converge to  
ations for  subsequent t ime levels with the same t ime increment.  
in for the first t ime increment and only a few i t e r -  
Running t ime on the UNIVAC 1107 computer was l e s s  than 0. 5 sec  
per  sweep, which is tolerable since on the average it does not require  
many sweeps for convergence. 
As was shown previously, an explicit finite-difference scheme 
was employed to solve the differential equations for the f r ee  surface.  
This could cause instability i f  time increments  were not properly chosen. 
With the presence of the mixed derivative and nonlinear t e r m s ,  there  is 
no known theoretical  method to obtain the stability c r i t e r i a  for such a 
difference scheme. 
numerical  experiments.  Since central  t ime differences were  used, i t  is 
conceivable that the difference scheme may always be unstable no mat te r  
how sma l l  the t ime increment is, a s  i n  the study performed by 
Richardson8. 
o r  DuFort and Frankel  difference scheme would have to be used. It i s  
unknown, however, what effects the mixed derivative and nonlinear t e r m s  
will  have on these difference schemes. 
The stability c r i t e r i a  a r e  to be determined only by 
If this should happen, either a forward t ime difference 
I 
15  
i 
i 
REFERENCES 
1. Forsythe,  G. E. and W. R. Wasow, Finite-Difference Methods for 
P a r t i a l  Differential Equations, John Wiley & Sons, Inc. , January  
1964 
2. Ingram. E. H . ,  "Equations Governing the Two-Dimensional Dynamic 
Behavior of a Liquid Surface in a Reduced Gravitational Field", Brown 
Engineering Company, Inc. , Technical Note R-165, October 1965 
3. Geiger,  F. W . ,  "Hydrostatics of a Fluid Between Pa ra l l e l  P la tes  at 
Low Bond Numbers", Brown Engineering Company, Inc. , Technical 
Note R-159, October 1965 
4. Beckenbach, E. F., Modern Mathematics for the Engineer,  University 
of California Engineering Extension Ser ies ,  McGraw-Hill Book 
Company, Inc. ,  1961 
5. Crandall ,  S. H . ,  Engineering Analysis,  McGraw-Hill Book Company, 
Inc . ,  1956 
6. Scat ter lee ,  H. M. and W. C. Reynolds, "The Dynamics of the F r e e  
Liquid Surface in Cylindrical  Containers under Strong Capi l lary and 
Weak Gravity Conditions", Stanford University, Technical Report  
LG-2, May 1964 
7. Harlow, F. H. ,  E.  L. Young and J. E. Welch, "Stability of Dif- 
f e rence  Equations, Selected Topics", Los Alamos Scientific 
Laboratory,  Report  No. LAMS-2452, July 28,  1960 
8. Richardson, L. F . ,  "The Approximate Arithmetical  Solution by 
Finite Differences of Physical Problems Involving Differential 
Equations, with a n  Application t o  the S t r e s ses  in a Masonry Dam", 
Philos. Trans .  Roy. SOC. ,  London, Ser .  A, Vol. 210, pp. 307-357, 
1910 
16 
I O R I G I N A T I N G  A C T I V I T Y  (Corporate author) 
Res ea r  ch Labor ator ie s 
Huntsville, Alabama 
Brown Engineering Company, Inc. 
2 0  R E P O R T  S E C U R I T Y  C L A I S I F I C A T I O N  
Unc las s i f  ie  d 
2 b  GROUP 
N/A 
Young, Car l  T. K. 
i REPORT DATE 7.. T O T A L  N O .  O F  P A G E S  
October 1965 22 
7 b .  NO.  O F  REPS 
8 
I.. C O N T R A C T  OR GRANT N O .  
NAS8-20073. . 
b. P R O j t C f  SO. 
d.  
9. ORIOINATOR'S REPORT NUMDLR(S) 
T N  R-164 I 
None 
c .  N/A Pb. OTHER R PORT NOfS) (Anyofhernumbare ha1 may ba aiaI&od I him reporJ 
11. SUPPLEMENTARY NOTES 12. SPONSORING MILITARY ACTIVITY 1 P & V E  Laboratory 
None I Marshall  Space Flight Center 
3. ABSTRACT 
This repor t  descr ibes  a numerical  procedure for 
solving a two-dimensional, unsteady flow problem. The 
fluid is in a tank, has a f r e e  surface, , a  periodic source 
(to be terminated) on the side of the tank, and is subject 
to a near  zero  gravitational field. The flow is assumed 
to  be incompressible and irrotational so that the prob- 
lem reduces to  a boundary value problem governed by 
the Laplace equation. 
14. KEY WORD8 
Numerical analysis 
Sloshing 
Hydrostatics 
Fluid mechanics 
F r e e  surface 
17 


Analysis of two-dimensional, unsteady flow in a propellant tank under low gravity by finite difference methods

Abstract

Similar works

Full text

Available Versions

NASA Technical Reports Server