A singular perturbation solution of viscous incompressible fluid flow between two eccentric rotating cylinders

The flow inside two concentric cylinders is one dimensional and an exact solution for quantities is easily found. However, when the cylinders axes are displaced by a small distance, two dimensional effects become obvious. In this research, the equations governing an incompressible viscous flow between two rotating cylinders are considered in polar coordinates that can be simplified by introducing vorticity and stream functions. By taking the curl of the vector form of the momentum equation, the pressure term is omitted. Because of the boundary conditions being in terms of perturbation parameter, a modified bi-polar coordinate system is introduced. This transforms the two eccentric cylinders into two concentric ones. By expanding the quantities in terms of the perturbation parameter up to second-order accuracy and by substitution into the vorticity-stream function, sets of differential equations are obtained to be solved for these functions. At the end, the closed-form velocity components are determined.

Many researchers tried to investigate the classical problem of lubrication theory [8,9]. They obtained an exact solution of the problem of two-dimensional creeping flow between two eccentrically arranged cylinders. These works describe the idealized flow in the gap of journal bearings which is filled by lubricating fluid. Chernyavsky [1] studied the problem of flow between circular cylinders with a free internal cylinder, in which, the dowel is under the action of gravity, the force from the fluid filling the space between the dowel and the bearing, and the action of the moment of forces with respect to the center of the dowel. Wannier [17] proved that the Reynolds equation may be obtained by expanding quantities in Stokes equation in terms of powers of very small thickness of film. Sommerfeld [15] solved the approximate equation of Reynolds in 1904, and Kamal [10] presented a brief of his solution. Wannier also represented an exact solution of Stokes equations for flow between eccentric cylinders which is in case of small Reynolds numbers. Wannier's solution does not include inertia effects. Kamal considered inertia effects, and expanded flow quantities in terms of inverse powers of kinematic viscosity. Kulinski and Ostrach [11] obtained velocity components for small eccentricity and moderate Reynolds numbers using perturbation method. Diprima and Stuart [4] worked on this and also got pressure distribution and torque. All these works are valid for small Reynolds numbers. Frene and Godet [5] studied the criteria for flow transition in a journal bearing. Wood [19] considered solutions for all ranges of Reynolds number, and expanded flow quantities in terms of eccentricity in polar coordinates. He claimed that if one of the cylinders is stationary, this expansion seems not to be valid for large Reynolds numbers, because inside the boundary layer attached to the stationary cylinders, flow quantities would not be analytic functions of eccentricity. In other words, the expansion is not valid inside the inner region near stationary cylinder, and secular terms appear in the regular expansion. Wood's solution is in Bessel functions and imaginary orders and due to its complexity it is not usually used to study flow stability. Meanwhile, it only includes first-order approximations. Diprima and Stuart proved that the Taylor number in which flow changes to vortex, is proportional to square eccentricity; hence, the second-order approximations are essential. Petrov [13] investigated the motion of an incompressible viscous fluid in a thin layer between two circular cylinders. He solved the equations of hydrodynamic lubrication for the unsteady plane parallel flow. Dai et al. [3] investigated the effect of eccentricity on Couette-Taylor transition for flow between infinite rotating cylinders. Their analysis method is Fourier expansion of the conservation equations in the axial direction, followed by projection onto a polynomial subspace.
In this research, a perturbation technique is used to solve the non-linear Navier-Stokes equations in a twodimensional region. Second-order approximations together with first-order ones are represented in terms of simple functions. To do this, the boundary layer theory is applied, in which Reynolds number will be considered a large parameter [14]. Up to now, expanding the flow quantities in terms of inverse Reynolds number, or a function of it, has not been considered. If so, we deal with a singular perturbation problem.

Problem formulation
The governing equations in cylindrical coordinates for steady-state conditions and constant density and viscosity are given by: The boundary conditions, according to the geometry given in Fig. 1, are: where, e is the eccentricity, a and b are the inner and outer cylinders' radii, respectively, and r * 0 is the distance of outer cylinder to the origin.
These cylinders are rotating at the same direction with constant angular velocities Ω 1 and Ω 2 and the origin of the coordinates system coincides with the center of the inner cylinder. If the following dimensionless variables are introduced: then, Eqs. (1)-(3) will be changed to the following dimensionless form: 1 r where, the Reynolds number is defined as: The dimensionless form of the boundary conditions (4) is: where, and ε = e a (17) Now, if the velocity components are defined by stream function [18]: the continuity equation (11) will automatically be satisfied. Then by taking derivative of Eqs. (12) and (13) with respect to θ and r , respectively, and combining them properly, the pressure term will be omitted and the following equation will be obtained in terms of vorticity function [6]: where, Boundary conditions (15) are now written in terms of ψ as follows: Equations (19) and (20) with boundary conditions (21) are the mathematical model of motion between two eccentric rotating cylinders.

Equations of motion in modified bi-polar coordinates system
As the coordinate describing outer cylinder surface is a function of ε, finding approximate solution by expanding the quantities in terms of ε is not possible. Hence, the modified bi-polar coordinates system is introduced to overcome this problem [19]. Using conformal mapping, we have: where, z = r e iθ and ζ = e iϕ and Consequently, the inner and outer cylinders will be changed to the coordinate lines = 1 and = β, respectively, in which [2]: Figure 2 illustrates the superposition of modified bi-polar coordinates system on the cylindrical one. These new coordinates' axes are exactly the same as the two-dimensional bi-polar ones; however, variables and ϕ have been modified in such a way that when the cylinders are concentric, they will be changed to r and θ . The stream function is related to the velocity components in this new coordinates system as follows [16]: The equations of stream and vorticity functions are: where, Then, boundary conditions (21) become: For these boundary conditions, we have q 1 = aΩ 1 and q 2 = bΩ 2 . Equations (27) and (28) with boundary conditions (30) are the governing equations of flow between two rotating cylinders in modified bi-polar coordinates system.

Linearization of the motion equations
As γ is a measure of eccentricity ε, it can be used to linearize the motion equations by asymptotic expansions of the stream and vorticity functions [12]: Also, J may be expanded as follows: where, Substituting Eqs. (31)-(33) into (27) and (28), and equating powers of γ to zero, gives: When the parameter γ equals zero, the cylinders are concentric, and its flow is one-dimensional. Hence, flow quantities of zero order, that is, ψ 0 , and ω 0 , are independent of ϕ. In this case, vorticity is constant everywhere in the flow field. By considering these facts, Eqs. (35)-(37) are simplified as: Expanding boundary conditions (30) in terms of powers of γ gives:

Solution of zeroth-order equations
Solution of Eq. (38) with boundary conditions (44) is: where,

Solution of first-order equations
Substituting ψ 0 into Eqs. (40) and (41) gives: As can be seen in the above equations and boundary conditions (45), the flow quantities must be functions of ϕ, hence, ψ 1 and ω 1 must take the following form: Substituting these functions into Eqs. (49) and (50), and rearranging them, gives: Boundary conditions (45) are then changed for f c and f s as: Equations (53) and (54) are now solved for large Reynolds numbers. So we deal with a boundary layer problem, or singular perturbation in which a small parameter, named perturbation parameter, is multiplied by the greatest derivative of the equation. In this case, first, using an ordinary series expansion, an outer solution is obtained that is only valid outside the boundary layer. Then, by a proper variable change, an inner solution will be obtained that is valid inside the boundary layer. Finally, these inner and outer solutions will be matched to get a valid composite solution throughout the flow field.

Outer solution
If the following expansions: are replaced into Eqs. (53) and (54), then by equating the powers of λ = Re −1/2 to zero, we get the following equations: Their solutions are: where, the constants c and d will be determined by matching outer solution with inner one. So the outer expansions for f c and f s are as follows: Constants c and d can be determined by matching these with inner solutions.

Inner solution near = 1
To obtain a solution near = 1, the following variable change is introduced: Note that, here the thickness of boundary layer is of order of λ. To find inner solution, we assume the following expansions: Substituting these expressions into Eqs. (53) and (54), and equating the powers of λ to zero, gives: The boundary conditions (55), relating to = 1 or ξ = 0, become: Solving Eqs. (65) and (66) with boundary conditions (67) leads to the following results for the inner solution: Constants c 0 , d 0 and c 1 , d 1 can be determined by matching with outer solution.

Inner solution near = β
To obtain a solution near = β, the following variable change is introduced: If the following expansions: are substituted into Eqs. (53) and (54), and equating the powers of λ to zero, the following equations will be obtained: in which, As can be seen, for q 2 = 0, δ equals zero. In this case, solutions of each of the Eqs. (72) or (73) will not be decreased exponentially. Hence, these solutions cannot be matched with outer solution. Because of this, when one of the cylinders is stationary, expanding stream function will not be analytic. So this study will be limited to q 2 > 0. The boundary conditions (55), relating to = β or η = 0, become: Constants C 0 , D 0 and C 1 , D 1 can be determined by matching with outer solution.

Composite solution
The above three solutions are now combined to get the composite solution which is given as follows: in which, ( f o ) i is inner expansion of outer solution, and the same for ( f o ) I . According to the matching principle of these solutions, inner expansion of outer solution must be equal to outer expansion of inner solution. So we have: The unknown constants resulted in the outer and inner solutions can be determined by the above equalities. The composite solutions are: Finally, approximate solution of ψ 1 is:

Solution of second-order equations
The second-order equations of problem are: in which, subject to the following boundary conditions: It is clear that ψ 2 and ω 2 vary as follows: Substituting these relations into Eqs. (84) and (85), we have:

Outer expansion
If the following single term approximations for X and Y X = X 0 + · · · and Y = Y 0 + · · · (98) are substituted into Eqs. (96) and (97), and collecting equal terms, the zeroth-order equations results: General solutions of these equations are: In which, constants a 0 , b 0 , c 0 , d 0 may be determined by matching with inner solutions.

Inner expansion near = 1
By introducing the variable change: and then, substituting the following single term expansions: into Eqs. (96) and (97), and equating similar powers of λ to zero, we have: Solutions of these equations subject to the related boundary conditions are: in which, constants A 2 and A 3 may be determined by matching with outer solution.

Inner expansion near = β
By introducing the variable change: and then, substituting the following single term expansions: into Eqs. (96) and (97), and equating similar powers of λ to zero, we have: in which, δ = Aβ 2 + B. Solutions of these equations subject to the related boundary conditions are: in which, constants B 2 and B 3 may be determined by matching with outer solution.

Composite solution
The composite solutions for X and Y are: Using the following matching conditions: the unknown constants in Eqs. (100), (104) and (108) are obtained, and so the composite solutions become: Therefore ψ 2 becomes:

Velocity components
Stream function ψ can be determined by Eqs. (47), (83) and (112) as follows: Now, by considering velocity components in terms of ψ in Eq. (18), and also having conformal mapping in Eq. (22), we can compute the velocity field using the chain rule as follows: In the case of concentric rotating cylinders, the only non-zero velocity component is v, and the simplified non-dimensional Navier-Stokes equations give the following dimensionless rotational velocity component: Our analysis results coincide exactly with this one-dimensional state as the eccentricity vanishes.

Results and conclusions
In this study, a closed-form mathematical solution of flow between two eccentric rotating cylinders is presented which is of great importance in calculating velocity profiles inside the flow field. The velocity distributions for ε = 0.1, b/a = 1.2, q 2 /q 1 = 0.5 and three different Reynolds number at four various angles (θ = 0, π/2, π, 3π/2) are represented in Fig. 3. As the gap size is decreased from θ = 0 to θ = π, it is observed that the velocity values are increased to compensate the same mass flow rate. Figure 4 also illustrates the velocity contours of the whole domain. As mentioned in [18], small-gap stability theory, for 0 ≤ Ω 2 /Ω 1 ≤ 1, predicts that instability occurs at Ta crit ≈ 1708. According to our notations, the Taylor number is defined as Ta = a(b−a) 3 (Ω 2 1 −Ω 2 2 ) ν 2 ; it may be re-written as Ta = Re 2 ( b a − 1) 3 [1 − Ω 1 and the Reynolds number is defined by (14). Hence, the critical Reynolds number for the case of b/a = 1.2 and q 2 /q 1 = 0.5 is Re crit ≈ 508. The results presented here are therefore in the range of Taylor stability criterion.
In this study, we deal with a singular perturbation problem arising from the special interesting behavior of the flow field near the solid boundaries. Another advantage of this solution method is to identify the physical behavior of flow, especially near the solid cylindrical walls where the boundary layers are formed. In fact, a multi-layer region is developed in the flow field at high Reynolds numbers; an outer layer as well as some inner  layers may be distinguished inside the domain. This solution is valid for small eccentricity at high Reynolds numbers, when the two cylinders are rotating at the same direction. So this approximate analytic solution can be a good reference to compare some other solutions.
Open Access This article is distributed under the terms of the Creative Commons Attribution License which permits any use, distribution, and reproduction in any medium, provided the original author(s) and the source are credited.