NOV. iiz
ARCHi Ef
SYMPOSIUM ON
"HYDRODYNAMICS OF SHIP AND OFFSHORE PROPULSION SYSTEMS"
H0VIK OUTSIDE OSLO, MARCH 20. - 25., 1977
"NUMERICAL METHODS APPLIED TO FULLY CAVITATING
FLOWS, WITH EMPHASIS ON THE FINITE ELEMENT METHOD"
By
R.L. Street and P.Y. Ko Department of Civil Engineering
Stanford University Stanford, California
SPONSOR: DET NORSKE VER ITAS
Ref.: PAPER 15/8 - SESSION 3
Lab. y. Scheepsbouwkund
Technische Ho9eschool
NUMERICAL METHODS APPLIED TO FULLY CAVITATINC FLOWS, WITH EMPHASIS ON THE FINITE ELEMENT METHOD1
R. L. Street and P. Y. Ko
Department of Civil Eninering
Stanford University Stanford, California 94305 U. S. A. A BSTRACT I Bibliotheek vn AfdIing Sclìee.s DOCUMENiÀTE
I: V3&- Qjt h
DATUMiA number of finite-difference and finite-element methods have been applied to free-surface flows. Two- and three-dimensional flow solutions by finite-difference methods are described, as are two- and
three-dimensional finite-element solutions. However, no fully non-linear solution for three-dimensional fully cavltating flows exists. A finite-element method for handling such flows is described. It is an iterative scheme for finding the f lowf leid velocity potential and the location of the free surface. Examples are given showing (a) quantitative agreement between the finite-element method and the finite-difference technique
of Brennen for axisyrninetric cavity flows and (b) simple three-dimensional flows past disks in a square water tunnel.
1,
INTRODUCTIONNo one has yet solved a fully non-linear three-dimensional fully-or super-cavitating flow. However, a number of authors have solved
fully non-linear two-dimensonai axisymmetric and near axisymmetric flows by numerical methods. This paper focuses on techniques that can be applied tò these problems and on the progress being made by the authors using a
finite element method to solve fully three-dimensional flows. We discuss
'This research was carried out under the Naval Sea Systems Command General Hydromechanics Research Program Subproject SR 023 01 01, administered by the David W. Taylor Naval Ship Research and Development Center, Contract
-2-Bibliotheek van de
Afdeling Stheepso_ en Schpvarunde
Tech -r-c!e
HOC
De!
DOCUMENTTIE
T-first the relevant literature and then the f ormulat on
particular problem, including first results and a comparison with another numerical solution of an axisymmetric cavity flow.
Our goal has been to develop a computer program for the simulation of three-dimensional fully cavitating flows about hydrofoils beginning with a problem formulation for two-dimensional flow past plates which we used with a finite difference method [1]. We have attacked the flow past disks in a water tunnel. The following review of recent work allows us
to discuss differences and similarities between various methods and to put our present studies in perspective.
2.
REVI W OF TECHNIQUESThere are several major categories of numerical techniques applied to free surface flows. There are finite-difference and finite-element methods. Amongst the finite-difference methods we also find formulations in which
the problem is solved in a stream function and velocity potential space or in which the formulation is based on the distrIbution of singularities on the free surface of the flow (normally called a boundary integral technique). The application of all these techniques to both two- and three-dimensional free surface flow problems is described below.
2.1. Two-Dimensional Methods
Our work (lJ in numerical simulation began with a study involving two-dimensional cavitating flows. We used a fInite-difference technIque and solved a problem of flow past a disk in a water tunnel by employing a Riabouchinsky model and imaging to achieve an exact formulation. Irregular finite-difference stars were employed along the curved free
grid in the neighborhood of the separation point on the flat plate.
The solution of the finite-difference equations for the velocity components was achieved by successive over-relaxation. Because the velocity components were employed as the unknowns, the solution required simultaneous
satis-faction of the Laplace equation for each velocity component and of the non-linear free surface boundary conditions. The solution converged and gave reasonable results. However, the cost for the two-dimensional case was equal to the cost for a solution of similar accuracy for a three-dimensional disk-In-water-tunnel flow by the finite-element method.
Brennen [2] combined a finite-difference technique with a mapping to the stream function--velocity potential space to solve the axisynunetric flow past a circular disk in a circular water tunnel using a Riabouchinsky model. A difficulty with the extension of this technique to more general bodies is that the solution is inverse because the stream
function--velocity
potential space solution must be mapped back to physical space to obtain the final geometry. In Section 3 below we will make a comparison of our finite-element solution with Brennen's results.White and Kline [3] developed a general method for the solution of turbulent separated axisymmetric flows. In their general development they employed a new boundary integral technique for the potential f lowfield outside the boundary layers developing in their diffuser flows. Unlike grid techniques, the boundary integral method requires computation of unknowns only over the flow boundary. The solution technique obtains the potential flow from the numerical solution of Green's third identity and iteratively solves for the correct location of an initially unknown free surface. In principal this technique can be extended to general three-dimensional flows. Interestingly, a very similar technIque was developed
ri
in 1953 by Armstrong and Dunham [4]. They applied a vortex sheet singularity to free streamline flows, thereby generating an integral equation whose
solution formed part of an iterative procedure to locate the initially unknown free surface. They solved several axisyinmetric cavity flows in an
infinite fluid field and obtained reasonable agreement with experiment. The application of the finite-element method to free streamline flows was first made by Chan, et al. [5,6]. They showed that the finite-element method could be used to solve for the velocity potential in two-dimensional
and axisyrnmetric ideal fluid flows involving a free surface. The technique involved is assumption of the location of the free surface, solution of the resulting well posed problem, and relocation of the free surface according to the free streamline boundary condition. Sarpkaya and Hiriart [7] applied a similar technique to solve the free streamline problem of flow in
curved jet deflectors.
2.2 Three-Dimensional Methods
Jeppson [8] developed an inverse formulation for
three-dimensional flows in which he defined a velocity potential and two additIonal functions which defined the flow paths. By changing the conventional roles played by the variables of the problem, he converted a free surface with an unknown position in physical space into a piane of known position in the
inverse stream function--velocity potential space. The technique is very similar to that of Brennen [2] sited above and suffers the saine disad-vantage, namely that the shape of curved solid bodIes cannot be prescribed
in advance in physical space.
Larock and Taylor [9] applied finite-element techniques to solve the jetfiow from a circular pipe and orifice under the influence of gravity.
4-The resulting flow is slightly three-dimensional, namely, the jet remains essentially circular but droops under the action of gravity creating a three-dimensional flowfield. They employed the velocity potential as the dependent variable and allowed it to vary quadratically within each
isoparametric hexahedron element. An iterative approach similar to that used in other free streamline problems mentioned above was used to
determine the free surface location whose precise location is initially unknown. The formulation of their problem is relatively straightforward. and they are able to prescribe the f lowrate in the pipe and then proceed with their iteration to find the location of the stream surface. We will
see below that the cavity flow problem is not as easily specified, leading to one of the major complications and differences between jet and cavity flows.
In sunnnary, the techniques of Mogel and Street [1], Brennen [2J, Sarpkaya and Hiriart [7], and Larock and Taylor [9] are of direct interest to us. Mogel and Street's work is of Interest because it is part of our
own beginnings and illustrates the problems associated with free boundarIes in finite-difference schemes. Brennen found good agreement between his solution and experiments. Hence, we duplicated his solution with our three-dimensIonal method as a check. Sarpkaya and Hirlart and Larock and Taylor illustrate the significant difference between jet and cavity flows as well as providing illustration of the application of the finite-element method to other problems. Our three-dimensional computer program bears a strong resemblance to that of Larock and Taylor because we both started with the same three-dimensional groundwater flow program developed at the Civil Engineering Department, University of Wales, Swansea, Wales.
3. THE PRESENT APPLICATION -- FULLY cAVITATING FLOWS
To demonstrate the use of a finite-element technique we wish to
formulate and solve the problem of a fully cavitating flow past a disk in
a water tunnel. For our analysis the flow la assumed to be incompressible, steady, and irrotational. Hence, it is governed by a velocity potential
(x,y,z) which satisfies the Laplace Equation = O within the fluid and generally
y = -vq (1)
so, in the usual notation and Cartesian coordInate system,
u
-/3x
V =
-3/y
(2)w = -a4/z
We have a free surface present so the streamline conditions
dx
u'
dx u(3)
are useful in describing the free surface. We neglect gravity influence so on the free surface
2 2 2 2
(4)
u +v +w
q constantAs a consequence, in terms of a coordinate system oriented along the streaxalire, the change of potential on the free surface is given by
- q (5)
3.1 The Finite Element
Method
We can now proceed with the formulation of the problem, but first should answer why are we using the finite element method. There are at least three reasons. First, f mIte elements routinely handle irregular
-6-(curved) boundaries. Second, finite elements routinely handle the
variation in mesh size so necessary to the solution of free surface flow problems. Third, finite elements produce a computationally efficient computer code.
In finite differences we represent a problem solution as a discrete set of values defined at a set of node points; partial differential equations are replaced with difference equations. In the finite-element method we represent a problem solution as a continuous functIon in terms of values at the edges of finite volumes in the problem domain and with a specified variation locally within such volumes as functions of the edge value. Partial differential equations are replaced by integrals which are subsequently replaced by a set of algebraic equations which can be solved directly and exactly.
There are three parts to a finite-element method formulation: (a)
the varIational principal.
Functional: x
f
1fr2
()2
(\2l
V
2
L\0x1
\y
z) j
in which the functional x is min:Lmized with respect to the velocity potential values specified at finite element nodes 4. . (b) The
functIonal representation within an element.
Within element: 4 = [N1, N2, --- N20] 2 = [N) {4}e (7)
20
We employ a quadratIc variation within the hexahedral 20-node element. In Equation 7 the are the values of the velocity potential at t1e
corner and mid-side nodes of the element, while the N1 are the well-known shape functions. (c) The isoparametric element. We have chosen,as noted
above, a 20-node hexahedron element. Figure 1 shows a schematic of our isoparametric element in the physical Cartesian space and in its trans-f orined space.
Now the functional in Equatton 6 is minimized with respect to the nodal values within each element (10]. Thus,
e
3N N 3N N. N N. e=
f
__t
+{} av -
0 (8)or for each
element
e
____ = [h]f}
o (9)The global sum over all elements produces a linear algebraic equation set
(H][fl
0 (10)A unique solution is guaranteed by specification of the or derivatives along the flow boundaries.
For the isoparametric element of Figure 1 we write a transform from the physical space to a set of natural local coordinates (,n,ç) such that each element occupies the cube bounded by i . Generally
then [10,11] =
Ni(,n,)x
1=1 = i= 1-8-z(,n,)
i= i
where
X1, y, z, and
are the nodal values In each element. The key idea is that the mathematical formulation is easily made on a set of cubes while physically we can handle complex curved boundary surfaces.The essence of the transformation (Equations 11) is in the Jacobian
ax/ac ay/a az/ag
[J] = x/ari Z3y/a az/an
3x/3 ay/3 az/ac
which
relates, for example,dxdydz = IJIddnd
SV
3.2 The Problem Formulation
Given the basic mathematics of the finit element method, we can turn to the physical problem. Figure 2 is a schematic of our initial
three-dimensional problem for flow past a disk in a square water tunnel. [The axisytumetric case is obtained by using a circular water tunnel
boundary.] Typical boundary conditions and geometry are also shown. and aN1 = [J] (14) aN1 3x an 3N1 ay aN1
Figure 3 shows a cross section along the flow to illustrate the imagIng used In the Riabouchinsky model flow. Note from the symmetry that we are required to solve only one-eighth of the total flowfjeld for this
simple case.
The boundary conditions for a fully cavitating flow in a water tunnel are both straight forward and restrictive. We see in Figure 2 that with a Riabouchinsky model, the upstream and downstream boundaries are equipotential lines. If the potential downstream is taken to be zero, then when a free surface shape is known (assumed or calculated) the potential can be computed on the free surface by Equation 5. The remainder of the flow boundaries are solid or no flow boundaries. The key question then becomes, "How do we select
F which is the upstream potential value?" It cannot be selected arbitrarily because in our confined flow there is only one physically-allowable relationshIp
between the upstream f lowrate or velocity (U) and the free surface
speed (q) when the flow geometry is specified (i.e., for given P, W, F and L in Figure 1). In other words, if P, W, F and L are given, there Is a unique value of the cavitation number
/ p
-p
-1=
UstU 1 1 2
\ USI Pl.Tus
in which p is the upstream pressure,
PC is.the pressure in the cavity and p is the fluid density.
In our case we have elected to select through an iteration procedure in which is adjusted until the velocity at separation from the disk is equal to q . [We actually set
q E i
without loss ofgenerality.]
(15)
'-lo-3.3 The Problem Solution Scheme
With the given boundary conditions (cf., Figure 2), an assumed free surface location, and specified geometry of disk and tunnel, the finite element method can be used to find the nodal values
of the velocity potential and to determine the fluid velocities u, y and w . We next move the free stream surface to a better location and adjust as needed.
The total process can be summarized as follows:
Establish the geometry. Input all data including locations of all corner nodes, assumed form of free surface, assumed etc. The computer program then locates all the mid-nodes.
Generate "stIffness" matrix H (Equation 10). Account for various boundary conditions where or c derivatives are specified
and for isoparametrlc mappings.
e. Solve for by Gaussian elimination and back substitution.
Calculate fluid velocities.
Move free surface if the separation velocity is close enough
to unity. Movement is accomplished by integration along free streamlines beginning at element corner nodes on the edge of the disk. This estab-lishes a set of lines which are tangent to the velocity in accordance with the streamline equations
dx u
A cubic Lagrangian interpolation is employed in planes x = constant to locate the element transverse mid-nodes which are missed by the integration process.
dz
w
and
-1.2-Establish new values of and of on the free surface. If surface was not moved in Step e, only the new
us is calculated; is decreased if the separation velocity is greater than unity.
Modify the "stiffness" matrix H if the free surface was moved in Step e. Only elements near the free surface are changed by
the move so only a small part of H is changed.
Return to Step c unless free surface movement (if Step e is executed) is less than an arbitrary, preset, small amount. If this happens the solution is complete and the final results are printed. Intermediate printout is made for diagnostic purposes as the solution proceeds.
The solution program has been implemented on the IBM 360/91 and 370/168 Triplex System at the Stanford Linear Accelerator Center. With 1325 nodes and 224 elements in a typical case, one iteration (Steps a through h) takes about 45 seconds and solution is achieved in about 15 iterations.
3.4 Results
We show results for two geometrics. The first is the case shown in Figure 2 for a square water
tunnel.
The second is the axi-synmietric flow past a disk in a circular water tunnel.Figure 2 shows the layout of our first example. Let the half-length L 5 , the half-width W = 10 , F 10 and the disk radius P = 2
Figure 4 gives cross-sectional plots for the cavity shape and the velocity field. FIgure 5 shows the effect on the cavity cross-section shapes as the tunnel half-width is reduced. The points for each case are the cavity shape at various distances downstream from the dIsk. Interestingly, the cavity cross section remains almost circular, except
for W/P = 2.5 , in spite of the square tunnel cross section. However,
the wail effect is substantial as illustrated in Figure 6. There we use the centerline cavity width B (at y O). It is hardly affected for WfP > 2.5
Finally, we examined a purely axisymmetr5.c case, whose solution is well known [2], in order to validate our method. Figure 7-shows
reproductions of two plots of cavity-tunnel characteristics as functions of a as given by Brennen [2]. We ran two test cases, viz., WfP = 7.5,
F = 20, P = 2 and L/P = 5 and 3. Our results which agree quantitatively with Brennen's are shown by the cross and circled dot respectively. The variables (calculated, not prescribed values) ¿n the plots are B and a for which we agree with Brennen's results to within 2 percent.
4.
PROGNOSISOur immediate goal is to simulate the cavity flow past an elliptic flat plate hydrofoil used in water tunnel experiments by Leehey at the Massachusetts Institute of Technology. However, not all our tests with fully three-dimensional flows have been successful.
The key problem seems to lie in the difference between jet-like flows, such as those solved in Ref s. [6,7, and 9], and bounded-interior
(cavity-like) flows, such as solved by White and Kline [3] and the present authors. In a jet flow one specifies theflowrate and upstream velocity; the only iteration then is to move the free surface until convergence is achieved. In bounded-interior (or cavity-like) flows, two parameters must be iterated and a feedback ioop is involved, e.g., where moving the free surface alters the relationship between
Us and the disk separation velocity in our cases. As a consequence of
this feedback, 'White and Kline [3] discovered that they had to grossly undershift their free streamline to prevent instability, i.e., they moved it only 10 to 20 percent of the distance indicated by the shift integration (See Equations 16) at each iteration.
We are just beginning to explore these problems in our cases, but we believe that the prognosis is good for obtaining a method applicable
to general three-dimensional fully cavitating flows.
'-14-REFERENCES
MOGEL, T.R, and STREET, R.L.. A Numerical Method for Steady-State
Cavity Flows.
J. Ship Res.
l8(l974):l,pp 22-31.BRENNEN, C.
A numerical solution of axisyetric cavity flows.
J. Fluid Mech.
37(l969):4,pp 671-688.WHITE, J.W., and KLINE, S.J.
A Calculation Method for Incompressible
Axisyminetric Flows, Including Unseparated, Fully Separated, and
Free Surface Flows.
Stanford, Dept. Mach. Eng., Thermosci. Div.,
Rept.
MD-35(1975).ARMSTRONG, A.H., and DUNhAM, J.H.
Axisynnnetric Cavity Flow.
Kent,
England, Arm. Res. Estab., Rep.
12/53(1953).CHAN, S.T.K., and LAROCK, B.E.
Fluid Flows from Axisymnietric
Orifices and Valves.
J. Hydr. Div., Proc. ASCE
99(1973):HY1, pp 81-97.CRAN, S.T.K., LAROCK, B.E., and HERRMANN, L.R.
Free-Surface Ideal
Fluid Flows by Finite Elements.
J. Hydr. Div., Proc. ASCE
99 (1973) :HY6,pp 959-974.
SARPKAYA, T., and
HIRIART, G.Finite Element Analysis of Jet
Impingement on Axisymxnetric Curved Deflectors.
Finite Elements in
Fluids, V. 1, Wiley, 1975.
JEPPSON, R.W.
Inverse Solution to Three-Dimensional Potential
Flows.
J. Engin. Mech. Dlv., Proc. ASCE
98(l972):EM4,pp 789-812.LAROCK, B.E., and TAYLOR, C.
Computing Three-Dimensional Free
Surface Flows.
mt. J. Num. Meth. Engin.
lO(1976),pp 1143-1152.HIJEBNER,
ICH.
The Finite Element Method for Engineers.
New York,
ERGATOUDIS, I., IRONS, B.M., and ZIENKIEWICZ, O.C. Curved, Isoparametric, "Quadrilateral" Elements for Finite Element Analysis. mt. J. Solids Structures 4(l968),pp 31-42. ZIENKIEWICZ, O.C. The Finite Element Method in Engineering
Science. London, McGraw-Hill, 1971.
'-16--Figure Captions
FIGURE 1. The Isoparametric Element.
FIGURE 2. Schematic of Initial Cavity Flow Problem. FIGURE 3. Imaging for Riabouchinsky Model.
FIGURE 4. Longitudinal Cross-Sections for Flow in a Square Water
Tunnel. L=5,F=lO,W=L0,andP=2.
FIGURE 5. Variation of Cavity Shape with Tunnel Width and Plate Diameter Ratio. P = 2; cross-sectIons just off disk and at cavity mid-section.
FIGURE 6. Variation of and Cavity Width with Relative Tunnel Width--the Wall Effect.
FIGURE 7. Comparison of Simulation of Axisynmietric Cavity Flow by Three-Dimensional Model and Brennen's Results [2]. For Three-DimensIonal Model: c 0.33: W = 15, L = 10,
F 15, P = 2 and o 0.43: W = 15, L 6, F = 15, P = 2.
FIGURE 1.
'A
ò,/òz:O
(nfl
fIti Aò4/äz=O
(8
<
Yx$NH0
coii
rn
k n w'\ \\\\
ò4/òz =0
(no flux)
Disk
FIGURE 2.Schematic of Initial Cavity Flow Problem.
Typical
finite element
(this is a quadratic element:
¡t has 20 nodes and general
parallelepiped shape)
F= IO
Forward
IDisk
(Image
Cavity
Region i//it
/
Solved in 3-D
/ FEM Method
RIABOUCHINSKV MODELFIGURE 3.
r
7 4
X
Longitudinal Cross-Sections for Flow in a Square Water Tunnel.
L 5, F = 10, w 10, and P Io II 12 13 14 5
o
2
4
3 I I t FIGURE 5.Variation of Cavity Shape with Tunnel Width and Plate Diameter Ratio.
P = 2; cross-sections
ju3t
off disk
and at cavity mid-section.
o
2
34
5
¿
W=5
uW: 7.5
O WIO
1.8 1.6 1.4 1.2
i0
0.8
0.6
0.4
0.2
o
0
2
4
6
8
Io
12w/P
FIGURE 7.
o
0.8
0.7
0.6
0.5
P/L
0.4
0.3
0.2
0.I
a t IQ c0.43
0.33
Curves are from Brennen [2].
0.6
0.5
0.4
I I I I I I E5
IO 1520
wfp
0.10.2
0.3
0.4
o.Comparison of Simulation of Axisymmetric Cavity Flow by Three-Dimensional Model and Brennen's Results [2]. For Three-Dimensional Model: o 0.33: W = 15, L = 10,
F 15, P ' 2 :niI 0.43: W IS. L= 6, F =
