A Finite Element Method for Incompressible Fluid Flow in a Noninertial Frame of Reference and a Numerical Study of Circular Couette Flow with Varying Angular Speed

Abstract: In this work, we develop a numerical strategy for analyzing the flows of an incompressible fluid in the gap between an arbitrarily shaped inner boundary that rotates inside a circular outer boundary. Such flows occur very commonly in turbomachinery applications. The numerical strategy is based on a noninertial frame of reference that is fixed to the rotating inner boundary so that Coriolis and angular acceleration effects have to be accounted for in its development. Since this strategy is based on a fixed mesh, it is much more economical and accurate than a general arbitrary Eulerian–Lagrangian strategy, which would typically require remeshing. In addition, we also conduct a numerical study for circular Couette flow with varying angular speed of the inner cylinder in an inertial frame of reference; such a study may prove useful in validating a theoretical stability analysis which currently seems to have been carried out only for the case of constant angular speed.


Introduction
This work develops a velocity-pressure integrated finite element formulation in a non-inertial frame for the general case of three-dimensional incompressible fluid flow without any simplifying assumptions such as potential flow.Most of the existing work has focused on inertial frame formulations.Kim [1] presented a velocity-pressure integrated mixed interpolation method for the analysis of fluid flows.Kim and Decker [2] compared solutions obtained by the continuous-pressure formulation and penalty formulation, and concluded that both have comparable accuracy.However, they showed that the penalty formulation is computationally efficient, as the pressure is eliminated at the element level and treated implicitly.
Jog and Rakesh Kumar [3] compared these formulations and concluded that the penalty formulation can lead to inaccurate solutions on a class of transient problems (although on some other problems it can yield accurate solutions as well).They also showed that the velocity-pressure formulation is more accurate and robust, and hence is to be preferred in spite of its slightly higher cost.The design of marine propellers and turbines involves problems with moving/rotating boundaries.For fluid problems with obstacles or moving boundaries, the fluid mesh may distort and has to be updated accordingly.Such problems can be solved by the arbitrary Lagrangian-Eulerian (ALE) formulation with automatic and continuous rezoning of the fluid mesh [4][5][6].Besides the problem of remeshing, projection errors also occur in projecting the solution from the old mesh to the new one.
Another approach to dealing with such moving boundary problems is to use a blade-fixed rotating frame of reference.In this formulation, there is no need to update the mesh, but body forces such as centrifugal and Coriolis forces have to be included in the formulation.Young [7] used this approach to solve the fluid-structure interaction problem for marine propellers.The finite element and boundary element methods (FEM and BEM) were used to model the solid and fluid parts, respectively, with nodal displacements and pressure as the variables.The hydrodynamic pressure due to the rigid rotation was obtained via the boundary element method (BEM), and this was applied as an external traction on the solid structure.In [8,9], a time-dependent hydroelastic analysis was conducted where the pressure obtained by the BEM was used to find the hydrodynamic mass and damping matrices, which were superimposed onto the structural mass and dynamic matrices.Potential flow theory (where the fluid is assumed to be incompressible and inviscid and the flow assumed to be irrotational) was used to find the pressure field.In [10], non-uniform B-splines are used to accurately represent circular geometries.The domain was subdivided into rotating and stationary domains and velocity compatibility was enforced at the interface.In [11], large eddy simulations were performed on the flow in a baffled stirred tank using a lattice-Boltzmann scheme for discretizing the Navier-Stokes equations.
In this work, we attempt to solve the three-dimensional Navier-Stokes equations (without any simplifying assumptions such as potential flow) in a blade-fixed frame of reference.In this approach, non-inertial body forces such as centrifugal and Coriolis forces have to be incorporated into the formulation.As a first approximation, we assume the structure to be rigid.In this approach, since we use a fixed mesh, the outer boundary is constrained to be circular.Since the outer boundary is circular in most turbomachinery applications, this approach is applicable to a large class of problems.This strategy thus offers a computationally efficient method to solve fluid flows with an arbitrarily shaped inner rotating boundary.The proposed method can be generalized by considering the flexibility of the blades, in which case, it has to be solved as a fluid-structure interaction problem.
Another problem that we address in this work is as follows.In circular Couette flow, the flow becomes unstable after a certain threshold Reynolds number.The stability analysis of this class of flows was carried out by Taylor [12].He found the critical Reynolds number Re cr for the onset of instabilities, both theoretically and experimentally, for different radii ratio.Meyer and Keller [13] carried out a numerical stability analysis of this class of flows.They used a Fourier expansion and a centered finite difference technique to find Re cr for different radii ratio.Moser, Moin, and Leonard [14] used a spectral method for solving the Navier-Stokes equations.They presented values of Re cr for a narrow as well as a wide gap between cylinders.The results obtained using our formulation are in agreement with those of [12][13][14].In all the above theoretical or numerical studies, the analysis has been restricted to a constant angular speed of rotation of the inner cylinder.In this work, we conduct a numerical study of circular Couette flow for varying angular speeds of rotation of the inner cylinder with the outer cylinder fixed; such a study is not only of interest in its own right, but could also prove useful in validating a theoretical stability analysis (which typically involves many simplifying assumptions).

Strong Form of the Governing Equations
Let Ω denote the domain whose boundary Γ is composed of two regions, Γ u and Γ t .We are interested in finding an approximate numerical solution to the following initial-boundary value problem.
Find the velocity u, stresses τ, rate of deformation D, and tractions t such that where u represents the velocity, p is the pressure, µ is the dynamic viscosity coefficient, n is the outward normal to Γ, b is the prescribed body force on Ω, t and ū are the prescribed tractions and velocities on Γ t and Γ u , respectively, and u 0 is the initial velocity prescribed on the domain.In the velocity-pressure formulation described below, we need that either the tractions or the velocities be prescribed on the boundary.

Noninertial Frame Velocity-Pressure Integrated Formulation
The noninertial frame formulation presented here is a modified version of the inertial frame formulation presented in detail in [3].In the velocity-pressure integrated formulation, Equation ( 1) is implemented in a weak sense as where the subscript δ denotes the variation of the quantity in question.Similarly, we also enforce Equations ( 2) and ( 6) in a weak sense as In the rotating (noninertial) frame of reference, the additional body forces that have to be considered are [15] Centrifugal force = ω × (ω × x), where ω denotes the angular velocity vector of the rotating (noninertial) frame of reference, and x denotes the position vector.Since it can be written as the gradient of a scalar field, the centrifugal force can be absorbed into the pressure field, and hence we need to consider only the latter two forces.The variable p thus denotes the modified pressure that incorporates the centrifugal force effects.By using the identity ∇ • (τ T u δ ) = u δ • (∇ • τ) + ∇u δ : τ, the constitutive relation given by Equation ( 3) and the divergence theorem, the above equation in a noninertial frame simplifies to where b now denotes the body force in the inertial frame of reference.For the problem at hand, the angular velocity vector is of the form ω = [0, 0, ω] T , W represents the skew-symmetric tensor of which ω is an axial vector, D is the rate of deformation tensor, and C is the fourth-order constitutive tensor which can be written in engineering form [15] as In the continuous-pressure formulation, we use different interpolation functions for the pressure and velocity fields.Let the velocity field and pressure field, their variation, and increments (denoted by subscript ∆) be interpolated as The shape functions N are the standard higher-order Lagrange shape functions used in an isoparametric formulation, and N p are the shape functions used for pressure field.Using the above interpolations, we have The semi discrete form is given by where By using the generalized trapezoidal rule, we get the fully discretized form of the equations as where Based on stability considerations, one typically uses α ∈ [1/2, 1], where α = 1/2 and 1 correspond to the Crank-Nicolson and backward Euler methods, respectively.At each time step, we solve for the nodal pressure and velocity degrees of freedom.
Once the velocity u and vorticity ∇ × u fields in the noninertial frame of reference are found using the above formulation, the angular velocity and nonzero component of the vorticity vector (denoted by γ) in the inertial frame of reference are recovered using where r represents the radial distance from the center.

Numerical Example
We use 27-node brick elements for meshing the domain in all of the following examples.In all examples, mesh and time step refinement were carried out to ensure convergence of the results with respect to mesh and time step refinement.The noninertial frame solution strategy presented in the previous section is first validated by comparing it against the analytical laminar flow solution presented in [16] for the case where the inner cylinder rotates with a sinusoidally varying speed.Since the analytical solution is presented in an inertial frame of reference, the numerical solution in the noninertial frame of reference is transformed to an inertial frame of reference using Equation (13a) before carrying out the comparison.A numerical solution is also directly found in the inertial frame by imposing a variable velocity on the inner cylinder.We assume the angular speed of the inner cylinder to be ω(t) = ω 0 sin( f t) and that the outer cylinder is fixed.The top and bottom surfaces of the cylinder are traction free.The radii of the inner and outer cylinders are 0.2 and 0.5, respectively.The fluid properties used are µ = 10 −3 and ρ = 1.The mesh used for finding the solution is n r × n θ × n z = 15 × 18 × 1, and values of the various parameters used are α = 1/2, f = 1.0, ω 0 = 0.1, and t = 10.An almost exact match is found between the numerical and analytical solutions as shown in Figure 1.Having validated the noninertial frame formulation, we now present solutions for a problem where no analytical or numerical solutions are available.The problem is to find the flow in the annulus between two concentric cylinders with a rigid plate attached to the inner rotating cylinder (see Figure 2).The angular speed of rotation of the inner cylinder is again assumed to be ω(t) = ω 0 sin( f t).We study the effect of various parameters controlling the flow characteristics-namely, the ratio of plate length along the radial direction to the radial gap between the cylinders, the frequency of oscillation of the inner cylinder, etc.We assume the inner and outer radii, r i and r 0 , to be 0.2 and 0.5, respectively.The shape of plate used is wedge-like with central angle 5 • .A 15 × 36 × 1 mesh is used to model the domain.The properties used are µ = 10 −3 and ρ = 1.In the subsequent discussion, the velocity flow field is plotted in the mid-plane perpendicular to the axis of the cylinder.It will be seen that vortices form near to the tip of the plate, while the velocity field is almost uniform sufficiently far away from the plate.
The velocity fields for a width/gap ratio of 0.4 between the plate and outer cylinder are obtained for ( f , ω 0 ) = (0.5, 0.05), (1, 0.1), and (2, 0.2). Figure 3 shows the velocity fields for the different cases at different times.Similarly, a typical vortex formation at different frequencies at a fixed instant of time is shown in Figure 4. is seen that a vortex forms on the trailing side of the plate, and is clockwise if the plate is moving anti-clockwise at that instant, and vice versa.

Instability in Couette Flow
Most of the literature on instabilities in the case of circular Couette flow (see Figure 5) has focused on the case where the inner cylinder rotates with a constant angular speed.We now present a numerical study to find the critical values of the parameters for the onset of instabilities in case the inner cylinder has a sinusoidally varying speed.It is hoped that these results will prove useful for validating an analytical instability analysis which so far seems to have been confined to the case of constant angular speed.
Taylor [12] was probably the first to conduct such a stability analysis.He found the critical Reynolds number at which the flow becomes unstable for different radii ratios, both analytically and experimentally.At low Reynolds number the flow is azimuthal.As we keep increasing the angular speed of the inner cylinder, at a certain critical Reynolds number Re cr the flow becomes unstable, and the radial and axial velocity components become nonzero.Taylor studied the effect of the radii ratio on this critical Reynolds number.

Constant Angular Speed
Before presenting our results for the varying angular speed case, we first numerically verify that the flow does indeed become unstable at the predicted Re cr (see Table 1).Such a numerical verification is important since Taylor's work is based on linearized stability analysis, which ignores both the nonlinear and transient terms.Since the radial velocity component starts growing (from zero) with the onset of the instability, one can find Re cr by monitoring this radial and azimuthal velocity component as shown in Figures 6 and 7.The contour plots for the radial and azimuthal velocity components are similar to those presented in [13], and are shown in Figure 8.

Sinusoidally Varying Angular Speed
We now find the critical Reynolds number Re cr for the case where the angular speed of rotation varies with time as ω(t) = ω 0 sin( f t).The radius of the inner cylinder used is r i = 0.0475 and the length of the cylinder is chosen to be five times the radial gap.The radius of the outer cylinder r o is varied to study the dependence of the gap on the flow.The fluid properties used are µ = 2.5 × 10 −3 and ρ = 100.An 8 × 8 × 5 mesh is used to model the domain.For variable angular speed, the Reynolds number is defined based on the amplitude ω 0 as Re = ρr i ω 0 (r o − r i )/µ.The effect of various parameters, such as radii ratio η and a nondimensional number based on the frequency f , on the critical Reynolds number are studied.
As mentioned, if there is no instability in the fluid flow, it will remain azimuthal with zero radial velocity.However, beyond a certain critical Reynolds number, we clearly observe a non-zero radial velocity component, which is an indication of a centrifugal instability in the fluid.In Figures 9 and 10, the radial and azimuthal velocities at r = (r i + r o )/2 for varying angular speed are plotted against the number of oscillations.It is seen that for Re = 100, the radial velocity u r is zero, while the azimuthal velocity u θ follows a sinusoidal pattern.Around Re = 140, the u r component starts growing, and u θ starts deviating after five cycles, while at Re = 180, this deviation occurs after two cycles.This implies that an instability is induced in the fluid at around Re = 140.At larger Re, the instability occurs even earlier.In a similar manner, Re cr can be obtained for different frequencies keeping the geometry fixed, and for different radii ratio keeping the frequency constant.For a geometry with radii ratio η = 0.7 (narrow gap) and η = 0.4 (wide gap), we have studied the dependence of the critical Reynolds number (Re cr ) on the dimensionless frequency parameter,

Figure 1 .
Figure 1.Comparison of the analytical and numerical solutions in the inertial and noninertial frames of reference.

Figure 2 .
Figure 2. Flow due to oscillation of inner cylinder with a rigid plate attached to it.

Figure 6 .
Figure 6.Variation of the radial velocity with radius for different Re for the constant angular speed case.

Figure 7 .
Figure 7. Variation of the azimuthal velocity with radius for different Re for the constant angular speed case.

Figure 8 .
Figure 8. Contours of the radial and azimuthal velocity components for the constant angular speed case.

Figure 9 .
Figure 9. Variation of radial velocity with time for different Re for the varying angular speed case.

Figure 10 .
Figure 10.Variation of azimuthal velocity with time for different Re for the varying angular speed case.

Table 1 .
Critical Reynolds number for the onset of instabilities for the case of constant angular speed.