NUMERICAL STUDIES OF VISCOUS INCOMPRESSIBLE FLOW BETWEEN TWO ROTATING CONCENTRIC SPHERES

The problem of determining the induced steady axially symmetric motion of an incompressible viscous ﬂuid conﬁned between two concentric spheres, with the outer sphere rotating with constant angular velocity and the inner sphere ﬁxed, is numerically investigated for large Reynolds number. The governing Navier-Stokes equations expressed in terms of a stream function-vorticity formulation are reduced to a set of nonlinear ordinary di ﬀ erential equations in the radial variable, one of second order and the other of fourth order, by expanding the ﬂow variables as an inﬁnite series of orthogonal Gegenbauer functions. The numerical investigation is based on a ﬁnite-di ﬀ erence technique which does not involve iterations and which is valid for arbitrary large Reynolds number. Present calculations are performed for Reynolds numbers as large as 5000. The resulting ﬂow patterns are displayed in the form of level curves. The results show a stable conﬁguration consistent with experimental results with no evidence of any disjoint closed curves.


Introduction
Confined rotating flows are usually solved using numerical techniques that involve iteration. This often reduces the region of applicability of the numerical method as a result of the instabilities arising from the iterative process on which the numerical method is based. In this paper, we apply a noniterative finite-difference technique to the study of the induced steady axially symmetric motion of an incompressible viscous fluid confined between two concentric spheres with the outer sphere rotating with constant angular velocity and the inner sphere fixed. The application of a new numerical technique to the rotating spheres problem is appropriate for the following reasons. First, the motion of viscous fluids in rotating containers is among the most fundamental problems in classical fluid dynamics and is of basic importance in finding a general hydrodynamic stability theory for rotating fluids in an enclosed cavity volume. Understanding the fluid dynamics of the motion is relevant in geophysical studies such as the study of atmospheric and oceanic circulations and also engineering design, particularly the study of gyroscopes and centrifuges. Second, the need to apply a new numerical technique to the rotating sphere problem is necessitated by the fact that there is a need to check the results in Greenspan [2], which differ substantially from those of Pearson [11], for Re = 1000. While Greenspan's result suggests the existence of two disjoint closed curves at the central portion of the flow region, Pearson's result does not predict the existence of such a phenomenon. Furthermore, Greenspan [2,Figure 14] shows the appearance of several disjoint closed curves for Re = 3000 while Schultz and Greenspan [13,Figure 6] show that there are only two sets of closed curves in the central region. Finally, researchers are always interested in numerical methods that are accurate and stable for a wide range of flow parameters. Many papers have been published on the numerical solution of steady, viscous, incompressible flow between two rotating concentric spheres. However, these papers have used iterative techniques and have produced results only for Re ≤ 3000. We present a noniterative method which is valid for both small and large Reynolds number and obtain results as large as Re = 5000.

Governing equations of motion
We consider the motion of an isothermal incompressible viscous fluid contained in a spherical annulus between two rotating concentric spheres. The spheres are assumed to be rigid and rotate freely about a vertical axis through their centre, which is taken to be the origin of the coordinate system. The inner sphere and outer sphere are constrained to rotate with prescribed angular velocities ω 1 and ω 2 , respectively. We start with the full incompressible Navier-Stokes equation with no body force [14], together with the incompressibility condition, ∇ · q = 0, (2.2) and the no-slip (rigid) boundary conditions on the spherical surfaces, where q = (u,v,w) is the velocity field, p is the kinematic pressure, and Re is the Reynolds number defined by ω 0 r 2 2 /v. In the two extreme cases, ω 0 = ω 1 if the outer sphere is stationary and ω 0 = ω 2 if the inner sphere is stationary. When both spheres rotate, the character of the flow is generally dominated by the motion of either the inner or outer sphere. The angular velocity of the dominant sphere is then taken as the characteristic angular velocity. The specific choice of ω 0 is indicated for each example considered. The governing equations have been nondimensionalized with the characteristic velocity scale ω 0 r 2 and the characteristic length scale r 2 . If the spherical polar coordinates (r,θ,ϕ) are chosen with the origin at the centre of the sphere and θ = 0, the axis of rotation, the motion is independent of the azimuthal angle ϕ and the Navier-Stokes equations can then be written in terms of a stream function in the meridian plane, Ψ, an angular velocity function Ω, as well as a vorticity function Φ. Thus, making the substitutions, the Navier-Stokes equations (2.1) and (2.2) are easily found to be where In view of the boundary conditions, we may assume the expansions where I n (µ) are the Gegenbauer functions of argument µ = cos θ. These functions form a complete orthogonal set in the range −1 ≤ µ ≤ 1. We have substituted r = e −x . Substituting (2.9)-(2.11) in the equations (2.5)-(2.8) and making use of the orthogonality properties of the functions I n (µ) give rise to infinite sets of ordinary differential equations. Substitution into (2.5) yields where R n (x) are nonlinear combinations of the various functions defined by (2.18) The quantities d 1 (n,m,k) and d 2 (n,m,k) which appear in (2.13), (2.17), and (2.18) are integrals involving products of Gegenbauer functions and derived functions and are given by The boundary conditions corresponding to (2.12) are where c i = −log e b i , b i = r i /r 2 are the surfaces of the spheres, h i = ω i /ω 0 , and δ n,1 is the Kronecker delta function. The boundary conditions for (2.14) are The functions Y n (x) are not known at the surface of the sphere and can be calculated as part of the solution. One way of doing this is to multiply (2.14) by e −(2n+1)x and e 2nx in turn, then integrate both sides of the equation with respect to x from c 1 to c 2 . After some integration by parts of the terms on the left-hand side and the use of the conditions The use of (2.22) poses additional computational difficulties as their values need to be determined along with the other flow variables during the computational process. Conditions (2.22) act as an initial guess for the boundary conditions for Y n (x) and are to be improved upon iteratively during the computational process. In principle, any quadrature formulae may be used but in practice; care needs to be taken about the choice of quadrature used as this may lead to inaccuracies and subsequent instability of the overall computational process. Since we are interested in developing a computational procedure that does not involve any iteration, at least in the conventional sense, it would be computationally more advantageous if the fourth-order equation is considered instead. Thus we substitute ψ n (x) = e −2x S n (x), so that (2.14) becomes and (2.15) becomes In the practical numerical integrations, the computations must be limited to the determination of a finite number of terms ns in each of the series (2.9), (2.10), and (2.11) by truncation of the series. We will determine later an appropriate choice for ns, for which terms with subscript greater than ns will be set to zero.

Computational method
Before we embark on the discretization of the governing equations in the computational domain between the two concentric spheres by a suitable finite-difference method, we observe from the previous numerical experiments that in finding numerical solutions to viscous incompressible flow problems, one discovers that for high Reynolds numbers, that is, Re > 400, one can use the solutions for lower Reynolds number flow as a starting point in an iterative process to get solutions for higher values of Re. Convergence of the overall iterative process may not be achieved using other starting values. Thus the results of the computations are quite sensitive to the choice of the initial guess approximations. This suggests that one can obtain a relationship between the initial guess and the final solution. This is possible using homotopy methods. Homotopy, continuation, or embedding methods are well known and have been used in the past for finding zeros of nonlinear functions [1]. Recently, the method has been given renewed attention and it is being used for the computation of matrix eigenvalues and eigenvectors [10] as well as for obtaining analytic solutions to nonlinear partial differential equations (p.d.e) [3,4]. The basic idea is to construct a homotopy or mapping from a trivial or known solution (initial guess) to the desired solution. This is achieved by first defining a homotopy or deformation in an embedding parameter λ ∈ [0,1] of the form where G : R n → R n is a trivial or known solution and H is assumed to be smooth. Typically, one chooses a homotopy such that and attempts to trace an implicitly defined curve c(s) in the parametric set of this homotopy map from a starting point to a solution point. In order to ensure that this curve intersects the homotopy level λ = 1 at finite points, it is sufficient to impose suitable boundary conditions which essentially prevent the curve from running to infinity before intersecting the homotopy level λ = 1 or from returning to the initial level. This basic concept may be extended to a system of nonlinear p.d.e. given by

∆F(x)
The starting problem ∆ 0 G(x) = 0 should meet two requirements. First, its solutions G(x) should be easy to compute. This is obviously best met by a p.d.e whose solutions are already known. Second, its solutions must satisfy the same initial and boundary conditions as those of the original problem. The corresponding homotopy equation can easily be derived from (3.4) by replacing ∆ with ∆(x,λ) and F(x) with H(x,λ). Thus the homotopy equation is or, using (3.6), We see from (3.8) that the process of λ increasing from zero to one is just the process of H going from G(x), the initial guess approximation, to F(x), the desired solution. This process is called deformation in topology. Equation (3.8) would be referred to as the zeroth-order deformation equation. We apply this concept of homotopy to the solutions of (2.12) and (2.24). For (2.12), we start with the linear problem subject to the boundary conditions where Q(x) is a nonzero arbitrary real-valued function. We then go on to solve the nonlinear equation subject to (2.20). Observe that Q(x) completely determines ∆ 0 and its selection becomes very important. We construct a family of ordinary differential equations in an embedding parameter λ such that (1 − λ)∆ 0 W n (x,λ) + λ∆W n (x,λ) = 0 (3.12) with the corresponding boundary conditions W n c 1 ,λ = 2h 1 b 2 1 δ 1n , W n c 2 ,λ = 2h 2 b 2 2 δ 1n . (3.13) When λ = 0, we see from (3.12) and (3.9) that W n (x,0) = Ω n (x). (3.14) Also, when λ = 1, we see from (3.12) and (3.11) that Thus the process of λ increasing from zero to one is just the process of W n varying from Ω n (x) to W n (x). Obviously, (3.14) and (3.15) give a relationship between the initial guess Ω n (x) and the solution W n (x). However, this kind of relationship is indirect; a direct relationship between Ω n (x) and W n (x) must be obtained. Following [3,4], we assume that the deformations W n (x,λ) are smooth enough about λ = 0 so that all of the mthorder deformation derivatives  The right-hand side term q m,n (x) is given by Thus, through this approach, one solves an infinite sequence of linear subproblems instead of the original nonlinear equations (3.11). The initial guess approximation can be obtained from (3.9) and (3.10) given the value of Q(x). Neglecting the nonlinear terms in (2.24) should give us the desired initial guess for (2.24). However, the last nonlinear term on the right-hand side of (2.24) provides the "driving force" for this equation and so cannot be neglected. One can therefore start with the problem where Z(x) is an arbitrary nonzero real-valued function and with boundary conditions φ n c 1 = φ n c 2 = dφ n dx c 1 = dφ n dx c 2 = 0, (3.27) then go on to solve the nonlinear problem subject to the boundary conditions (2.26). The choice of Z(x) is very important as it completely determines ∆ 0 . The zeroth-order deformation equation in an embedding parameter λ is then with the corresponding boundary conditions Thus, as λ increases from zero to one, Ψ n (x,λ) increases from the initial guess approximation to the desired solution S n (x). Therefore, (3.31) and (3.32) give a relationship between the initial guess φ n (x) and the solution S n (x). However, this kind of relationship is indirect. A direct relationship between φ n (x) and S n (x) can be obtained as follows. Assume that the deformations Ψ n (x,λ) are smooth enough about λ = 0 so that all of the mth-order deformation derivatives Assuming the above Taylor series is convergent at λ = 1, we have from (3.33) that The above expression gives a direct relationship between the initial guess approximation φ n (x) and the desired solution S n (x) through the unknown terms ψ m,n (x) whose governing equation can be obtained from (3.30) by differentiating the zeroth-order deformation equation (3.30) m times with respect to λ, then setting λ = 0, and finally dividing by m!.
Thus we have where the linear differential operator ∆ m,n is given as and the right-hand side term h m,n (x) is defined as e −x ψ k,l (x)d 1 (2n + 1,2i + 1,2l + 1), × d 1 2n + 1,2i + 1,2l + 1 − d 2 (2n + 1,2l + 1,2i + 1) − 2ψ k,l (x)d 1 (2n + 1,2i + 1,2l + 1) , with the corresponding boundary conditions The basic problem now is to determine a suitable choice for Q(x) as well as for Z(x), then to go on to solve the set of ordinary differential equations (3.21) and (3.39) subject to conditions (3.22) and (3.43) using an appropriate finite-difference method. In the practical numerical integration of these equations, the computations must be limited to the determination of a finite number of terms in each of the series (2.9), (2.10), and (2.11) by truncation of the series. A truncation of order ns is defined to be the process of setting to zero all the terms with subscript greater than ns in (2.9), (2.10), and (2.11) and likewise in (2.12) and (2.25) and solving the resulting set of 2ns differential equations subject to the prescribed boundary conditions. Also, in the practical evaluation of the series (3.20) and (3.38), we assume we have summed up enough terms when the term we have just added is smaller than some small ε times the magnitude of the sum thus far accumulated. If the terms do not get small fast enough, the series becomes unsuitable for numerical computations. The convergence of a series such as (3.20) and (3.38) as given by the homotopy analysis method depends upon the choice of the initial approximation and the choice of the auxiliary linear operator and the auxiliary functions Q(x), Z(x). As long as the initial approximation is not too bad, one can get convergent results by choosing a proper value for Q(x), Z(x), although a worse initial approximation often corresponds to a slowly convergent series. There are many ways to choose a best initial approximation; also, one can apply a generalized Taylor theorem to greatly enlarge convergence regions of approximation series (3.20) and (3.38) [5,6]. The numerical method uses a novel form of finite-difference representations, which involves exponential terms and is superior to standard truncated Taylor representation and which was originally developed by Allen and Southwell for second-order ordinary differential equations [12]. The finite-difference representations of the fourth-order equation (3.39) as well as the second-order equations (3.21) are derived following a procedure that involves considering the coefficients of the differential equations as constant coefficient equations and obtaining for each equation an exact form of difference equations (exact in the sense that if the equations were constant coefficient ordinary differential   (3.45). An implementation of the TDMA is given in [15].

Results and discussion
The calculations are carried out over a range of Re from 10 to 5000. Truncations of higher order were carried out for a given value of Re by using the results of the previous truncation as a starting assumption. In this way, the number of terms in the solution was built up. For small Re, convergence was quite rapid and only a few terms were required to describe the properties of the solution quite accurately. As the Reynolds number is increased, the number of terms required to give good accuracy also increases. Solutions were carried out for several values of ns for each value of Re and some judgment was made as to the maximum value of ns required. As Re increases, more terms are required for convergence. It was found that the initial guess approximation for Re = 100 can also be used as a starting value for higher values of Re. A suitable choice of Q(x) and Z(x) needs to be made to ensure convergence of the series (3.20) and (3.38), respectively. Typically, one chooses Q(x) = Z(x) = 1 for values of Re < 1000 and Z(x) = Re for higher values. A close examination of the computer output shows that, in the exception of end points, the first term used in the series representation of the stream functions, namely, S 1 (x), is positive definite and other terms, namely, S n (x), n > 1, can change sign. Also, we observe that for the cases considered, S 1 (x) |S ns (x)|. This ensures the overall convergence of our computed solutions. The radius of the inner sphere was taken as half that of the outer sphere. The essence of this is to check the results of Greenspan [2] as well as Schultz and Greenspan [13] which differ substantially from those of Pearson [11] for Re = 1000. A plot of the stream functions for Re < 1000 agrees qualitatively with those of Pearson and are consistent with the results of Munson and Joseph [7,8] and other published results. Secondary flow patterns were clearly visible in the meridional plane and a recirculation region developed near the equator. It was observed that the secondary flow    Schultz and Greenspan. For Re = 5000, however, no comparison could be made with existing numerical results since existing methods are confined to flow smaller than or equal to 3000 due to instabilities in the iterative procedures in which they are based. Contour plots of the stream function for Re = 5000 are given in Figure 4.4. Our results are consistent with the flow visualization patterns obtained from experiments by Munson and Menguturk [9] as well as those of Wimmer [16].