- Email: [email protected]
Finite difference solution to the flow between two rotating spheres
Computers& FluidsVol. 12, No. 2, pp. 77-92, 1984 Printed in Great Britain.
0045-7930/84 $3.00+.00 © 1984PergamonPress Ltd,
FINITE DIFFERENCE SOLUTION TO THE FLOW BETWEEN TWO ROTATING SPHERES S. C. R. DENNIS Department of Applied Mathematics, University of Western Ontario, London, Canada
and L. QUARTAPELLE Istituto di Fisica, Politecnico di Milano, Milano, Italy (Received 5 January 1983; in revised form 4 May 1983) AImraet--This paper describes a numerical method for calculating incompressible viscous flows between two concentric rotating spheres. The dependent variables describing the axisymmetric flow field are the azimuthal components of the vorticity, of the velocity vector potential and of the velocity. The coupled set of governing partial differential equations is written as a system of strictly second-order equations by introducing vorticity conditions of an integral character in a meridional plane. Such conditions generalize the one-dimensional integral conditions employed by Dennis and Singh to calculate steady-state solutions of the same problem using Gegenbauer polynomials and finite differences. The basic equations are discretized in space and in time by means of the finite-difference method. A fourth-order accurate centred-difference approximation of the advection terms is employed and a nonlinearly implicit scheme for the discrete time integration is here considered. A general finite-difference algorithm for steady-state and time-dependent problems is obtained which has no relaxation parameter and makes extensive use of fast elliptic solvers. The numerical results obtained by the present method are found to be in good agreement with the literature and confirm the nonuniqueness of the steady-state solution in a narrow spherical gap at certain regimes. 1. I N T R O D U C T I O N
When the axisymmetric motion of a viscous incompressible fluid in a spherical annulus is described in terms of the scalar variables such as the vorticity, the stream function and the circumferential velocity, the no-slip conditions on the surface of the spheres gives rise to a double boundary specification for the stream function and its normal derivative. This means that the system of two second-order partial differential equations for the vorticity ~ and the stream function ¢ must be more properly considered as a single fourth-order equation for ~, and that ( must be regarded as a shorthand notation for the Laplacian of ~ (at least, as far as the radial variable is concerned). This fact is typical of the ~-¢ equations whenever no-slip boundary conditions are present and it was clearly recognized by Pearson, who based the first finite-difference calculation of the flow between two concentric rotating spheres [1, 2] on the fourth-order equation for the stream function (see also Ref. [3]). On the other hand, the strictly second-order character of the vorticity-stream function equations can actually be attained provided that conditions of an integral type for the vorticity are introduced in the formulation of the continuum problem. Vorticity conditions of this type have been employed since more than a decade by Dennis et al. in the solution of two-dimensional and axisymmetric problems (see [4] and the references therein). These authors have used spatial discretizations of mixed type employing the truncated series method to represent the dependence on one spatial variable and the fiv.ite-difference method for the second variable. Then, they have introduced one-dimensional integral conditions for the vorticity to supplement the finite-difference equations for the vorticity expansion coefficients. Recently, the integral character of the vorticity conditions in incompressible viscous flows in two [5, 6] and three [6] dimensions has been recognized and algorithms employing the two-dimensional integral conditions for plane problems have been introduced and numerically investigated [7, 8]. The aim of this work is to formulate a finite-difference method employing the axisymmetric integral conditions for the vorticity to calculate stationary as well as transient flows in a spherical annulus. By virtue of such vorticity conditions, the numerical solution of the governing equations is reduced to a cascade of elliptic problems. These are solved by means 77
78
s.c.R. DENNISand L. QUARTAPELLE
of fast cyclic-reduction algorithms which make the present approach particularly efficient. Furthermore, the present study presents some novelties with respect to existing numerical work for plane problems using vorticity integral conditions which make such an investigation worthwhile. Firstly, owing to the existence of a rotation axis and to the symmetry of the flow with respect to the equatorial plane, the vorticity has a prescribed (zero) value on a part of the boundary of the two-dimensional integration domain. This implies that the vorticity equation is supplemented by conditions of both boundary (local) and integral (nonlocal) type. Secondly, in the considered problem, the vorticity integral conditions are homogeneous since the rotation of the spheres affects the meridional flow only through a nonlinear term in the vorticity transport equation. Therefore, the new conditions control the distribution of the vorticity within the domain rather than the generation of the vorticity near the boundaries, as typically occurs in plane flows. Thirdly, the set of coupled equations governing the fluid motion comprises a third equation for the azimuthal velocity. A final advantage of the finite-difference discretization in both spatial directions is the possibility of evaluating the head and the pressure by a simple procedure of line integration similar to the one devised for plane problems. 2. BASIC EQUATIONS Let us consider the flow of a viscous incompressible fluid occupying the region between two concentric spheres which can rotate about a common axis through the centre. The flow is assumed to be symmetrical about the axis of rotation which is taken as the axis 0 = 0 of the spherical coordinates (r, 0, ~o) with the origin placed at the centre of the spheres, so that all quantities are independent of the azimuthal angle tp. It is then possible to express the vectors velocity u and vorticity ~ in terms of their azimuthal and meridional components by introducing the azimuthal component of the vector potential (of the velocity field) denoted here by ~Oand called stream function for simplicity. We, therefore, have u = ~u + V × (~$),
(2.1)
= ~( + V x (~u),
(2.2)
where ~ is the unit vector normal to a meridional plane ~b = const, and V is the twodimensional gradient operator in a meridional plane, i.e. V - (O/Or, r-~O/O0). Let Q~ and ~92 be the angular velocity of the inner and outer spheres of radius rl and r2, respectively. The angular velocity t20 = max {(2~, t22}and the radius r2 of the outer sphere are taken as the units of angular velocity and length. The units of velocity, stream function and vorticity are u' - rzQo, ~' ~ r2u' -- rz2Q0 and (' - u ' / r 2 = (20, respectively, and the Reynolds number is defined by R = r22~o/V, where v is the coefficient of kinematic viscosity of the fluid. The Navier-Stokes equations for the dependent variables (, ~ and u read, in dimensionless form, 1
2
0(
k E ~ +~-~ +f(~, ~O,u) = 0, -
E2~
= ~,
(2.3) (2.4)
1 2 Ou -RE u + -~ + g(u, ~) = 0,
(2.5)
where E2 -= V2 -
1
0'
(2.6)
Ou Ofl O0 Or'
(2.7)
r 2 sin 2
J(~t, fl) = &t Ofl - - Or O0
Finite difference solution to the flow between two rotating spheres
79
Jl(~t) = J(a, r sin 0),
(2.8)
f((, ~,, u) = 1j((, ~k) + [~bJl(O + (Jt(~k) - 2uJl(u)l/(r 2 sin 0),
(2.9a)
g(u, ~ ) = 1j(u, d/ ) + [tpJl(u) - uJl($ )l/(r2sin 0).
(2.9b)
r
r
By virtue of the symmetry of the flow with respect to the equatorial plane, the equations can be solved in the reduced domain
v = [r~2<<.r <~ l] x IO < O < 21, where rt2 =_rl/r 2 (Fig. 1). The boundary conditions are thus stated in the form ~=$=u=0,
Ou
0=0,
(2.10)
rr
( = ~h = ~-~ = O, 0 = ~ ,
a¢
~k = ~ r =0'
r=rlzandr=
u = --t20r sin 0, ~'~2
u=--t20r sin0,
(2.1 l)
l,
(2.12)
r = rl2,
(2.13)
r = 1,
(2.14)
where, in general, f2~ = f2t(t) and 02 = I22(t). Conditions (2.10) on the axis of rotation are required to rule out the singularity of the eqns (2.3)-(2.5) at 0 = 0. Conditions (2.11) impose the symmetry ofu and the antisymmetry of~ and ~b with respect to the equatorial line 0 = ~/2. The no-slip conditions on the surface of the rotating spheres translate into eqns (2.12) for the motion in the meridional plane and into eqns (2.13) and (2.14) for the motion in the azimuthal direction. The mathematical
C_~fG
Fig. 1. Sphericalannulus and reduced integration domain for problems with symmetrywith respect to the equatorial plane
80
S . C . R . DENNIS and L. QUARTAPELLE
statement of the problem is made complete by the initial conditions ((t = 0 ) = (0 and u (t = 0) = u0, where (0 and u0 are the prescribed initial fields of the (azimuthal) vorticity and velocity, respectively. 3. INTEGRAL CONDITIONS ON THE VORTICITY IN AXISYMMETRIC PROBLEMS
The double specification (2.12) for ~k at the spherical boundaries can be replaced by an equivalent conditioning for the vorticity if use is made of the vorticity integral conditions for axisymmetric incompressible flows described in Ref. [6]. Let V denote an axisymmetric three-dimensional domain with boundary S. Moreover, let v be the two-dimensional domain resulting from the section of V with a meridional plane and s the boundary of v. Assume that an axisymmetric vector velocity is prescribed on the entire S, i.e. Uls = b = b(r, 0, t). Then, the azimuthal component ( of the vorticity vector is subjected to the integral conditions [6]
;
dvw(tl = -
f [ dsw
l v,w ,l
br/- a- u
(3.1)
w
with E2r/=0
in
v,
(3.2)
where w = w (r, 0) = r sin 0 and the second integral is evaluated along the entire closed curve s. In eqn (3.1) a and b are boundary data which depend on the prescribed velocity b according to the following relationships
1;
a = a(s, t) = w
d s w n " b,
(3.3a)
o
b = b(s, t) = - 3 • b,
(3.3b)
where n and 3 are normal and tangential unit vectors. In the problem considered here, we have v = [rl2, 1] x [0, rr/2] and there are two "internal" boundaries 0 = 0 and 0 = n/2 where the vorticity is prescribed. Therefore, the integral conditions (3.1) to be considered are only those required to fix the value of the vorticity on the two spheres. On the other hand, the velocity prescribed at the boundary is such that n" b = 0 and 3" b = 0 on the entire boundary of v except on the equatorial line where 3 • b is not known. It follows that the functions r/ to be used in the integral conditions (3.1) assume the following boundary values: arbitrary on the surfaces of the spheres, zero on the equatorial line and zero on the rotation axis, to rule out the singularity of the operator E z. The set of linearly independent functions r/ to be used in eqn (3.1) is, therefore, characterized by the problem (3.4a)
E2r] = 0
r/50,
r=rl2
q=0,
0=0
and
r=l,
(3.4b)
and
0=~,
(3.4c)
and henceforth the symbol q will denote any function which is a solution of eqns (3.4). The resulting vorticity integral conditions have the homogeneous form
f
dvw~rl = 0.
(3.5)
By introducing the usual notation (~, fl) =
dvwotfl = f'~/2 d0 t I d r r 2 s i n O ~ f l d0
dq2
(3.6)
Finite difference solution to the flow between two rotating spheres
81
for the Hilbert scalar product, eqn (3.5) can be rewritten as (~, r/) = 0.
(3.7)
The introduction of the integral conditions for the vorticity enables us to rewrite the set of eqns (2.3)-(2.5) and conditions (2.10)-(2.14) as a system of three strictly second-order, nonlinear partial differential equations
E2~ + N +f(~, ~, u) = 0,
(~, ,~) = 0, ~=0,
8=0
(3.8b) 7~
and
8 = 2,
- E2~k = ~', Is = O;
(3.8c) (3.9a) (3.9b)
1 2 Ou - - ~ E u +~-~ + g ( u , ~O) = 0,
u=0,
(3.8a)
(3.10a)
0=0,
(3.10b)
o-~ = o, o = ~ ,
(3.lOci
0u
[21r sin 8, u = --O0
r = r12,
(3.1 0d)
f22 u=--/20 r sin0,
r = 1.
(3.10e)
We notice that each equation is provided with its own independent conditioning so that the equations can be solved in sequence, apart from the coupling due to the nonlinear advection terms f and g. 4. TIME DISCRETIZATION SCHEME Starting from the formulation of the continuum problem provided by eqns (3.8)-(3.10), time discretized approximations and numerical algorithms for the calculation of timedependent and steady-state flows can be derived. Let us discretize the problem in time by the finite-difference method assuming a two-level scheme which allows us to solve the equations for ~ and ~k independently from the equation for u. We write ( _ E 2 + e)~.+l + Rf(~,+,, d/,+ 1, u") = e~",
(4.1)
__ E 2 ~ n +l = ~. + 1
(4.2)
( _ E 2 + e)un+l + Rg(u n, $n+l) = eu",
(4.3)
where e = R / A t and t n+l -~ ( n + 1)At. For simplicity, we have not rewritten the conditions for ~" ÷ l, $ , +1 and u" +l. They are obtained from the conditions of the continuum problem by appending the superscript index n + 1 to the unknown variables in the case of homogeneous conditions and by replacing the nonhomogeneous conditions (3.10d) and
S. C. R. DENNISand L. QUARTAPELLE
82 (3.10e), respectively, by
u "+1 = - - r f2o
sinO,
r
=
rl2 ,
(4.3d)
and 12~+1 u" + l = - r s i n 0 , O0
r=l,
(4.3e)
where ~27+l ---Ol(t "+1) and ~¢'2~+ 1 ~ - O2(tn+l). The steady-state solution is considered to be reached when
for some norm[I. I[ and a small positive constant a. Due to the presence of the term f(~, +l, ~b"+1, u") in the vorticity equation, eqns (4.1) and (4.2) are nonlinearly coupled and an iterative method of solution must be considered. The solution (~, ~O)"+ ~ is obtained as the limit of the sequence (~, ~b)", m = 0 , 1, 2 . . . . . defined as follows. Start from (~, ~b)m=0 = (¢, ~O)".Then, from u", ~" and (~, qJ)% calculate ~m+ 1 and ~Om+ 1 as the solutions of the equations ( - E 2 + e){ " + ' = e{" - Rf({ m, ~Om, u"),
(4.4)
_ E2@,,+ l = ~m+,,
(4.5)
supplemented by the homogeneous integral and boundary conditions as before. The iterative process is terminated when [1(~, ~ / ) m + l _ _ (~, I//)m[] ~ E
for some norm 1[. I[ and a small positive constant ¢. By setting e = 0, in eqn (4.4), the iterative algorithms (4.4) and (4.5) can also be used to calculate the solution of the nonlinear steady-state equations provided the equation for u is properly inserted in the iteration cycle. The resulting iterative method is defined by the following sequence of three independent elliptic equations -- E2~ m +1 = __ Rf(~',
~k", u"),
__ E21~tm+ 1 ~_- ~ m + l
(4.6) (4.7)
_ E2u,"+l = - R g ( u m, I~m+l).
(4.8)
5. SPATIAL DISCRETIZATION BY CENTRED DIFFERENCES Let the domain v be discretized by a uniform mesh of points defined by O,=(i-l)AO, rj=(j-1)Ar+r,2,
i=1,2 ..... M+I, j=l,2
.... ,N+I,
(5.1) (5.2)
where AO = (rr/2)/M and Ar =_ (1 - rn)/N. We assume the classical second-order accurate centred-difference approximation of the
Finite difference solution to the flow between two rotating spheres
83
elliptic operator E 2 defined in eqn (2.6), namely,
Ed2o~ij_m(V 2 '
1 )
1
2
r 2 sin 2 0 ~t,.j = (rjAr)2[rj+½(gij+l -- otij) -- r]_½(otij-- o~ij_l)] 1
+ (rjA 0)2 sin 0/sin 0~+½(ge+ ij - ~j) - sin 0i_ ½(% - ~_ lj)]
~t~j (rj sin 0~)2'
(5.3)
where rj+_½= rj +__½Ar and 0i±½= 0g -t- ½AO. If N is taken equal to an integer power of 2, the discrete counterpart of the elliptic eqns (4. I)-(4.3) can be solved by means of fast elliptic solvers [9] appropriately modified to deal also with the operator - E 2 + e, occurring in time dependent problems. In the spatially discrete case, the manifold of the functions r/defined by eqns (3.4) is characterized by the following set of equations
Ed2rlij=O, 2<~i <~M, 2<~j <~N,
(5.4a)
qij¢O,
(5.4b)
r/ij = 0 ,
2<~i<~M, j = l , i=l,M+l,
N+I,
I~
(5.4c)
The number of distinct boundary points considered in the nonhomogeneous boundary condition (5.4b) is equal to 2(M - 1). Therefore, there are L = 2(M - 1) linearly independent functions r/iJ, which will be denoted by r/~ij, ~ = 1, 2 , . . . , L. The two-dimensional integral given in eqn (3.6) can be evaluated according to the standard discrete approximation M+I
(o:, fl)~ = ~ i=1
N+I
~ Hi jr/sin O,~ijfl,j,
j=l
where Hij is equal to 1, ½and ¼for internal, boundary and corner points, respectively. The Jacobian (2.7) is discretized using standard fourth-order accurate differences throughout. The advection terms are evaluated with a higher accuracy than the diffusion terms since the presence of R in the transport equation makes the advection truncation error more important than the diffusion one at the intermediate Reynolds numbers we are interested in. For internal points not adjacent to the boundary the five-point centred approximation
(~/c~x)k = 1 ( _ ~k+2 + 8~k+ ~- 8ak_ 1+ a~-2)//Ix is used, whereas for internal points adjacent to the boundary the five-point uncentred approximation
(t~o:/~X)k =2 = 1-~(2Ctk+ 3 -- 6gk +2 + 18Ctk+ 1-- 1Oak + 3~tk_1)/Ax is used (the boundary is located at k = 1). For points near to the equator, i.e. for i = M, centred differences are used by introducing image points below the equatorial plane and by exploiting the symmetry of the flow. The same applies to evaluate J and Ji at i = M + 1 which are required due to the Neumann boundary condition for u "+1 at 0 = ~/2. 6. D I S C R E T I Z E D
VORTICITY INTEGRAL
CONDITIONS
In the spatially discrete case the vorticity integral conditions (3.8b) provide L linear algebraic equations. The typical discretized vorticity equations to be solved in the above
84
S. C. R. DENNIS and L. QUARTAPELLE
schemes are therefore ( - E 2 + e)~ = h, (~, q , ) = 0,
(6,1a)
~ = 1 , 2 . . . . . L, 7~
= 0,
0 = 0 and 0 = ~ ,
(6.1b) (6.1c)
to be interpreted in the sense of the assumed spatial discretization. In eqn (6.1a) ( = (,+1 and h =e~"-Rf(~",d/",u"), with possibly e = 0 . The linear system of algebraic equations (6.1) can be solved by means of the following decomposition scheme [8]. The solution ( is sought in the form L
= (+
~ 2,~,,
(6.2)
.'t=l
where ~"is the solution of the problem ( - E 2 + e)(" = h, ~'= arbitrary,
r
r12 and r = 1,
(6.3b)
and
(6.3c)
=
~=0, 0 = 0
(6.3a)
O = ~, ~
and the functions ~,, ~ = 1, 2 , . . . , L, are the solutions of the homogeneous Helmholtz problems ( - E 2 + e)~ = 0, ~t ~" ~t~t',
r = /'12
and
(6.4a) r
= 1,
(6.4b)
7~
~,=0,
0=0
and
0=~.
(6.4c)
If the arbitrariness in eqn (6.3b) is exploited by taking ~ = 0, it =- {21,..., 2L} is the discrete approximation of the (unknown) value of the vorticity on the spheres. The unknown it is found by imposing that ( satisfies the integral conditions (6.1b). We obtain the linear system of algebraic equations Ait=e
(6.5)
where the matrix A and the vector e are defined, respectively, by A,,, = (~,, ~/,),
(6.6)
c~ = - (~', t/,).
(6.7)
By virtue of Green's theorem, it can be shown that matrix A is symmetric. The calculation of A and e directly through expressions (6.6) and (6.7) is computationally not convenient. However, the use of Green's theorem shows that A and c can be evaluated in a different manner using a method suggested by Glowinski and Pironneau in a paper on the biharmonic problem [5]. The extension of this method to the case of axisymmetric problems in spherical coordinates of interest here is briefly summarized in the Appendix.
Finite differencesolution to the flow between two rotating spheres 85 7. EVALUATION OF THE PRESSURE The pressure, like other flow variables, is independent of the azimuthal angle tp. Therefore, it can be calculated from the scalar fields (, ~ and u, by a procedure of line integration in a meridional plane, similar to the method used by Burggraf in the case of plane stationary flows [10]. Here, the (dimensionless) head q and pressure p = q - llul2 are obtained from the Navier-Stokes equations for the velocity vector u, written as
Ou
t3t+~ xu=-Vq+~
1 V2
(7.1)
u,
where u and ~ are defined in eqns (2.1) and (2.2), respectively. Since all the variables are independent of ~p, Vq can be expressed in terms of the component of eqn (7.1) in a meridional plane, i.e. Vq=
1V w
-w
~
+
xt~-~
w
V(w~)+
u,- V ( w u ) , w
(7.2)
where w = r sin 0, as previously. In the time discretized case, eqn (7.2) becomes
Vq"+'= --wlVLw/
,tt
"[-~n-~R) l X ~ -- (n+ I I~-V(W~Cn+I) "I- un+ I!w V(wun+
(7.3)
If the r.h.s, of eqn (7.3) is known, q"+~ can be obtained by means of a line integration calculation assuming that q"+~ = 0 for 0 = 0 and r = rts. In the spatially discretized case, the gradients in the r.h.s, of eqn (7.3) are evaluated by second-order accurate centred differences. Then the line integration is performed by means of the trapezoidal rule averaged in the two directions. Once q,+l has been determined, p , + l is evaluated through the relationship p, + J = q, + 1 _
[u" + 1is +
V(w~b"+ 1)
(7.4)
using again centred differences. 8. RESULTS AND DISCUSSION The results given in this section were obtained taking equal to 10 -4 both the convergence parameter E for the iterative solution of the nonlinear ( - ~ equations and the threshold parameter a for the achievement of the steady state. Furthermore, the convergence criteria are based on the infinite norm. The calculations have been performed using the singleprecision arithmetic of the Amdahl 470/V8 machine ( ~ seven significant figures) except for double-precision accumulation in the LLT factorization and substitutions for the symmetric matrix A. A quantity of interest is the torque acting on the spheres. A dimensionless coefficient M may be defined as in [4] to give M =2jo3f"/2ra(OU~r- u ) sins 0 dO' where the integral is to be evaluated on the sphere concerned. The quantity Ou/Or is approximated by the one-sided, three-point, second-order accurate formula Ou/Or = 1(- 3ul + 4u2 - u3)/Ar. In Table 1 we give the calculated values of M1 and ME when rE = 2rl, f21 = t20, 02 = 0 and t21 = 0, ~2 = t20 for different mesh sizes. At the steady state, Mt and ME must be equal and opposite. The coarse radial mesh (N = 8) is found to be inadequate in both problems. A finer radial resolution (N = 16) provides values of - M~ and M s which differ by < 19/0at R ~< 500 and by ~ 3% at R = 1000, for the case with the inner sphere rotating and the outer sphere at rest. The values of - M ~ at R = 500 and at R = 1000 are in a satisfactory agreement with CAF 12:2-B
86
S. C. R. DENNIS and L. QUARTAPELLE Table 1. Nondimensional torque for rotating spheres Mesh size Mox Nr
0 2=
0
f21 = 0
R
-M,
M2
-Mi
M2
32x8
100 500 1000
0.435 0.728 0.903
0.446 0.951 1.321
0.500 0.814 0.983
0.540 0.656 0.611
2 0 x 16
100 500 1000
0.445 0.770 1.039
0.444 0.770 1.007
0.497 0.749 0.962
0.517 0.781 0.928
32 x 16
100 500 1000
0.445 0.770 1.039
0.444 0.768 1.006
0.497 0.749 0.966
0.517 0.733 0.920
the values 0.738 and 0.978 obtained by Dennis and Singh [4]. In the second case with the inner sphere at rest and the outer sphere rotating, the accuracy is lower but the results are in substantial agreement with the literature and the finest mesh (M = 32, N = 16) was assumed to be adequate for R > 500 in both cases. In order to present the numerical results obtained by the proposed method, we report some contour plots for the variables ~, O, u, co = u/(r sin 0) (angular velocity), q andp. Since
Vx(~0)=
rsin0
&0
'
'
the lines 0 = const are not the stream lines of the meridional flow, but for the case 0 = 0. The stream lines can be derived from the "true" stream function variable ¢ ' - Or sin 0. Nevertheless, the contour plots of 0 can be compared qualitatively with previously published plots of the stream lines because the modulation factor r sin 0 is a monotonic and relatively smooth function of both independent variables. In any case, the function Or sin 0 at 0 = 45 ° for the cases R = 10, 50 and 100, with r2 = 2q, ~r~1 = 0, ~'~2 = ~r'~0, has been compared with the corresponding solutions obtained by Munson and Joseph by a perturbation method [11]. The maximum local difference is found to be never > 3~. The vorticity and stream function fields for the case R = 100, ~r'~1 = ~'~0, ~2 "~ 0 are shown in Fig. 2. It is found that the region near the inner rotating sphere is one of intense vorticity generation. The pattern of the meridional flow agrees in character with the one calculated by Greenspan [12, 13]. In Fig. 3 the steady solution at R = 10 for a case with both spheres rotating at the same rate but in opposite directions is reported. This solution has been obtained by m = 6 iterations of the purely steady-state algorithm defined by eqns (4.6)-(4.8). The contour plot o f p shows that the pressure assumes maximum values on both spheres at the equator. Figure 4 contains the contour plots of the stream function and the angular velocity co = u/(r sin 0) at the steady state for the case R = 500, ~71 = 0, Q2 = O0 .
Fig. 2. Vorticity and stream function for R = 100, t71 = D0, g22 = 0. The inner sphere rotates and the outer sphere is at rest.
Finite difference solution to the flow between two rotating spheres
87
Fig. 3. Vorticity, stream function, head and pressure for R = 10, fit = -t20, ~2 = t20. Both spheres rotate at the same rate but in opposite directions.
Fig. 4. Stream function and angular velocity co = u/(r sin 0) for R = 500, I2t = 0, I22= t20. The inner sphere is at rest and the outer sphere rotates. calculated by the transient algorithm with At = 1 (n = 56, m = 244). There is a very reasonable agreement between the details o f Fig. 4 and the corresponding solution by Dennis and Singh [4, Fig. 3]. C o n t o u r s o f constant ~,, u and to for the case R = 1000, t2~ = 0, 02 = t20 at times t = 10, 20, 30 and 80 are given in Fig. 5. These m a y be c o m p a r e d with the transient solution obtained by Pearson [1]. The c o n t o u r s o f the azimuthal velocity clearly show the f o r m a t i o n o f a region o f almost solid b o d y rotation outside a cylindrical sheath o f radius approximately equal to that o f the inner sphere. As an indication o f the local accuracy o f the present numerical solution, we notice that the m i n i m u m o f ~,r sin 0 in the recirculatory zone near to the e q u a t o r is equal to - 0 . 0 4 5 , very close to the value - 0 . 0 5 calculated by Dennis and Singh. The final test o f the m e t h o d is the following problem o f a fairly n a r r o w gap: r~ = 0.85r2, I2t = O0 and 02 = 0. This problem has been extensively investigated both experimentally [14,
88
S. C. R. DENNIS and L. QUARTAPELLE
Fig. 5. Stream function, velocity and angular velocity ~ = u/(r sin 0) for R = 1000, t2~ = 0, 02 = t20. T h e inner sphere is at rest and the outer sphere is started impulsively at t = 0. Transient solution at t = 10, 20, 30 and 80 (At = 0.5).
16] and numerically [15-17]. Three different values of the Reynolds number Rw = _ florl2/v = R(ri/r2) 2 based on the radius of the inner sphere have been considered: R w = 700, 900 and 1500. Figure 6(a) shows the steady-state solutions for R w = 700 and Rw = 900. They have been obtained by the transient algorithm starting from rest and using At = 0.25. The total number of time steps and of iterations are n = 501, rn = 1535 for Rw = 700 and n = 450, m = 1934 for Rw = 900. For instance, the calculation of the case Rw = 900 requires the solution of about 10,000 elliptic equations. In the present case the size of the mesh is M = 50, N = 8 and the entire calculation including the graphical representation of all variables requires less than 4'. For Rw = 1500, two different stationary solutions have been calculated and the plots of the stream function and pressure are shown in Figs. 6(b, c), respectively. The solution with no vortex cell is obtained starting from rest impulsively whereas the solution with two counter rotating cells is obtained from the steadystate solution for Rw = 900. There is an excellent agreement between these solutions and the corresponding solutions calculated by different numerical methods [15-17]. The dimensionless size ~;t/(r2 - r0 of the first vortex cell near the equator is equal to 1.2, very close to the experimental value 1.3 found by Wimmer [14] (2 is the dimensional size of the vortex defined in [14]). Therefore, the present work confirms the nonuniqueness of the flow for Rw = 1500. The proposed numerical method has also been employed to attempt to calculate the stable flow at Rw = 1500 with a single vortex cell detected experimentally by Wimmer [14]. The rotating sphere is set in motion slowly according to the law flu(t) = max (1, t/z), z being a conveniently large time, so as to allow the development of the single-vortex flow [14]. Using z = 200 with At = 0.15 or z = 500 with At = 0.5, the fluid motion has no vortex initially but the two vortices appear at t --, 150 or t ,-, 250, respectively. A more gentle start up procedure is probably required to obtain the third stable solution at Rw = 1500 for this problem.
Finite difference solution to the flow between two rotating spheres
89
\
;I
ig
!!~1;
Fig. 6(a). Stream function for Rw = 700 (upper) and R w = 900 (lower), D~ = D0, D2 = 0. T
r
v ~ r ~
x ~,
r
~T~&qM~UNPTTgN
!&
-C
Fig. 6(b). Stream function for R w = 1500, D~ = D 0 and D 2 = 0. Upper: steady-state solution starting impulsively from rest. Lower: steady-state solution starting from the solution for R w = 900.
"
S. C. R. DENNISand L. QUARTAPELLE
90
.wm~
~
9.
1
Fig. 6(c). Pressure for Rw= 1500, flj = O0 and f12 =0. Upper: steady-state solution starting impulsively from rest. Lower: steady-state solution starting from the solution for Rw = 900.
9. C O N C L U S I O N
This work has demonstrated the applicability of two-dimensional vorticity integral conditions to the calculation of rotationally symmetric flows in spherical regions. A parameter-free algorithm for solving steady as well as unsteady problems has been formulated which is characterized by the following distinctive features: the noniterative determination of the vorticity values on the spheres, an implicit scheme for the time integration and the use of fast algorithms to solve the cascade of elliptic equations. The proposed method has been applied to a few test problems and has been found rapid and accurate in all the computed cases. The general formulation and the numerical method described is applicable to a number of other problems involving rotating spheres such as, e.g. the determination of steady secondary streaming induced by periodic torsional oscillations of the spheres [18, 19]. REFERENCES I. C. E. Pearson, A numerical study of the time-dependent viscous flow between two rotating spheres. J. Fluid Mech. 18, 323 (1967). 2. C. E. Pearson, in Studies in Numerical Analysis I. SIAM, 65 0968). 3. D. D. Joseph, Stability o f Fluid Motions I. Springer, Berlin (1976). 4. S. C. R. Dennis and S. N. Singh, Calculation of the flow between two rotating spheres by the method of series truncation. J. Comp. Phys. 18, 297 0978). 5. R. Glowinski and O. Pironneau, Numerical methods for the first biharmonic equation and for the two-dimensional Stokes problem. S I A M Rev. 21, 167 (1979). 6. L. Quartapelle and F. Valz-Gris, Projection conditions on the vorticity in viscous incompressible flows. Int. J. Num. Meth. Fluids 1, 129 (1981). 7. L. Quartapelle, Vorticity conditioning in the computation of two-dimensional viscous flows. J. Comp. Phys. 40, 453 (1981). 8. L. Quartapelle and M. Napolitano, A. method for solving the factorized vorticity-stream function equations by finite elements. Int. J. Num. Meth. Fluids 0984). 9. P. N. Swartztrauber and R. Sweet, Efficient FORTRAN subprograms for the solution of elliptic partial differential equations. N C A R - T N / I A - 1 0 9 (1975).
Finite difference solution to the flow between two rotating spheres
91
10. O. R. Burggraf, Analytical and numerical studies of the structure of steady separated flows. J. Fluid Mech. 24, 113 (1966). 11. B. R. Munson and D. D. Joseph, Viscous incompressible flow between concentric rotating spheres. Part 1. Basic flow. J. Fluid Mech. 49, 289 (1971). 12. D. Greenspan, Numerical studies of steady, viscous, incompressible flow between two rotating spheres. Computers and Fluids 3, 69 (1975). 13. D. Schultz and D. Greenspan, Improved solution of steady, viscous, incompressible flow between the rotating spheres. Computers and Fluids 7, 175 (1979). 14. M. Wimmer, Experiments on a viscous fluid flow between concentric rotating spheres. J. Fluid Mech. 7g, 317 (1976). 15. J.-P. Bonnet et T. Alziary de Roquefort, Ecoulement entre deux sphrre concentriques en rotation. J. Meca. 15, 373 (1976). 16. E. Krause and F. Bartels, in Approximation Methods for Navier-Stokes Problems, R. Rautmann, Editor. Springer, Berlin (1980). 17. G. Schrauf, Branching of Navier-Stokes equations in a spherical gap. Proc. 8th Int. Conf. Num. Meth. Fluid Dyn. (Edited by E. Krause). Lecture Notes in Physics 170, 474 (1982). 18. B. R. Munson and R. W. Douglass, Viscous flow in oscillatory spherical annuli, Phys. Fluids 22, 205 (1979). 19. Z. Zapryanov and C. Christov, Numerical study of the viscous flow in oscillatory spherical annuli. Comp. Meth. Appl. Mech. Engng 22, 49 (1980). APPENDIX The evaluation of A and c directly from the defining expressions (6.6) and (6.7) presents two disadvantages: (i) the calculation of the integral for each matrix coefficient or vector component requires O(MN) arithmetic operations; (ii) the storage of the functions r&, ~ = 1. . . . . L, which can be determined once and for all at the beginning of the calculation, requires L M N words of computer memory. It is possible to circumvent both inconveniences by a method proposed by Glowinski and Pironneau [5] at the expense of doubling the number of elliptic equations to be solved. Their method is here adapted to the case of axisymmetric problems in spherical coordinates. Let us consider two vector fields A and B with the following spherical components
A = (0, O, A¢~= A(r, 0)), a = (0, O, B~ = B(r, 0)).
For the special case of these rotationally symmetric vector fields the vector analogue of Green's theorem in three dimensions assumes the simple form
f d v w ( A E 2 B - BE2A)=
fdswn • (AVB -
BVA),
(Ai)
where, as previously, w = w(r, 0) = r sin 0 and n is the normal unit vector. Instead of considering the functions ~, defined through eqns (5.4), let us introduce the functions W~ defined as follows W~ = arbitrary in v,
r=rl2
W~=6~,,
and
(A2a) r=l.
(A2b)
For simplicity, in this Appendix we do not indicate explicitly that all functions are prescribed to be zero at 0 = 0 and 0 = n/2; this will be assumed implicitly. Let us now consider the following elliptic problems ( - E 2 + e)(~ = 0,
(~=W~,
r=rl2
and
(A3a) r=l;
- - E 2 ~ = (~, ~b~=0,
r=q2
(A3b) (A4a)
and
r=l.
(A4b)
With the aid of Green's theorem (A1), it is a simple matter to show that
(A5)
Similarly, considering the two elliptic problems ( - E 2 + e)~"= h, (=arbitrary,
r=rl2
(A6a)
and
r=l;
(A6b) (A7a)
~=0,
r=rl2
and
r=l;
(A7b)
92
S . C . R . DENNIS and L. QUARTAPELLE
it is easily obtained that c~ --- - C , .o)
(AS)
If the arbitrariness of the functions W~ at the internal points is exploited by taking W~ = 0 at all the internal points, the matrix coefficients A~,, and the vector components c~ can be evaluated through eqns (A5) and (A8) in 0(1) arithmetic operations. Furthermore, due to the boundary conditions (A2b), the functions W, can be generated explicitly whenever needed and there is no storage requirement.
0045-7930/84 $3.00+.00 © 1984PergamonPress Ltd,
FINITE DIFFERENCE SOLUTION TO THE FLOW BETWEEN TWO ROTATING SPHERES S. C. R. DENNIS Department of Applied Mathematics, University of Western Ontario, London, Canada
and L. QUARTAPELLE Istituto di Fisica, Politecnico di Milano, Milano, Italy (Received 5 January 1983; in revised form 4 May 1983) AImraet--This paper describes a numerical method for calculating incompressible viscous flows between two concentric rotating spheres. The dependent variables describing the axisymmetric flow field are the azimuthal components of the vorticity, of the velocity vector potential and of the velocity. The coupled set of governing partial differential equations is written as a system of strictly second-order equations by introducing vorticity conditions of an integral character in a meridional plane. Such conditions generalize the one-dimensional integral conditions employed by Dennis and Singh to calculate steady-state solutions of the same problem using Gegenbauer polynomials and finite differences. The basic equations are discretized in space and in time by means of the finite-difference method. A fourth-order accurate centred-difference approximation of the advection terms is employed and a nonlinearly implicit scheme for the discrete time integration is here considered. A general finite-difference algorithm for steady-state and time-dependent problems is obtained which has no relaxation parameter and makes extensive use of fast elliptic solvers. The numerical results obtained by the present method are found to be in good agreement with the literature and confirm the nonuniqueness of the steady-state solution in a narrow spherical gap at certain regimes. 1. I N T R O D U C T I O N
When the axisymmetric motion of a viscous incompressible fluid in a spherical annulus is described in terms of the scalar variables such as the vorticity, the stream function and the circumferential velocity, the no-slip conditions on the surface of the spheres gives rise to a double boundary specification for the stream function and its normal derivative. This means that the system of two second-order partial differential equations for the vorticity ~ and the stream function ¢ must be more properly considered as a single fourth-order equation for ~, and that ( must be regarded as a shorthand notation for the Laplacian of ~ (at least, as far as the radial variable is concerned). This fact is typical of the ~-¢ equations whenever no-slip boundary conditions are present and it was clearly recognized by Pearson, who based the first finite-difference calculation of the flow between two concentric rotating spheres [1, 2] on the fourth-order equation for the stream function (see also Ref. [3]). On the other hand, the strictly second-order character of the vorticity-stream function equations can actually be attained provided that conditions of an integral type for the vorticity are introduced in the formulation of the continuum problem. Vorticity conditions of this type have been employed since more than a decade by Dennis et al. in the solution of two-dimensional and axisymmetric problems (see [4] and the references therein). These authors have used spatial discretizations of mixed type employing the truncated series method to represent the dependence on one spatial variable and the fiv.ite-difference method for the second variable. Then, they have introduced one-dimensional integral conditions for the vorticity to supplement the finite-difference equations for the vorticity expansion coefficients. Recently, the integral character of the vorticity conditions in incompressible viscous flows in two [5, 6] and three [6] dimensions has been recognized and algorithms employing the two-dimensional integral conditions for plane problems have been introduced and numerically investigated [7, 8]. The aim of this work is to formulate a finite-difference method employing the axisymmetric integral conditions for the vorticity to calculate stationary as well as transient flows in a spherical annulus. By virtue of such vorticity conditions, the numerical solution of the governing equations is reduced to a cascade of elliptic problems. These are solved by means 77
78
s.c.R. DENNISand L. QUARTAPELLE
of fast cyclic-reduction algorithms which make the present approach particularly efficient. Furthermore, the present study presents some novelties with respect to existing numerical work for plane problems using vorticity integral conditions which make such an investigation worthwhile. Firstly, owing to the existence of a rotation axis and to the symmetry of the flow with respect to the equatorial plane, the vorticity has a prescribed (zero) value on a part of the boundary of the two-dimensional integration domain. This implies that the vorticity equation is supplemented by conditions of both boundary (local) and integral (nonlocal) type. Secondly, in the considered problem, the vorticity integral conditions are homogeneous since the rotation of the spheres affects the meridional flow only through a nonlinear term in the vorticity transport equation. Therefore, the new conditions control the distribution of the vorticity within the domain rather than the generation of the vorticity near the boundaries, as typically occurs in plane flows. Thirdly, the set of coupled equations governing the fluid motion comprises a third equation for the azimuthal velocity. A final advantage of the finite-difference discretization in both spatial directions is the possibility of evaluating the head and the pressure by a simple procedure of line integration similar to the one devised for plane problems. 2. BASIC EQUATIONS Let us consider the flow of a viscous incompressible fluid occupying the region between two concentric spheres which can rotate about a common axis through the centre. The flow is assumed to be symmetrical about the axis of rotation which is taken as the axis 0 = 0 of the spherical coordinates (r, 0, ~o) with the origin placed at the centre of the spheres, so that all quantities are independent of the azimuthal angle tp. It is then possible to express the vectors velocity u and vorticity ~ in terms of their azimuthal and meridional components by introducing the azimuthal component of the vector potential (of the velocity field) denoted here by ~Oand called stream function for simplicity. We, therefore, have u = ~u + V × (~$),
(2.1)
= ~( + V x (~u),
(2.2)
where ~ is the unit vector normal to a meridional plane ~b = const, and V is the twodimensional gradient operator in a meridional plane, i.e. V - (O/Or, r-~O/O0). Let Q~ and ~92 be the angular velocity of the inner and outer spheres of radius rl and r2, respectively. The angular velocity t20 = max {(2~, t22}and the radius r2 of the outer sphere are taken as the units of angular velocity and length. The units of velocity, stream function and vorticity are u' - rzQo, ~' ~ r2u' -- rz2Q0 and (' - u ' / r 2 = (20, respectively, and the Reynolds number is defined by R = r22~o/V, where v is the coefficient of kinematic viscosity of the fluid. The Navier-Stokes equations for the dependent variables (, ~ and u read, in dimensionless form, 1
2
0(
k E ~ +~-~ +f(~, ~O,u) = 0, -
E2~
= ~,
(2.3) (2.4)
1 2 Ou -RE u + -~ + g(u, ~) = 0,
(2.5)
where E2 -= V2 -
1
0'
(2.6)
Ou Ofl O0 Or'
(2.7)
r 2 sin 2
J(~t, fl) = &t Ofl - - Or O0
Finite difference solution to the flow between two rotating spheres
79
Jl(~t) = J(a, r sin 0),
(2.8)
f((, ~,, u) = 1j((, ~k) + [~bJl(O + (Jt(~k) - 2uJl(u)l/(r 2 sin 0),
(2.9a)
g(u, ~ ) = 1j(u, d/ ) + [tpJl(u) - uJl($ )l/(r2sin 0).
(2.9b)
r
r
By virtue of the symmetry of the flow with respect to the equatorial plane, the equations can be solved in the reduced domain
v = [r~2<<.r <~ l] x IO < O < 21, where rt2 =_rl/r 2 (Fig. 1). The boundary conditions are thus stated in the form ~=$=u=0,
Ou
0=0,
(2.10)
rr
( = ~h = ~-~ = O, 0 = ~ ,
a¢
~k = ~ r =0'
r=rlzandr=
u = --t20r sin 0, ~'~2
u=--t20r sin0,
(2.1 l)
l,
(2.12)
r = rl2,
(2.13)
r = 1,
(2.14)
where, in general, f2~ = f2t(t) and 02 = I22(t). Conditions (2.10) on the axis of rotation are required to rule out the singularity of the eqns (2.3)-(2.5) at 0 = 0. Conditions (2.11) impose the symmetry ofu and the antisymmetry of~ and ~b with respect to the equatorial line 0 = ~/2. The no-slip conditions on the surface of the rotating spheres translate into eqns (2.12) for the motion in the meridional plane and into eqns (2.13) and (2.14) for the motion in the azimuthal direction. The mathematical
C_~fG
Fig. 1. Sphericalannulus and reduced integration domain for problems with symmetrywith respect to the equatorial plane
80
S . C . R . DENNIS and L. QUARTAPELLE
statement of the problem is made complete by the initial conditions ((t = 0 ) = (0 and u (t = 0) = u0, where (0 and u0 are the prescribed initial fields of the (azimuthal) vorticity and velocity, respectively. 3. INTEGRAL CONDITIONS ON THE VORTICITY IN AXISYMMETRIC PROBLEMS
The double specification (2.12) for ~k at the spherical boundaries can be replaced by an equivalent conditioning for the vorticity if use is made of the vorticity integral conditions for axisymmetric incompressible flows described in Ref. [6]. Let V denote an axisymmetric three-dimensional domain with boundary S. Moreover, let v be the two-dimensional domain resulting from the section of V with a meridional plane and s the boundary of v. Assume that an axisymmetric vector velocity is prescribed on the entire S, i.e. Uls = b = b(r, 0, t). Then, the azimuthal component ( of the vorticity vector is subjected to the integral conditions [6]
;
dvw(tl = -
f [ dsw
l v,w ,l
br/- a- u
(3.1)
w
with E2r/=0
in
v,
(3.2)
where w = w (r, 0) = r sin 0 and the second integral is evaluated along the entire closed curve s. In eqn (3.1) a and b are boundary data which depend on the prescribed velocity b according to the following relationships
1;
a = a(s, t) = w
d s w n " b,
(3.3a)
o
b = b(s, t) = - 3 • b,
(3.3b)
where n and 3 are normal and tangential unit vectors. In the problem considered here, we have v = [rl2, 1] x [0, rr/2] and there are two "internal" boundaries 0 = 0 and 0 = n/2 where the vorticity is prescribed. Therefore, the integral conditions (3.1) to be considered are only those required to fix the value of the vorticity on the two spheres. On the other hand, the velocity prescribed at the boundary is such that n" b = 0 and 3" b = 0 on the entire boundary of v except on the equatorial line where 3 • b is not known. It follows that the functions r/ to be used in the integral conditions (3.1) assume the following boundary values: arbitrary on the surfaces of the spheres, zero on the equatorial line and zero on the rotation axis, to rule out the singularity of the operator E z. The set of linearly independent functions r/ to be used in eqn (3.1) is, therefore, characterized by the problem (3.4a)
E2r] = 0
r/50,
r=rl2
q=0,
0=0
and
r=l,
(3.4b)
and
0=~,
(3.4c)
and henceforth the symbol q will denote any function which is a solution of eqns (3.4). The resulting vorticity integral conditions have the homogeneous form
f
dvw~rl = 0.
(3.5)
By introducing the usual notation (~, fl) =
dvwotfl = f'~/2 d0 t I d r r 2 s i n O ~ f l d0
dq2
(3.6)
Finite difference solution to the flow between two rotating spheres
81
for the Hilbert scalar product, eqn (3.5) can be rewritten as (~, r/) = 0.
(3.7)
The introduction of the integral conditions for the vorticity enables us to rewrite the set of eqns (2.3)-(2.5) and conditions (2.10)-(2.14) as a system of three strictly second-order, nonlinear partial differential equations
E2~ + N +f(~, ~, u) = 0,
(~, ,~) = 0, ~=0,
8=0
(3.8b) 7~
and
8 = 2,
- E2~k = ~', Is = O;
(3.8c) (3.9a) (3.9b)
1 2 Ou - - ~ E u +~-~ + g ( u , ~O) = 0,
u=0,
(3.8a)
(3.10a)
0=0,
(3.10b)
o-~ = o, o = ~ ,
(3.lOci
0u
[21r sin 8, u = --O0
r = r12,
(3.1 0d)
f22 u=--/20 r sin0,
r = 1.
(3.10e)
We notice that each equation is provided with its own independent conditioning so that the equations can be solved in sequence, apart from the coupling due to the nonlinear advection terms f and g. 4. TIME DISCRETIZATION SCHEME Starting from the formulation of the continuum problem provided by eqns (3.8)-(3.10), time discretized approximations and numerical algorithms for the calculation of timedependent and steady-state flows can be derived. Let us discretize the problem in time by the finite-difference method assuming a two-level scheme which allows us to solve the equations for ~ and ~k independently from the equation for u. We write ( _ E 2 + e)~.+l + Rf(~,+,, d/,+ 1, u") = e~",
(4.1)
__ E 2 ~ n +l = ~. + 1
(4.2)
( _ E 2 + e)un+l + Rg(u n, $n+l) = eu",
(4.3)
where e = R / A t and t n+l -~ ( n + 1)At. For simplicity, we have not rewritten the conditions for ~" ÷ l, $ , +1 and u" +l. They are obtained from the conditions of the continuum problem by appending the superscript index n + 1 to the unknown variables in the case of homogeneous conditions and by replacing the nonhomogeneous conditions (3.10d) and
S. C. R. DENNISand L. QUARTAPELLE
82 (3.10e), respectively, by
u "+1 = - - r f2o
sinO,
r
=
rl2 ,
(4.3d)
and 12~+1 u" + l = - r s i n 0 , O0
r=l,
(4.3e)
where ~27+l ---Ol(t "+1) and ~¢'2~+ 1 ~ - O2(tn+l). The steady-state solution is considered to be reached when
for some norm[I. I[ and a small positive constant a. Due to the presence of the term f(~, +l, ~b"+1, u") in the vorticity equation, eqns (4.1) and (4.2) are nonlinearly coupled and an iterative method of solution must be considered. The solution (~, ~O)"+ ~ is obtained as the limit of the sequence (~, ~b)", m = 0 , 1, 2 . . . . . defined as follows. Start from (~, ~b)m=0 = (¢, ~O)".Then, from u", ~" and (~, qJ)% calculate ~m+ 1 and ~Om+ 1 as the solutions of the equations ( - E 2 + e){ " + ' = e{" - Rf({ m, ~Om, u"),
(4.4)
_ E2@,,+ l = ~m+,,
(4.5)
supplemented by the homogeneous integral and boundary conditions as before. The iterative process is terminated when [1(~, ~ / ) m + l _ _ (~, I//)m[] ~ E
for some norm 1[. I[ and a small positive constant ¢. By setting e = 0, in eqn (4.4), the iterative algorithms (4.4) and (4.5) can also be used to calculate the solution of the nonlinear steady-state equations provided the equation for u is properly inserted in the iteration cycle. The resulting iterative method is defined by the following sequence of three independent elliptic equations -- E2~ m +1 = __ Rf(~',
~k", u"),
__ E21~tm+ 1 ~_- ~ m + l
(4.6) (4.7)
_ E2u,"+l = - R g ( u m, I~m+l).
(4.8)
5. SPATIAL DISCRETIZATION BY CENTRED DIFFERENCES Let the domain v be discretized by a uniform mesh of points defined by O,=(i-l)AO, rj=(j-1)Ar+r,2,
i=1,2 ..... M+I, j=l,2
.... ,N+I,
(5.1) (5.2)
where AO = (rr/2)/M and Ar =_ (1 - rn)/N. We assume the classical second-order accurate centred-difference approximation of the
Finite difference solution to the flow between two rotating spheres
83
elliptic operator E 2 defined in eqn (2.6), namely,
Ed2o~ij_m(V 2 '
1 )
1
2
r 2 sin 2 0 ~t,.j = (rjAr)2[rj+½(gij+l -- otij) -- r]_½(otij-- o~ij_l)] 1
+ (rjA 0)2 sin 0/sin 0~+½(ge+ ij - ~j) - sin 0i_ ½(% - ~_ lj)]
~t~j (rj sin 0~)2'
(5.3)
where rj+_½= rj +__½Ar and 0i±½= 0g -t- ½AO. If N is taken equal to an integer power of 2, the discrete counterpart of the elliptic eqns (4. I)-(4.3) can be solved by means of fast elliptic solvers [9] appropriately modified to deal also with the operator - E 2 + e, occurring in time dependent problems. In the spatially discrete case, the manifold of the functions r/defined by eqns (3.4) is characterized by the following set of equations
Ed2rlij=O, 2<~i <~M, 2<~j <~N,
(5.4a)
qij¢O,
(5.4b)
r/ij = 0 ,
2<~i<~M, j = l , i=l,M+l,
N+I,
I~
(5.4c)
The number of distinct boundary points considered in the nonhomogeneous boundary condition (5.4b) is equal to 2(M - 1). Therefore, there are L = 2(M - 1) linearly independent functions r/iJ, which will be denoted by r/~ij, ~ = 1, 2 , . . . , L. The two-dimensional integral given in eqn (3.6) can be evaluated according to the standard discrete approximation M+I
(o:, fl)~ = ~ i=1
N+I
~ Hi jr/sin O,~ijfl,j,
j=l
where Hij is equal to 1, ½and ¼for internal, boundary and corner points, respectively. The Jacobian (2.7) is discretized using standard fourth-order accurate differences throughout. The advection terms are evaluated with a higher accuracy than the diffusion terms since the presence of R in the transport equation makes the advection truncation error more important than the diffusion one at the intermediate Reynolds numbers we are interested in. For internal points not adjacent to the boundary the five-point centred approximation
(~/c~x)k = 1 ( _ ~k+2 + 8~k+ ~- 8ak_ 1+ a~-2)//Ix is used, whereas for internal points adjacent to the boundary the five-point uncentred approximation
(t~o:/~X)k =2 = 1-~(2Ctk+ 3 -- 6gk +2 + 18Ctk+ 1-- 1Oak + 3~tk_1)/Ax is used (the boundary is located at k = 1). For points near to the equator, i.e. for i = M, centred differences are used by introducing image points below the equatorial plane and by exploiting the symmetry of the flow. The same applies to evaluate J and Ji at i = M + 1 which are required due to the Neumann boundary condition for u "+1 at 0 = ~/2. 6. D I S C R E T I Z E D
VORTICITY INTEGRAL
CONDITIONS
In the spatially discrete case the vorticity integral conditions (3.8b) provide L linear algebraic equations. The typical discretized vorticity equations to be solved in the above
84
S. C. R. DENNIS and L. QUARTAPELLE
schemes are therefore ( - E 2 + e)~ = h, (~, q , ) = 0,
(6,1a)
~ = 1 , 2 . . . . . L, 7~
= 0,
0 = 0 and 0 = ~ ,
(6.1b) (6.1c)
to be interpreted in the sense of the assumed spatial discretization. In eqn (6.1a) ( = (,+1 and h =e~"-Rf(~",d/",u"), with possibly e = 0 . The linear system of algebraic equations (6.1) can be solved by means of the following decomposition scheme [8]. The solution ( is sought in the form L
= (+
~ 2,~,,
(6.2)
.'t=l
where ~"is the solution of the problem ( - E 2 + e)(" = h, ~'= arbitrary,
r
r12 and r = 1,
(6.3b)
and
(6.3c)
=
~=0, 0 = 0
(6.3a)
O = ~, ~
and the functions ~,, ~ = 1, 2 , . . . , L, are the solutions of the homogeneous Helmholtz problems ( - E 2 + e)~ = 0, ~t ~" ~t~t',
r = /'12
and
(6.4a) r
= 1,
(6.4b)
7~
~,=0,
0=0
and
0=~.
(6.4c)
If the arbitrariness in eqn (6.3b) is exploited by taking ~ = 0, it =- {21,..., 2L} is the discrete approximation of the (unknown) value of the vorticity on the spheres. The unknown it is found by imposing that ( satisfies the integral conditions (6.1b). We obtain the linear system of algebraic equations Ait=e
(6.5)
where the matrix A and the vector e are defined, respectively, by A,,, = (~,, ~/,),
(6.6)
c~ = - (~', t/,).
(6.7)
By virtue of Green's theorem, it can be shown that matrix A is symmetric. The calculation of A and e directly through expressions (6.6) and (6.7) is computationally not convenient. However, the use of Green's theorem shows that A and c can be evaluated in a different manner using a method suggested by Glowinski and Pironneau in a paper on the biharmonic problem [5]. The extension of this method to the case of axisymmetric problems in spherical coordinates of interest here is briefly summarized in the Appendix.
Finite differencesolution to the flow between two rotating spheres 85 7. EVALUATION OF THE PRESSURE The pressure, like other flow variables, is independent of the azimuthal angle tp. Therefore, it can be calculated from the scalar fields (, ~ and u, by a procedure of line integration in a meridional plane, similar to the method used by Burggraf in the case of plane stationary flows [10]. Here, the (dimensionless) head q and pressure p = q - llul2 are obtained from the Navier-Stokes equations for the velocity vector u, written as
Ou
t3t+~ xu=-Vq+~
1 V2
(7.1)
u,
where u and ~ are defined in eqns (2.1) and (2.2), respectively. Since all the variables are independent of ~p, Vq can be expressed in terms of the component of eqn (7.1) in a meridional plane, i.e. Vq=
1V w
-w
~
+
xt~-~
w
V(w~)+
u,- V ( w u ) , w
(7.2)
where w = r sin 0, as previously. In the time discretized case, eqn (7.2) becomes
Vq"+'= --wlVLw/
,tt
"[-~n-~R) l X ~ -- (n+ I I~-V(W~Cn+I) "I- un+ I!w V(wun+
(7.3)
If the r.h.s, of eqn (7.3) is known, q"+~ can be obtained by means of a line integration calculation assuming that q"+~ = 0 for 0 = 0 and r = rts. In the spatially discretized case, the gradients in the r.h.s, of eqn (7.3) are evaluated by second-order accurate centred differences. Then the line integration is performed by means of the trapezoidal rule averaged in the two directions. Once q,+l has been determined, p , + l is evaluated through the relationship p, + J = q, + 1 _
[u" + 1is +
V(w~b"+ 1)
(7.4)
using again centred differences. 8. RESULTS AND DISCUSSION The results given in this section were obtained taking equal to 10 -4 both the convergence parameter E for the iterative solution of the nonlinear ( - ~ equations and the threshold parameter a for the achievement of the steady state. Furthermore, the convergence criteria are based on the infinite norm. The calculations have been performed using the singleprecision arithmetic of the Amdahl 470/V8 machine ( ~ seven significant figures) except for double-precision accumulation in the LLT factorization and substitutions for the symmetric matrix A. A quantity of interest is the torque acting on the spheres. A dimensionless coefficient M may be defined as in [4] to give M =2jo3f"/2ra(OU~r- u ) sins 0 dO' where the integral is to be evaluated on the sphere concerned. The quantity Ou/Or is approximated by the one-sided, three-point, second-order accurate formula Ou/Or = 1(- 3ul + 4u2 - u3)/Ar. In Table 1 we give the calculated values of M1 and ME when rE = 2rl, f21 = t20, 02 = 0 and t21 = 0, ~2 = t20 for different mesh sizes. At the steady state, Mt and ME must be equal and opposite. The coarse radial mesh (N = 8) is found to be inadequate in both problems. A finer radial resolution (N = 16) provides values of - M~ and M s which differ by < 19/0at R ~< 500 and by ~ 3% at R = 1000, for the case with the inner sphere rotating and the outer sphere at rest. The values of - M ~ at R = 500 and at R = 1000 are in a satisfactory agreement with CAF 12:2-B
86
S. C. R. DENNIS and L. QUARTAPELLE Table 1. Nondimensional torque for rotating spheres Mesh size Mox Nr
0 2=
0
f21 = 0
R
-M,
M2
-Mi
M2
32x8
100 500 1000
0.435 0.728 0.903
0.446 0.951 1.321
0.500 0.814 0.983
0.540 0.656 0.611
2 0 x 16
100 500 1000
0.445 0.770 1.039
0.444 0.770 1.007
0.497 0.749 0.962
0.517 0.781 0.928
32 x 16
100 500 1000
0.445 0.770 1.039
0.444 0.768 1.006
0.497 0.749 0.966
0.517 0.733 0.920
the values 0.738 and 0.978 obtained by Dennis and Singh [4]. In the second case with the inner sphere at rest and the outer sphere rotating, the accuracy is lower but the results are in substantial agreement with the literature and the finest mesh (M = 32, N = 16) was assumed to be adequate for R > 500 in both cases. In order to present the numerical results obtained by the proposed method, we report some contour plots for the variables ~, O, u, co = u/(r sin 0) (angular velocity), q andp. Since
Vx(~0)=
rsin0
&0
'
'
the lines 0 = const are not the stream lines of the meridional flow, but for the case 0 = 0. The stream lines can be derived from the "true" stream function variable ¢ ' - Or sin 0. Nevertheless, the contour plots of 0 can be compared qualitatively with previously published plots of the stream lines because the modulation factor r sin 0 is a monotonic and relatively smooth function of both independent variables. In any case, the function Or sin 0 at 0 = 45 ° for the cases R = 10, 50 and 100, with r2 = 2q, ~r~1 = 0, ~'~2 = ~r'~0, has been compared with the corresponding solutions obtained by Munson and Joseph by a perturbation method [11]. The maximum local difference is found to be never > 3~. The vorticity and stream function fields for the case R = 100, ~r'~1 = ~'~0, ~2 "~ 0 are shown in Fig. 2. It is found that the region near the inner rotating sphere is one of intense vorticity generation. The pattern of the meridional flow agrees in character with the one calculated by Greenspan [12, 13]. In Fig. 3 the steady solution at R = 10 for a case with both spheres rotating at the same rate but in opposite directions is reported. This solution has been obtained by m = 6 iterations of the purely steady-state algorithm defined by eqns (4.6)-(4.8). The contour plot o f p shows that the pressure assumes maximum values on both spheres at the equator. Figure 4 contains the contour plots of the stream function and the angular velocity co = u/(r sin 0) at the steady state for the case R = 500, ~71 = 0, Q2 = O0 .
Fig. 2. Vorticity and stream function for R = 100, t71 = D0, g22 = 0. The inner sphere rotates and the outer sphere is at rest.
Finite difference solution to the flow between two rotating spheres
87
Fig. 3. Vorticity, stream function, head and pressure for R = 10, fit = -t20, ~2 = t20. Both spheres rotate at the same rate but in opposite directions.
Fig. 4. Stream function and angular velocity co = u/(r sin 0) for R = 500, I2t = 0, I22= t20. The inner sphere is at rest and the outer sphere rotates. calculated by the transient algorithm with At = 1 (n = 56, m = 244). There is a very reasonable agreement between the details o f Fig. 4 and the corresponding solution by Dennis and Singh [4, Fig. 3]. C o n t o u r s o f constant ~,, u and to for the case R = 1000, t2~ = 0, 02 = t20 at times t = 10, 20, 30 and 80 are given in Fig. 5. These m a y be c o m p a r e d with the transient solution obtained by Pearson [1]. The c o n t o u r s o f the azimuthal velocity clearly show the f o r m a t i o n o f a region o f almost solid b o d y rotation outside a cylindrical sheath o f radius approximately equal to that o f the inner sphere. As an indication o f the local accuracy o f the present numerical solution, we notice that the m i n i m u m o f ~,r sin 0 in the recirculatory zone near to the e q u a t o r is equal to - 0 . 0 4 5 , very close to the value - 0 . 0 5 calculated by Dennis and Singh. The final test o f the m e t h o d is the following problem o f a fairly n a r r o w gap: r~ = 0.85r2, I2t = O0 and 02 = 0. This problem has been extensively investigated both experimentally [14,
88
S. C. R. DENNIS and L. QUARTAPELLE
Fig. 5. Stream function, velocity and angular velocity ~ = u/(r sin 0) for R = 1000, t2~ = 0, 02 = t20. T h e inner sphere is at rest and the outer sphere is started impulsively at t = 0. Transient solution at t = 10, 20, 30 and 80 (At = 0.5).
16] and numerically [15-17]. Three different values of the Reynolds number Rw = _ florl2/v = R(ri/r2) 2 based on the radius of the inner sphere have been considered: R w = 700, 900 and 1500. Figure 6(a) shows the steady-state solutions for R w = 700 and Rw = 900. They have been obtained by the transient algorithm starting from rest and using At = 0.25. The total number of time steps and of iterations are n = 501, rn = 1535 for Rw = 700 and n = 450, m = 1934 for Rw = 900. For instance, the calculation of the case Rw = 900 requires the solution of about 10,000 elliptic equations. In the present case the size of the mesh is M = 50, N = 8 and the entire calculation including the graphical representation of all variables requires less than 4'. For Rw = 1500, two different stationary solutions have been calculated and the plots of the stream function and pressure are shown in Figs. 6(b, c), respectively. The solution with no vortex cell is obtained starting from rest impulsively whereas the solution with two counter rotating cells is obtained from the steadystate solution for Rw = 900. There is an excellent agreement between these solutions and the corresponding solutions calculated by different numerical methods [15-17]. The dimensionless size ~;t/(r2 - r0 of the first vortex cell near the equator is equal to 1.2, very close to the experimental value 1.3 found by Wimmer [14] (2 is the dimensional size of the vortex defined in [14]). Therefore, the present work confirms the nonuniqueness of the flow for Rw = 1500. The proposed numerical method has also been employed to attempt to calculate the stable flow at Rw = 1500 with a single vortex cell detected experimentally by Wimmer [14]. The rotating sphere is set in motion slowly according to the law flu(t) = max (1, t/z), z being a conveniently large time, so as to allow the development of the single-vortex flow [14]. Using z = 200 with At = 0.15 or z = 500 with At = 0.5, the fluid motion has no vortex initially but the two vortices appear at t --, 150 or t ,-, 250, respectively. A more gentle start up procedure is probably required to obtain the third stable solution at Rw = 1500 for this problem.
Finite difference solution to the flow between two rotating spheres
89
\
;I
ig
!!~1;
Fig. 6(a). Stream function for Rw = 700 (upper) and R w = 900 (lower), D~ = D0, D2 = 0. T
r
v ~ r ~
x ~,
r
~T~&qM~UNPTTgN
!&
-C
Fig. 6(b). Stream function for R w = 1500, D~ = D 0 and D 2 = 0. Upper: steady-state solution starting impulsively from rest. Lower: steady-state solution starting from the solution for R w = 900.
"
S. C. R. DENNISand L. QUARTAPELLE
90
.wm~
~
9.
1
Fig. 6(c). Pressure for Rw= 1500, flj = O0 and f12 =0. Upper: steady-state solution starting impulsively from rest. Lower: steady-state solution starting from the solution for Rw = 900.
9. C O N C L U S I O N
This work has demonstrated the applicability of two-dimensional vorticity integral conditions to the calculation of rotationally symmetric flows in spherical regions. A parameter-free algorithm for solving steady as well as unsteady problems has been formulated which is characterized by the following distinctive features: the noniterative determination of the vorticity values on the spheres, an implicit scheme for the time integration and the use of fast algorithms to solve the cascade of elliptic equations. The proposed method has been applied to a few test problems and has been found rapid and accurate in all the computed cases. The general formulation and the numerical method described is applicable to a number of other problems involving rotating spheres such as, e.g. the determination of steady secondary streaming induced by periodic torsional oscillations of the spheres [18, 19]. REFERENCES I. C. E. Pearson, A numerical study of the time-dependent viscous flow between two rotating spheres. J. Fluid Mech. 18, 323 (1967). 2. C. E. Pearson, in Studies in Numerical Analysis I. SIAM, 65 0968). 3. D. D. Joseph, Stability o f Fluid Motions I. Springer, Berlin (1976). 4. S. C. R. Dennis and S. N. Singh, Calculation of the flow between two rotating spheres by the method of series truncation. J. Comp. Phys. 18, 297 0978). 5. R. Glowinski and O. Pironneau, Numerical methods for the first biharmonic equation and for the two-dimensional Stokes problem. S I A M Rev. 21, 167 (1979). 6. L. Quartapelle and F. Valz-Gris, Projection conditions on the vorticity in viscous incompressible flows. Int. J. Num. Meth. Fluids 1, 129 (1981). 7. L. Quartapelle, Vorticity conditioning in the computation of two-dimensional viscous flows. J. Comp. Phys. 40, 453 (1981). 8. L. Quartapelle and M. Napolitano, A. method for solving the factorized vorticity-stream function equations by finite elements. Int. J. Num. Meth. Fluids 0984). 9. P. N. Swartztrauber and R. Sweet, Efficient FORTRAN subprograms for the solution of elliptic partial differential equations. N C A R - T N / I A - 1 0 9 (1975).
Finite difference solution to the flow between two rotating spheres
91
10. O. R. Burggraf, Analytical and numerical studies of the structure of steady separated flows. J. Fluid Mech. 24, 113 (1966). 11. B. R. Munson and D. D. Joseph, Viscous incompressible flow between concentric rotating spheres. Part 1. Basic flow. J. Fluid Mech. 49, 289 (1971). 12. D. Greenspan, Numerical studies of steady, viscous, incompressible flow between two rotating spheres. Computers and Fluids 3, 69 (1975). 13. D. Schultz and D. Greenspan, Improved solution of steady, viscous, incompressible flow between the rotating spheres. Computers and Fluids 7, 175 (1979). 14. M. Wimmer, Experiments on a viscous fluid flow between concentric rotating spheres. J. Fluid Mech. 7g, 317 (1976). 15. J.-P. Bonnet et T. Alziary de Roquefort, Ecoulement entre deux sphrre concentriques en rotation. J. Meca. 15, 373 (1976). 16. E. Krause and F. Bartels, in Approximation Methods for Navier-Stokes Problems, R. Rautmann, Editor. Springer, Berlin (1980). 17. G. Schrauf, Branching of Navier-Stokes equations in a spherical gap. Proc. 8th Int. Conf. Num. Meth. Fluid Dyn. (Edited by E. Krause). Lecture Notes in Physics 170, 474 (1982). 18. B. R. Munson and R. W. Douglass, Viscous flow in oscillatory spherical annuli, Phys. Fluids 22, 205 (1979). 19. Z. Zapryanov and C. Christov, Numerical study of the viscous flow in oscillatory spherical annuli. Comp. Meth. Appl. Mech. Engng 22, 49 (1980). APPENDIX The evaluation of A and c directly from the defining expressions (6.6) and (6.7) presents two disadvantages: (i) the calculation of the integral for each matrix coefficient or vector component requires O(MN) arithmetic operations; (ii) the storage of the functions r&, ~ = 1. . . . . L, which can be determined once and for all at the beginning of the calculation, requires L M N words of computer memory. It is possible to circumvent both inconveniences by a method proposed by Glowinski and Pironneau [5] at the expense of doubling the number of elliptic equations to be solved. Their method is here adapted to the case of axisymmetric problems in spherical coordinates. Let us consider two vector fields A and B with the following spherical components
A = (0, O, A¢~= A(r, 0)), a = (0, O, B~ = B(r, 0)).
For the special case of these rotationally symmetric vector fields the vector analogue of Green's theorem in three dimensions assumes the simple form
f d v w ( A E 2 B - BE2A)=
fdswn • (AVB -
BVA),
(Ai)
where, as previously, w = w(r, 0) = r sin 0 and n is the normal unit vector. Instead of considering the functions ~, defined through eqns (5.4), let us introduce the functions W~ defined as follows W~ = arbitrary in v,
r=rl2
W~=6~,,
and
(A2a) r=l.
(A2b)
For simplicity, in this Appendix we do not indicate explicitly that all functions are prescribed to be zero at 0 = 0 and 0 = n/2; this will be assumed implicitly. Let us now consider the following elliptic problems ( - E 2 + e)(~ = 0,
(~=W~,
r=rl2
and
(A3a) r=l;
- - E 2 ~ = (~, ~b~=0,
r=q2
(A3b) (A4a)
and
r=l.
(A4b)
With the aid of Green's theorem (A1), it is a simple matter to show that
(A5)
Similarly, considering the two elliptic problems ( - E 2 + e)~"= h, (=arbitrary,
r=rl2
(A6a)
and
r=l;
(A6b) (A7a)
~=0,
r=rl2
and
r=l;
(A7b)
92
S . C . R . DENNIS and L. QUARTAPELLE
it is easily obtained that c~ --- - C , .o)
(AS)
If the arbitrariness of the functions W~ at the internal points is exploited by taking W~ = 0 at all the internal points, the matrix coefficients A~,, and the vector components c~ can be evaluated through eqns (A5) and (A8) in 0(1) arithmetic operations. Furthermore, due to the boundary conditions (A2b), the functions W, can be generated explicitly whenever needed and there is no storage requirement.