- Email: [email protected]
The numerical solution of planar and axisymmetric cavitational flow problems
Computers & Fluids Vol. 12, No. 1, pp. 55-65, 1984 Printed in Great Britain.
0045-7930/84 $3.00+ .00 Pergamon Press Ltd.
THE NUMERICAL SOLUTION OF PLANAR AND AXISYMMETRIC CAVITATIONAL FLOW PROBLEMS JOYCE M. AITCHISON Oxford University Computing Laboratory, 19 Parks Road, Oxford OXl 3PL, England (Received 30 October 1981; in revisedform 6 September 1982) Abstract--This paper describes a method for the numerical solution of planar and axisymmetric Riabouchinsky cavity flows. The problem is formulated as the minimization of a functional over a variable region, which is solved numerically by the method of variable finite elements. The convergence difficulties experienced by some other authors do not appear here, and the method is found to be successful for both planar and axisymmetric problems.
1. INTRODUCTION 1.1 The model In this paper we consider the cavitational flow of an incompressible inviscid fluid past a plate or disk in a channel of finite width and infinite length. The fluid flows from left to right (Fig. la) and has a uniform velocity q0 far upstream. It meets a plate CK placed symmetrically in the channel at right angles to the flow. The basic problem is to find the flow pattern behind the plate. Immediately behind the plate a cavity is formed containing air at a constant pressure: our primary objective is to find the shape and size of this cavity. It is reasonably clear from experiment what sort of shape the cavity should have at its left hand end. However, it is not at all clear (from experiment or theory) what form the cavity should take at the right hand end. This phenomenum, known as the closure condition, is still much discussed. Various models have been put forward (see, for example, Birkoff and Zarantello [1], Gilbarg[2] or Wu [3]). They include an infinite cavity length and a re-entrant jet. Here we consider the model due to Riabouchinsky[4], in which we suppose that an image plate C'K" can be placed in the flow at some point downstream from the original plate and that the flow geometry will be symmetric about the line midway between CK and C'K" (Fig. lb). Since the flow is also symmetric about the axis of the channel, we need only consider the quarter of the flow region shown in Fig. 2. The lines AG and EG are axes of symmetry. Previous numerical solutions using this model include those by Mogel and Street [5] for planar flow and by Young et al. [6], Arms and Gates [7], Brennen [8] and Fox and Sankar [9] for the axisymmetric problem. We will first describe the case of planar flow and then give the corresponding equations for axisymmetric flow in Section 1.4. 1.2 The differential equations for planar flow In forming the governing differential equations of the flow we will make use of Bernoulli's equation
p + ½q2= constant where p is the pressure and q the velocity of the fluid. Along the free surface CD, p = 0 and so ½q2 = constant = ½qc2 (say). Thus p + ½q2= lqc2 everywhere. Inside the fluid p i> 0 and so q2 ~
56
J.M.
AITCHISON
(a)
I
I I
(bJ
Fig. 1. (a) A general cavity flow. (b) The Riabouchinsky model.
E
I
C
I I i
I
B.
L
,G
Fig. 2. The region of solution.
number, tr, which is defined by qc 2 - -
O-~.
qo2
q0 2
The stream function ~k(x, y ) satisfies the following differential equation, where I) is the flow region A B C D E F (Fig. 2)
~2~ a2~ dx---~+tgy--~=0 in~
The numerical solution of planar and axisymmetriccavitational flow problems
57
subject to: qoY
d~ =
on A F ( x - - ~ - o o )
= 0
on AB, BC onFE
~b = qoH
t~b ?x
0
on ED
and q;=O on CD, - -
-
-
qc
~n where n is the direction of the outward normal from Q. Note the double boundary condition on the free surface CD. The lengths H and d will be specified together with the upstream velocity q0, but there is some freedom in the specification of another parameter to fully define the problem. The natural choice from a physical point of view is to fix qc and let the half-length of the cavity, L, be unknown. This is equivalent to specifying the pressure far upstream. However, small changes in the pressure lead to large changes in the cavity length and so this formulation is numerically unstable. Hence we prefer to fix the length L and let q~ be an unknown constant to be determined as part of the solution. The problem can be made non-dimensional by putting x*=x/H,
y*=y/H,
~k*=~b/qoH
d* = d/H,
L* = L / H ,
qc* = qc/qo.
The equations become (dropping the *'s) in f~
3x z +~-~-yZ= 0
(1)
subject to =y
on A F
= 0
on AB, BC
= 1
on FE
t~-x= 0
on ED
ag,
(2)
and
O~_ On
on CD qc
where the region t2 is still that shown in Fig. 2, but with A F The cavitation number is now given by a = qJ - 1.
= 1.
(3)
The numerical problem is now to solve eqns (1) and (2) for ~b(x, y), qc and the free surface C A F Vol. 12, No. I - - E
58
J . M . AITCHISON
CD. In fact eqns (1) and (2) can be solved numerically for a range of values of qc. Only one of these corresponds to the true flow pattern. The technique for finding this particular value of q,. will be described in Section 4, for the present we will assume that q, is fixed. 1.3 The associated variational principle Let
+ \ O y / i d + - ~ q c 2 dxdy.
E[~,, ~] =
(4)
The function ~b(x,y) and the region f~ which satisfy eqns (1) and (2) can be found by minimizing the functional E with respect to both ~b and f~. f~ is allowed to vary by moving the line CD, where the point C is fixed and the point D moves along the line x = L. is subject to the constraints ~k = 1 on FE
(5)
on AF
=y ~b=O
o n A B , BC, CD.
Under these conditions it is easily shown (Courant and Hilbert[10], pp. 260-262) that E has a stationary value when • ~
t3x
= 0
in f~
- 0
on ED
and (V~b)2 = q,,2
on CD.
Thus eqns (1) and (2) are' fully satisfied. It can be further shown that this stationary value is in fact a minimum. This minimization problem will be solved numerically by the method of variable finite elements. 1.4 Axisymmetric flow For the case of axisymmetric flow we take coordinates (z, r) replacing (x, y). The non-dimensional equations corresponding to (1) and (2) are ~2~
1 O~O 02~/ r Or + ~ = 0
Or2
in f~
(6)
subject to = r 2 on AF ~b=O
a~O Oz
o n A B , BC
= 1
on FE
0
on ED
and #/=0 ~ r On 1
on CD. qc
(7)
The numerical solution of planar and axisymmetrie eavitational flow problems
59
In this case qc 2
a=~
1.
(8)
The solution to eqns (6) and (7) corresponds to a stationary value of E[~h,D]=
f ~ f {12 r LFfto~h~ \ t o r ] 2 + (to@)2] ~z + ~,q c 2r } drdz.
is allowed to vary by moving
(9)
CD as before, and ~h is constrained by
FE ~ = r 2 on AF ~O= 0 on AB, BC, CD. = 1
on
(1o)
Again it is easily shown that E has a stationary value when tO (1 to~b'~
to (1tomb'] 1 tomb
r toz
in f~
- 0
on
EL)
on
CD;
and ~ (\~r j
\ t o z ] j = qc2
and so eqns (6) and (7) are fully satisfied. 2. T H E M E T H O D O F V A R I A B L E F I N I T E E L E M E N T S
The finite element method for the numerical solution of variational principles over fixed domains is well-known (see, for example, Strang and Fix[ll] or Zienkiewicz[12]). However, many problems in fluid dynamics can be expressed as a variational principle over a variable domain, typified by the example given here. The numerical solution of such problems can be obtained by the method of variable finite elements. The region ~ is covered by a net of elements: we use triangles here, although the arguments are equally applicable to quadrilaterals (but not rectangles). The function ~, is approximated, as usual, by piecewise polynomials, the coefficients of which are chosen to minimize the functional E[~, ~)]. The region 1) is varied by distorting the geometry of the elements. This is achieved by allowing the coordinates of the vertices of the triangles to vary. However, care must be taken that we do not allow too many degrees of freedom in this geometric variation. If we allow the (x, y) coordinates of all the vertices to be free then we solve not only for the position of CD but also for the optimal position of nodes within ft. Under certain conditions this might be worthwhile, but here the problem would become too large. The approach used here is that all x coordinates are kept constant and only a limited number of y coordinates are allowed to vary. The linkage structure between elements remains fixed.
2.1 Specification of the geometry The initial triangulation is performed on a region where the line CD has an arbitrary, simple form, for example the line y = d (Fig. 3). For the present work this region was triangulated using an interactive graphics program on the PRIME 400 computer at the Rutherford Laboratory. This program is described in detail in Aitchison[13]. The triangulation can be described by the following data:
60
J. M. AITCHISON F
E
A
B
Fig. 3. A basic grid.
NP:. number of vertices. NE: number of elements. X, Y[I:NP]: arrays giving the (x, y) coordinates of each vertex. NOP[I:NE, 1:3]: integer array giving the three vertices of each element. Each subsequent position of the free boundary is specified by a set of ordinates {h~}7=~. Figure 4 shows the example n = 3. The triangulation for this region is generated from the basic triangulation as follows. (i) The array NOP giving the linkage structure remains the same. (ii) The array X giving the x coordinates of the vertices remains the same. (iii) The entries in the array Y which refer to points to the left of the plate remain the same. (iv) The other entries in the array Y are scaled linearly as follows. Consider the vertex (x0, Y0) in the basic grid. Let h0 be the height of CD at x = x0 in the basic grid. Let h~ be the corresponding height in the current grid. Then the vertex (x0, Y0) becomes (x0, Y0 where Yl=l
( 1 - - hi) --(l--y0). (1 -- h0)
The most detailed grid used in this investigation is shown in Fig. 5.
.B
Fig. 4. A distorted grid.
LC"~
Fig. 5. A typical grid.
The numericalsolution of planar and axisymmetriecavitationalflowproblems
61
2.2 The discrete minimization problem The function ff is approximated by piecewise linear functions for the case of planar flow, nodes being taken at the vertices of the triangles. However, for the axisymmetric problem quadratic functions are more appropriate since ~b behaves like r 2 away from the plate. Nodes are taken at the vertices and at the midpoints of the sides of the triangles. An expression for the discrete approximation of the functional E[ff, f~] can now be written in terms of the geometric v a r i a b l e s {hi}~=1 and the nodal values {ui}~=l which approximate to ~ at those points where it is not directly prescribed by boundary conditions. The discrete functional Ed is of the form /
Ea(u, y) = ~ u A (y)u + b'(y) • u + c(y)
where the matrix A is symmetric and positive definite. For stationary values we require
~E~ Ou~
=0
i=l,...,s
(11)
i = 1. . . . . n.
(12)
and
~E~
--=0 Oh~ Equations (11) give
d(h)u= -b(h)
(13)
which is linear in u, but eqns (12) are difficult to derive analytically since they involve the differentiation of a complicated function of h. An alternative method is as follows: for given h, eliminate u from Ed by solving eqn (13). Then Ed(U, h) -- E*(h) = lU'A (h)u + b'(h) • u + c(h)
= ½b'(h) • u + c(h)
where u = - A (h)- ~b(h). We now find the value of h which minimizes E*(h) by a quasi-Newton method due to Fletcher[14]. This optimization problem is much smaller than the original one since the majority of the unknowns have now been eliminated. 2.3 Efficient solution of the linear equations The set of linear eqns (13) has to be solved for many different values ofh. Various steps are taken to ensure that this is done as efficiently as possible. The nodes on the finite element grid are numbered so that ni < nj only if xi <<,xj. Thus the equations relating to those nodes to the left of BC appear in the top part of the matrix A. Let these be the first n~ rows. The elements of A in these rows do not depend on h and so do not have to be recalculated at each stage. The equations (13) are solved by Cholesky decomposition, A = L L r. Since A is symmetric, only its lower half is stored and the matrix L overwrites this. On the first iteration of the minimization process, the whole of A is calculated and decomposed. Subsequently only those rows beyond row n 1 are calculated and decomposed. Thus the first n~ rows remain in decomposed form. The matrix A is banded, but has variable bandwidth due to the irregularity of the nodes. It is stored in a one-dimensional array in variable bandwidth form.
62
J . M . AITCItISON 3. T H E S I N G U L A R I T Y
AT THE SEPARATION
POINT
In both the planar and axisymmetric problems there is a mild singularity in the function at the separation point C. In the following formulae q -~ q,.
3.1 Planar flow Take a local coordinate system (4, r/) centered at C. Thus =x
and
r/=y-d.
Define local polar coordinates (p, 0) so that =psin0
and
q=pcos0.
Note that this is not the usual definition of 0. A local series solution for ff can be found by standard techniques. A first approximation to the local solution for ~ is: ~ =q
- p sinO
+Blp3/2cos3 / +-i-~B~Zp2sin20
+O(p 5'2)
where the free surface is = Blrl 3/2 + 0(r/2).
3.2 Axisymmetric flow A similar local coordinate system is taken. Let ~ = z, t / = r - d and define p, 0 so that = p sin 0 and t / = p cos 0. The techniques for obtaining local solutions for axisymmetric flow are less well-known. The method used to derive the following formula is similar to that described by Crank and Furzeland [15]. A first approximation for ¢ is
~b =qd
15B2l P 2 sin20 } --p sin0 + Blp3/2cos30/2 +~-~ _ ½qp2sin 20 + 0@ 5/2)
where the free boundary is = B it/3/2. 3.3
Application to numerical solution
The above formulae were considered with respect to their incorporation in the numerical solutions, but finally they were not included for two reasons. (i) The formulae show that the singularity is very mild and so a straightforward finite element solution will not be seriously upset by the omission of special treatment. (gi) The series are very slowly convergent and so need several terms to make them usable. 4. C A L C U L A T I O N
OF THE CONSTANT
q,
The problem, as so far described, can be solved for a range of values of q,.. However, it is known (Birkhoff and Zarantello[1]) that for fixed L and d only one value of q,. gives an analytic solution. Recall that 1_2 = ~qc 1~2 p + ~q
everywhere in the fluid. Hence we can calculate the pressure p at any point in the flow in
The numerical solution of planar and axisymmetric cavitational flow prqblems
63
terms of the constant qc which has been prescribed and the function ~b which has been calculated. If the value of qc is taken to be too low then there are areas of fluid where p is negative. This is obviously wrong since this would lead to the creation of further cavities. Thus the algorithm for determining qc is as follows: (i) Solve the problem for a value of qc which is known to be too low. (ii) Calculate the average pressure in each element {Pi}. (iii) If any of the {pi} are negative, increase q~ and repeat step (ii), otherwise stop. Thus we aim to find the lowest value of q~ for which the pressure is positive everywhere inside the flow.
5. R E S U L T S
5.1 P l a n a r f l o w
Referring to Fig. 2, the results are presented for A F = 1, B C = d, B G = L , A B = a and D G = b. In the original problem A B is infinite, but obviously for a numerical solution it must be given a finite value, A B = a. Experiments were carried out to find a suitable value of a so that this value had no significant effect on the flow. The main object of the computations was to find the variation of qc and b with L and d. The results are shown in Table 1. These results Show the correct qualitative behaviour in terms of increase or decrease of qc or b with respect to L or d. There are very few other results with which to make comparison. Mogel and Street[5] have solved the planar problem in terms of velocity components by finite difference methods. They quote results for the case d = 0.1 which are higher than those given here, and also show much greater dependence on the value of a than we found in the present investigation. The results for the case d = 0.1, L = 0.5 are compared in Table 2, where M&S represents the results of Mogel and Street and JMA represents the present investigation. There is one difference between the models used: Mogel and Street use the boundary condition ~b l a x = 0 on A F instead of ~b = y. However, on experiment, we found that this made very little difference to the values of qc and b, and did not account for the differences between the two sets of results. The only analytic results available are for a channel of infinite width, and are given in Birkhoff and Zarantello[1]. The formulae are given in terms of elliptic integrals of the first and second kind, but it is possible to evaluate these numerically. If we consider the case where L / d = 5, and study the effect of increasing the ratio of channel width to plate width, the results are seen to tend to the analytic solution, as shown in Table 3.
Table 1. Variation of qc and b with L and d d
L
a
qc
b
cr
0.05 0.1 0.1 0.2 0.2
0.25 0.5 1.0 0.5 1.0
0.6 1.0 1.4 1.0 1.4
1.370 1.429 1.346 1.744 1.696
0.1020 0.2013 0.2258 0.3346 0.3862
0.877 1.042 0.812 2.042 1.876
Table 2. Comparison of different methods Method JMA M&S M&S M&S
a a a a
= = = =
1 0.5 1 2
b
a
0.2013 0.242 0.245 0.245
1.04 2.11 1.44 1.22
64
J. M. AITCHISON Table 3. The effect of increasing channel width d
L
0.2 1.0 0.1 0.5 0.05 0.25 analytic
1/d
cr
5 10 20 oc
1.876 1.042 0.877 0.892
Table 4. Variation of q,. and b with d and L d
L
a
qc
b
a
0.125 O.125 0.25 0.25
0.5 1.0 0.5 1.0
1.5 2.0 1.5 2.0
2.258 2. ! 99 2.559 2.488
0.2033 0.2502 0.3579 0.4121
0.275 0.209 0.637 0.548
5.2 Axisymmetricflow The main results for this problem are shown in Table 4. Again the results show the correct qualitative behaviour. For this problem there are no analytic solutions available but there are some previous numerical solutions. Fox and Sankar[9] consider the case d = 0.125, L = 1.0, and produce the value of qc= 2.11 which is lower than the value obtained here. Brennen[8] has graphs showing his calculated dependence of a on L and d. By interpolation from his graphs we obtain the following
d=0.125,
L=I.0,
cr=0.28-0.29
d=0.!25,
L=0.5,
a=0.34-0.35.
Both of these are higher than those obtained in the present work. Calculations were performed using a variety of grids. The results quoted are for the finest grid used which seems to cope well with the singularity at the separation point. 6. C O N C L U S I O N
The method described above works well for both planar and axisymmetric problems. There are none of the difficulties of convergence of the free boundary position found in some other methods. In so far as comparison with other results is possible, the results here compare reasonably well with those from other methods. The method of variable elements is not limited to this problem and could be applied to any flow problem which can be described by a suitable minimization principle. Acknowledgements--Most of the work described here was performed while the author was the holder of a joint research fellowship with the Rutherford Laboratory, Chilton, and Pembroke College Oxford. The use of the computing facilities at the Rutherford Laboratory is gratefully acknowledged.
REFERENCES 1. G. Birkhoff and E. H. Zarantello, Jets, Wakes and Cavities. Academic Press, New York (1957). 2. D. Gilbarg, Jets and Cavities. In Handbuch der Physik (Edited by S. Flugge), Vol. IX. Springer, Berlin (1960). 3. T. Y. Wu, Inviscid Cavity and Wake Flows. In Basic Developments in Fluid Dynamics (Edited by M. Holt), Vol. 1. Academic Press, New York (1968). 4. D. Riabouchinsky, On steady fluid motion with free surfaces. Proc. London Math. Soc. 19, 206-215 (1919). 5. l'. R. Mogel and R. L. Street, A numerical method for steady state cavity flows. J. Ship Res. 18, 22-31 (1974). 6. D. M. Young, L. D. Gates, R. J. Arms and D. F. Elizer, The computation of an axisymmetric free boundary problem on NORC, Part I. NPG Report 1413 US Naval Proving Ground, Dahlgren, Virginia (1955). 7. R. J. Arms and L. D. Gates, The computation of an axisymmetric free boundary problem on NORC, Part II. NPG Report 1533. US Naval Proving Ground, Dahlgren, Virginia (1957). 8. C. Brennen, A numerical solution of axisymmetric cavity flows. J. Fluid Mech. 37, 671-688 (1969). 9. L. Fox and R. Sankar, The Regula-Falsi method for free boundary problems. J. Inst. Math. Applic. 12, 49-54 (1973).
The numerical solution of planar and axisymmetric cavitational flow problems
65
10. R. Courant and D. Hilbert, Methods of Mathematical Physics, Vol. I. Interscience, New York (1953). 11. G. Strang and G. J. Fix, An Analysis of the Finite Element Method. Prentice Hall, Englewood Cliffs, New Jersey (1973). 12. O. C. Zienkiewicz, The Finite Element Method. McGraw-Hill, London (1977). 13. J. M. Aitchison, The interactive generation of a 2-D finite element mesh. Report RL-78-034, Rutherford Laboratory, Chilton (1978). 14. R. Fletcher, Fortran subroutines for minimization by quasi-Newton methods. AERE Report 7125, Harwell (1972). 15. J. Crank and R. M. Furzeland, The treatment of boundary singularities in axially symmetric problems containing discs. J. Inst. Math. Applic. 20, 355-370 (1977).
0045-7930/84 $3.00+ .00 Pergamon Press Ltd.
THE NUMERICAL SOLUTION OF PLANAR AND AXISYMMETRIC CAVITATIONAL FLOW PROBLEMS JOYCE M. AITCHISON Oxford University Computing Laboratory, 19 Parks Road, Oxford OXl 3PL, England (Received 30 October 1981; in revisedform 6 September 1982) Abstract--This paper describes a method for the numerical solution of planar and axisymmetric Riabouchinsky cavity flows. The problem is formulated as the minimization of a functional over a variable region, which is solved numerically by the method of variable finite elements. The convergence difficulties experienced by some other authors do not appear here, and the method is found to be successful for both planar and axisymmetric problems.
1. INTRODUCTION 1.1 The model In this paper we consider the cavitational flow of an incompressible inviscid fluid past a plate or disk in a channel of finite width and infinite length. The fluid flows from left to right (Fig. la) and has a uniform velocity q0 far upstream. It meets a plate CK placed symmetrically in the channel at right angles to the flow. The basic problem is to find the flow pattern behind the plate. Immediately behind the plate a cavity is formed containing air at a constant pressure: our primary objective is to find the shape and size of this cavity. It is reasonably clear from experiment what sort of shape the cavity should have at its left hand end. However, it is not at all clear (from experiment or theory) what form the cavity should take at the right hand end. This phenomenum, known as the closure condition, is still much discussed. Various models have been put forward (see, for example, Birkoff and Zarantello [1], Gilbarg[2] or Wu [3]). They include an infinite cavity length and a re-entrant jet. Here we consider the model due to Riabouchinsky[4], in which we suppose that an image plate C'K" can be placed in the flow at some point downstream from the original plate and that the flow geometry will be symmetric about the line midway between CK and C'K" (Fig. lb). Since the flow is also symmetric about the axis of the channel, we need only consider the quarter of the flow region shown in Fig. 2. The lines AG and EG are axes of symmetry. Previous numerical solutions using this model include those by Mogel and Street [5] for planar flow and by Young et al. [6], Arms and Gates [7], Brennen [8] and Fox and Sankar [9] for the axisymmetric problem. We will first describe the case of planar flow and then give the corresponding equations for axisymmetric flow in Section 1.4. 1.2 The differential equations for planar flow In forming the governing differential equations of the flow we will make use of Bernoulli's equation
p + ½q2= constant where p is the pressure and q the velocity of the fluid. Along the free surface CD, p = 0 and so ½q2 = constant = ½qc2 (say). Thus p + ½q2= lqc2 everywhere. Inside the fluid p i> 0 and so q2 ~
56
J.M.
AITCHISON
(a)
I
I I
(bJ
Fig. 1. (a) A general cavity flow. (b) The Riabouchinsky model.
E
I
C
I I i
I
B.
L
,G
Fig. 2. The region of solution.
number, tr, which is defined by qc 2 - -
O-~.
qo2
q0 2
The stream function ~k(x, y ) satisfies the following differential equation, where I) is the flow region A B C D E F (Fig. 2)
~2~ a2~ dx---~+tgy--~=0 in~
The numerical solution of planar and axisymmetriccavitational flow problems
57
subject to: qoY
d~ =
on A F ( x - - ~ - o o )
= 0
on AB, BC onFE
~b = qoH
t~b ?x
0
on ED
and q;=O on CD, - -
-
-
qc
~n where n is the direction of the outward normal from Q. Note the double boundary condition on the free surface CD. The lengths H and d will be specified together with the upstream velocity q0, but there is some freedom in the specification of another parameter to fully define the problem. The natural choice from a physical point of view is to fix qc and let the half-length of the cavity, L, be unknown. This is equivalent to specifying the pressure far upstream. However, small changes in the pressure lead to large changes in the cavity length and so this formulation is numerically unstable. Hence we prefer to fix the length L and let q~ be an unknown constant to be determined as part of the solution. The problem can be made non-dimensional by putting x*=x/H,
y*=y/H,
~k*=~b/qoH
d* = d/H,
L* = L / H ,
qc* = qc/qo.
The equations become (dropping the *'s) in f~
3x z +~-~-yZ= 0
(1)
subject to =y
on A F
= 0
on AB, BC
= 1
on FE
t~-x= 0
on ED
ag,
(2)
and
O~_ On
on CD qc
where the region t2 is still that shown in Fig. 2, but with A F The cavitation number is now given by a = qJ - 1.
= 1.
(3)
The numerical problem is now to solve eqns (1) and (2) for ~b(x, y), qc and the free surface C A F Vol. 12, No. I - - E
58
J . M . AITCHISON
CD. In fact eqns (1) and (2) can be solved numerically for a range of values of qc. Only one of these corresponds to the true flow pattern. The technique for finding this particular value of q,. will be described in Section 4, for the present we will assume that q, is fixed. 1.3 The associated variational principle Let
+ \ O y / i d + - ~ q c 2 dxdy.
E[~,, ~] =
(4)
The function ~b(x,y) and the region f~ which satisfy eqns (1) and (2) can be found by minimizing the functional E with respect to both ~b and f~. f~ is allowed to vary by moving the line CD, where the point C is fixed and the point D moves along the line x = L. is subject to the constraints ~k = 1 on FE
(5)
on AF
=y ~b=O
o n A B , BC, CD.
Under these conditions it is easily shown (Courant and Hilbert[10], pp. 260-262) that E has a stationary value when • ~
t3x
= 0
in f~
- 0
on ED
and (V~b)2 = q,,2
on CD.
Thus eqns (1) and (2) are' fully satisfied. It can be further shown that this stationary value is in fact a minimum. This minimization problem will be solved numerically by the method of variable finite elements. 1.4 Axisymmetric flow For the case of axisymmetric flow we take coordinates (z, r) replacing (x, y). The non-dimensional equations corresponding to (1) and (2) are ~2~
1 O~O 02~/ r Or + ~ = 0
Or2
in f~
(6)
subject to = r 2 on AF ~b=O
a~O Oz
o n A B , BC
= 1
on FE
0
on ED
and #/=0 ~ r On 1
on CD. qc
(7)
The numerical solution of planar and axisymmetrie eavitational flow problems
59
In this case qc 2
a=~
1.
(8)
The solution to eqns (6) and (7) corresponds to a stationary value of E[~h,D]=
f ~ f {12 r LFfto~h~ \ t o r ] 2 + (to@)2] ~z + ~,q c 2r } drdz.
is allowed to vary by moving
(9)
CD as before, and ~h is constrained by
FE ~ = r 2 on AF ~O= 0 on AB, BC, CD. = 1
on
(1o)
Again it is easily shown that E has a stationary value when tO (1 to~b'~
to (1tomb'] 1 tomb
r toz
in f~
- 0
on
EL)
on
CD;
and ~ (\~r j
\ t o z ] j = qc2
and so eqns (6) and (7) are fully satisfied. 2. T H E M E T H O D O F V A R I A B L E F I N I T E E L E M E N T S
The finite element method for the numerical solution of variational principles over fixed domains is well-known (see, for example, Strang and Fix[ll] or Zienkiewicz[12]). However, many problems in fluid dynamics can be expressed as a variational principle over a variable domain, typified by the example given here. The numerical solution of such problems can be obtained by the method of variable finite elements. The region ~ is covered by a net of elements: we use triangles here, although the arguments are equally applicable to quadrilaterals (but not rectangles). The function ~, is approximated, as usual, by piecewise polynomials, the coefficients of which are chosen to minimize the functional E[~, ~)]. The region 1) is varied by distorting the geometry of the elements. This is achieved by allowing the coordinates of the vertices of the triangles to vary. However, care must be taken that we do not allow too many degrees of freedom in this geometric variation. If we allow the (x, y) coordinates of all the vertices to be free then we solve not only for the position of CD but also for the optimal position of nodes within ft. Under certain conditions this might be worthwhile, but here the problem would become too large. The approach used here is that all x coordinates are kept constant and only a limited number of y coordinates are allowed to vary. The linkage structure between elements remains fixed.
2.1 Specification of the geometry The initial triangulation is performed on a region where the line CD has an arbitrary, simple form, for example the line y = d (Fig. 3). For the present work this region was triangulated using an interactive graphics program on the PRIME 400 computer at the Rutherford Laboratory. This program is described in detail in Aitchison[13]. The triangulation can be described by the following data:
60
J. M. AITCHISON F
E
A
B
Fig. 3. A basic grid.
NP:. number of vertices. NE: number of elements. X, Y[I:NP]: arrays giving the (x, y) coordinates of each vertex. NOP[I:NE, 1:3]: integer array giving the three vertices of each element. Each subsequent position of the free boundary is specified by a set of ordinates {h~}7=~. Figure 4 shows the example n = 3. The triangulation for this region is generated from the basic triangulation as follows. (i) The array NOP giving the linkage structure remains the same. (ii) The array X giving the x coordinates of the vertices remains the same. (iii) The entries in the array Y which refer to points to the left of the plate remain the same. (iv) The other entries in the array Y are scaled linearly as follows. Consider the vertex (x0, Y0) in the basic grid. Let h0 be the height of CD at x = x0 in the basic grid. Let h~ be the corresponding height in the current grid. Then the vertex (x0, Y0) becomes (x0, Y0 where Yl=l
( 1 - - hi) --(l--y0). (1 -- h0)
The most detailed grid used in this investigation is shown in Fig. 5.
.B
Fig. 4. A distorted grid.
LC"~
Fig. 5. A typical grid.
The numericalsolution of planar and axisymmetriecavitationalflowproblems
61
2.2 The discrete minimization problem The function ff is approximated by piecewise linear functions for the case of planar flow, nodes being taken at the vertices of the triangles. However, for the axisymmetric problem quadratic functions are more appropriate since ~b behaves like r 2 away from the plate. Nodes are taken at the vertices and at the midpoints of the sides of the triangles. An expression for the discrete approximation of the functional E[ff, f~] can now be written in terms of the geometric v a r i a b l e s {hi}~=1 and the nodal values {ui}~=l which approximate to ~ at those points where it is not directly prescribed by boundary conditions. The discrete functional Ed is of the form /
Ea(u, y) = ~ u A (y)u + b'(y) • u + c(y)
where the matrix A is symmetric and positive definite. For stationary values we require
~E~ Ou~
=0
i=l,...,s
(11)
i = 1. . . . . n.
(12)
and
~E~
--=0 Oh~ Equations (11) give
d(h)u= -b(h)
(13)
which is linear in u, but eqns (12) are difficult to derive analytically since they involve the differentiation of a complicated function of h. An alternative method is as follows: for given h, eliminate u from Ed by solving eqn (13). Then Ed(U, h) -- E*(h) = lU'A (h)u + b'(h) • u + c(h)
= ½b'(h) • u + c(h)
where u = - A (h)- ~b(h). We now find the value of h which minimizes E*(h) by a quasi-Newton method due to Fletcher[14]. This optimization problem is much smaller than the original one since the majority of the unknowns have now been eliminated. 2.3 Efficient solution of the linear equations The set of linear eqns (13) has to be solved for many different values ofh. Various steps are taken to ensure that this is done as efficiently as possible. The nodes on the finite element grid are numbered so that ni < nj only if xi <<,xj. Thus the equations relating to those nodes to the left of BC appear in the top part of the matrix A. Let these be the first n~ rows. The elements of A in these rows do not depend on h and so do not have to be recalculated at each stage. The equations (13) are solved by Cholesky decomposition, A = L L r. Since A is symmetric, only its lower half is stored and the matrix L overwrites this. On the first iteration of the minimization process, the whole of A is calculated and decomposed. Subsequently only those rows beyond row n 1 are calculated and decomposed. Thus the first n~ rows remain in decomposed form. The matrix A is banded, but has variable bandwidth due to the irregularity of the nodes. It is stored in a one-dimensional array in variable bandwidth form.
62
J . M . AITCItISON 3. T H E S I N G U L A R I T Y
AT THE SEPARATION
POINT
In both the planar and axisymmetric problems there is a mild singularity in the function at the separation point C. In the following formulae q -~ q,.
3.1 Planar flow Take a local coordinate system (4, r/) centered at C. Thus =x
and
r/=y-d.
Define local polar coordinates (p, 0) so that =psin0
and
q=pcos0.
Note that this is not the usual definition of 0. A local series solution for ff can be found by standard techniques. A first approximation to the local solution for ~ is: ~ =q
- p sinO
+Blp3/2cos3 / +-i-~B~Zp2sin20
+O(p 5'2)
where the free surface is = Blrl 3/2 + 0(r/2).
3.2 Axisymmetric flow A similar local coordinate system is taken. Let ~ = z, t / = r - d and define p, 0 so that = p sin 0 and t / = p cos 0. The techniques for obtaining local solutions for axisymmetric flow are less well-known. The method used to derive the following formula is similar to that described by Crank and Furzeland [15]. A first approximation for ¢ is
~b =qd
15B2l P 2 sin20 } --p sin0 + Blp3/2cos30/2 +~-~ _ ½qp2sin 20 + 0@ 5/2)
where the free boundary is = B it/3/2. 3.3
Application to numerical solution
The above formulae were considered with respect to their incorporation in the numerical solutions, but finally they were not included for two reasons. (i) The formulae show that the singularity is very mild and so a straightforward finite element solution will not be seriously upset by the omission of special treatment. (gi) The series are very slowly convergent and so need several terms to make them usable. 4. C A L C U L A T I O N
OF THE CONSTANT
q,
The problem, as so far described, can be solved for a range of values of q,.. However, it is known (Birkhoff and Zarantello[1]) that for fixed L and d only one value of q,. gives an analytic solution. Recall that 1_2 = ~qc 1~2 p + ~q
everywhere in the fluid. Hence we can calculate the pressure p at any point in the flow in
The numerical solution of planar and axisymmetric cavitational flow prqblems
63
terms of the constant qc which has been prescribed and the function ~b which has been calculated. If the value of qc is taken to be too low then there are areas of fluid where p is negative. This is obviously wrong since this would lead to the creation of further cavities. Thus the algorithm for determining qc is as follows: (i) Solve the problem for a value of qc which is known to be too low. (ii) Calculate the average pressure in each element {Pi}. (iii) If any of the {pi} are negative, increase q~ and repeat step (ii), otherwise stop. Thus we aim to find the lowest value of q~ for which the pressure is positive everywhere inside the flow.
5. R E S U L T S
5.1 P l a n a r f l o w
Referring to Fig. 2, the results are presented for A F = 1, B C = d, B G = L , A B = a and D G = b. In the original problem A B is infinite, but obviously for a numerical solution it must be given a finite value, A B = a. Experiments were carried out to find a suitable value of a so that this value had no significant effect on the flow. The main object of the computations was to find the variation of qc and b with L and d. The results are shown in Table 1. These results Show the correct qualitative behaviour in terms of increase or decrease of qc or b with respect to L or d. There are very few other results with which to make comparison. Mogel and Street[5] have solved the planar problem in terms of velocity components by finite difference methods. They quote results for the case d = 0.1 which are higher than those given here, and also show much greater dependence on the value of a than we found in the present investigation. The results for the case d = 0.1, L = 0.5 are compared in Table 2, where M&S represents the results of Mogel and Street and JMA represents the present investigation. There is one difference between the models used: Mogel and Street use the boundary condition ~b l a x = 0 on A F instead of ~b = y. However, on experiment, we found that this made very little difference to the values of qc and b, and did not account for the differences between the two sets of results. The only analytic results available are for a channel of infinite width, and are given in Birkhoff and Zarantello[1]. The formulae are given in terms of elliptic integrals of the first and second kind, but it is possible to evaluate these numerically. If we consider the case where L / d = 5, and study the effect of increasing the ratio of channel width to plate width, the results are seen to tend to the analytic solution, as shown in Table 3.
Table 1. Variation of qc and b with L and d d
L
a
qc
b
cr
0.05 0.1 0.1 0.2 0.2
0.25 0.5 1.0 0.5 1.0
0.6 1.0 1.4 1.0 1.4
1.370 1.429 1.346 1.744 1.696
0.1020 0.2013 0.2258 0.3346 0.3862
0.877 1.042 0.812 2.042 1.876
Table 2. Comparison of different methods Method JMA M&S M&S M&S
a a a a
= = = =
1 0.5 1 2
b
a
0.2013 0.242 0.245 0.245
1.04 2.11 1.44 1.22
64
J. M. AITCHISON Table 3. The effect of increasing channel width d
L
0.2 1.0 0.1 0.5 0.05 0.25 analytic
1/d
cr
5 10 20 oc
1.876 1.042 0.877 0.892
Table 4. Variation of q,. and b with d and L d
L
a
qc
b
a
0.125 O.125 0.25 0.25
0.5 1.0 0.5 1.0
1.5 2.0 1.5 2.0
2.258 2. ! 99 2.559 2.488
0.2033 0.2502 0.3579 0.4121
0.275 0.209 0.637 0.548
5.2 Axisymmetricflow The main results for this problem are shown in Table 4. Again the results show the correct qualitative behaviour. For this problem there are no analytic solutions available but there are some previous numerical solutions. Fox and Sankar[9] consider the case d = 0.125, L = 1.0, and produce the value of qc= 2.11 which is lower than the value obtained here. Brennen[8] has graphs showing his calculated dependence of a on L and d. By interpolation from his graphs we obtain the following
d=0.125,
L=I.0,
cr=0.28-0.29
d=0.!25,
L=0.5,
a=0.34-0.35.
Both of these are higher than those obtained in the present work. Calculations were performed using a variety of grids. The results quoted are for the finest grid used which seems to cope well with the singularity at the separation point. 6. C O N C L U S I O N
The method described above works well for both planar and axisymmetric problems. There are none of the difficulties of convergence of the free boundary position found in some other methods. In so far as comparison with other results is possible, the results here compare reasonably well with those from other methods. The method of variable elements is not limited to this problem and could be applied to any flow problem which can be described by a suitable minimization principle. Acknowledgements--Most of the work described here was performed while the author was the holder of a joint research fellowship with the Rutherford Laboratory, Chilton, and Pembroke College Oxford. The use of the computing facilities at the Rutherford Laboratory is gratefully acknowledged.
REFERENCES 1. G. Birkhoff and E. H. Zarantello, Jets, Wakes and Cavities. Academic Press, New York (1957). 2. D. Gilbarg, Jets and Cavities. In Handbuch der Physik (Edited by S. Flugge), Vol. IX. Springer, Berlin (1960). 3. T. Y. Wu, Inviscid Cavity and Wake Flows. In Basic Developments in Fluid Dynamics (Edited by M. Holt), Vol. 1. Academic Press, New York (1968). 4. D. Riabouchinsky, On steady fluid motion with free surfaces. Proc. London Math. Soc. 19, 206-215 (1919). 5. l'. R. Mogel and R. L. Street, A numerical method for steady state cavity flows. J. Ship Res. 18, 22-31 (1974). 6. D. M. Young, L. D. Gates, R. J. Arms and D. F. Elizer, The computation of an axisymmetric free boundary problem on NORC, Part I. NPG Report 1413 US Naval Proving Ground, Dahlgren, Virginia (1955). 7. R. J. Arms and L. D. Gates, The computation of an axisymmetric free boundary problem on NORC, Part II. NPG Report 1533. US Naval Proving Ground, Dahlgren, Virginia (1957). 8. C. Brennen, A numerical solution of axisymmetric cavity flows. J. Fluid Mech. 37, 671-688 (1969). 9. L. Fox and R. Sankar, The Regula-Falsi method for free boundary problems. J. Inst. Math. Applic. 12, 49-54 (1973).
The numerical solution of planar and axisymmetric cavitational flow problems
65
10. R. Courant and D. Hilbert, Methods of Mathematical Physics, Vol. I. Interscience, New York (1953). 11. G. Strang and G. J. Fix, An Analysis of the Finite Element Method. Prentice Hall, Englewood Cliffs, New Jersey (1973). 12. O. C. Zienkiewicz, The Finite Element Method. McGraw-Hill, London (1977). 13. J. M. Aitchison, The interactive generation of a 2-D finite element mesh. Report RL-78-034, Rutherford Laboratory, Chilton (1978). 14. R. Fletcher, Fortran subroutines for minimization by quasi-Newton methods. AERE Report 7125, Harwell (1972). 15. J. Crank and R. M. Furzeland, The treatment of boundary singularities in axially symmetric problems containing discs. J. Inst. Math. Applic. 20, 355-370 (1977).