1.1 Problem Description
A scenario is assumed in which fluid is enclosed between two coaxial cylinders that can be considered mathematically infinite in length. In practice, however, these typically have a finite length. This mathematical approximation allows modeling the two-dimensional case without initially considering the effects of boundary conditions.
The fluid under consideration is incompressible and Newtonian, and the flow is assumed to be laminar. The cylinders each rotate with an individual angular velocity Ω1 and Ω2, respectively. The radii of the cylinders are designated R1 and R2, where the Z-axis represents the coaxial cylinder axis.
In the literature, the resulting fluid motion is known as Couette flow, a special form of flow between rotating cylinders. The investigation of this phenomenon provides insights into complex flow conditions and is of interest for various applications in fluid mechanics, engineering, and physics.
1.2 Analytical Solution
The analytical solution is based on solving the Navier-Stokes equations for incompressible fluids. It proves advantageous to consider the present geometry using cylindrical coordinates, which simplifies the calculations and exploits the symmetry of the problem.
From symmetry, it becomes evident that the axial velocity ($$u_z$$) and the radial velocity ($$u_r$$) of the fluid are each zero, $$u_z=u_r=0$$ while the circumferential velocity is considered as a function of radius $$r$$. Likewise, the pressure $$p$$ varies as a function of radius $$r$$ $$u_\varphi=u\left(r\right);\ \ p=p(r)$$These relationships facilitate the analytical calculation of the Navier-Stokes equations and thus form the basis for developing solution approaches to address the flow problem between the cylinders.
In the same way that the velocity of a solid body can be calculated by solving Newton’s second law,
$$\mathbf{F}=\frac{d\left(m\mathbf{v}\right)}{dt}=\frac{d\mathbf{I}}{dt}$$
the velocity of a fluid (liquid or gas) can be calculated by solving the Navier-Stokes equations.
$$\frac{D\left(m\mathbf{u}\right)}{Dt}=\mathbf{F}$$
The Navier-Stokes equations are analogous to Newton’s second law and state that the rate of change of momentum of a fluid equals the sum of external forces acting on the fluid. However, the Navier-Stokes equations are applied to a finite volume of fluid rather than to a solid body.
The following general notation of the equations is often found in the literature:
$$\frac{\partial}{\partial t}\left(\rho\mathbf{u}\right)+\nabla\cdot\left(\rho\mathbf{uu}\right)=-\nabla p+\nabla\cdot\mathbf{\tau}+\rho\mathbf{g}$$
The interested reader is referred to the extensive technical literature available on this subject. Without going into the mathematical transformations, the representation of the
Navier-Stokes equations in cylindrical coordinates is applied:
- Radial momentum equation: $$r:\rho\left(\frac{\partial u_r}{\partial t}+u_r\frac{\partial u_r}{\partial r}+\frac{u_\varphi}{r}\frac{\partial u_r}{\partial\varphi}-\frac{u_\varphi^2}{r}+u_z\frac{\partial u_r}{\partial z}\right)=-\frac{\partial p}{\partial r}+\mu\left(\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial u_r}{\partial r}\right)+\frac{1}{r^2}\frac{\partial^2u_r}{\partial\varphi^2}+\frac{\partial^2u_r}{\partial z^2}-\frac{u_r}{r^2}-\frac{2}{r^2}\frac{\partial u_\varphi}{\partial\varphi}\right)+\frac{1}{3}\mu\left(\frac{1}{r}\frac{\partial}{\partial r}\left(ru_r\right)+\frac{1}{r}\frac{\partial u_\varphi}{\partial\varphi}+\frac{\partial u_z}{\partial z}\right)+\rho g_r$$
- Axial momentum equation: $$\varphi:\rho\left(\frac{\partial u_\varphi}{\partial t}+u_r\frac{\partial u_\varphi}{\partial r}+\frac{u_\varphi}{r}\frac{\partial u_\varphi}{\partial\varphi}+u_z\frac{\partial u_\varphi}{\partial z}-\frac{u_ru_\varphi}{r}\right)=-\frac{1}{r}\frac{\partial p}{\partial\varphi}+\mu\left(\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial u_\varphi}{\partial r}\right)+\frac{1}{r^2}\frac{\partial^2u_\varphi}{\partial\varphi^2}+\frac{\partial^2u_\varphi}{\partial z^2}-\frac{u_\varphi}{r^2}-\frac{2}{r^2}\frac{\partial u_r}{\partial\varphi}\right)+\frac{1}{3}\mu r\left(\frac{1}{r}\frac{\partial}{\partial r}\left(ru_r\right)+\frac{1}{r}\frac{\partial u_\varphi}{\partial\varphi}+\frac{\partial u_z}{\partial z}\right)+\rho g_\varphi$$
- Azimuthal momentum equation: $$z:\rho\left(\frac{\partial u_z}{\partial t}+u_r\frac{\partial u_z}{\partial r}+\frac{u_\varphi}{r}\frac{\partial u_z}{\partial\varphi}+u_z\frac{\partial u_z}{\partial z}\right)=-\frac{\partial p}{\partial z}+\mu\left(\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial u_z}{\partial r}\right)+\frac{1}{r^2}\frac{\partial^2u_z}{\partial\varphi^2}+\frac{\partial^2u_z}{\partial z^2}\right)+\frac{1}{3}\mu\left(\frac{1}{r}\frac{\partial}{\partial r}\left(ru_r\right)+\frac{1}{r}\frac{\partial u_\varphi}{\partial\varphi}+\frac{\partial u_z}{\partial z}\right)+\rho g_z$$
From the two-dimensional consideration, exploitation of symmetry, and neglect of gravity, the equations can be massively simplified:
(1)
$$\frac{dp}{dr}=\rho\frac{{\ u}_\varphi^2}{r}$$
(2)
$$\frac{d^2u_\varphi}{dr^2}+\frac{1}{r}\cdot\frac{du_\varphi}{dr}-\frac{u_\varphi}{r^2}=0$$
Equations (1) and (2) can be solved independently of each other. The second has a solution of the type
$$u_\varphi=ar+\frac{b}{r}$$
The constants a and b can be obtained from the boundary conditions at the cylinder surfaces. The fluid adheres to these due to the so-called “no-slip” condition and thus exhibits the same velocities as the cylinders themselves:
$$u_\varphi(r=R_1)=R_1\Omega_1$$ and $$u_\varphi(r=R_2)=R_2\Omega_2$$
$$a=\frac{\Omega_2R_2^2-\Omega_1R_1^2}{R_2^2-R_1^2};\ \ b=\frac{\left(\Omega_1-\Omega_2\right)R_1^2R_2^2}{R_2^2-R_1^2}$$
Substituting the constants yields the velocity distribution over the circumference:
(3)
$$u_\varphi=\frac{\Omega_2R_2^2-\Omega_1R_1^2}{R_2^2-R_1^2}\cdot r\ +\frac{\left(\Omega_1-\Omega_2\right)R_1^2R_2^2}{R_2^2-R_1^2}\cdot\frac{1}{r}$$
To obtain the pressure distribution in the fluid, equation (3) is substituted into (1):
$$\frac{dp\left(r\right)}{dr}=\frac{\rho}{r}\left(\frac{R_2^2\Omega_2-R_1^2\Omega_1}{r\left(R_2^2-R_1^2\right)}+\frac{r-R_1^2\frac{\Omega_1+\Omega_2}{R_2^2+R_1^2}}{-R_1^2+R_2^2}\right)^2$$
The solution to this differential equation can be determined by separation of variables:
(4)
$$p\left(r\right)=\frac{R_1^4R_2^4\rho\left(\Omega_1-\Omega_2\right)^2}{2r^2\left(-R_1^2+R_2^2\right)^2}-\frac{R_1^2R_2^4\rho\left(\Omega_1-\Omega_2\right)^2}{2\left(-R_1^2+R_2^2\right)^2}\ln{\left(R_1\right)}+\frac{2R_1^2R_2^2\rho\left(\Omega_1-\Omega_2\right)\left(-\Omega_1R_1^2+\Omega_2R_2^2\right)}{\left(-R_1^2+R_2^2\right)^2}\ln{\left(r\right)}+\frac{R_1^2\rho\left(-\Omega_1R_1^2+\Omega_2R_2^2\right)^2}{2\left(-R_1^2+R_2^2\right)^2}+\frac{r^2\rho\left(-\Omega_1R_1^2+\Omega_2R_2^2\right)^2}{2\left(-R_1^2+R_2^2\right)^2}$$
Equation (4) describes the pressure distribution in the fluid.
All sought unknowns have now been determined and we can examine a specific example with numerical values.
1.3 Specific Example
Consider two counter-rotating cylinders with the following boundary conditions:
- R1 = 70 mm,
- R2 = 100 mm
- Ω1 = 180 RPM $$ \approx\ $$18.85 rad/s
- Ω2 = -Ω1
WD-40 is used as the medium between the cylinders:
- $$\rho$$ = 970 kg/m3 (density at room temperature)
- η = 0.0216019 kg/m/s (the dynamic viscosity is needed later for the numerical analysis)
1.4 Numerical Model
For the numerical solution of this problem, we use the finite volume method in a two-dimensional solution domain. The calculation is performed as steady-state, incompressible, and laminar. The geometry, boundary conditions, and material properties are adopted as described above. Only the mesh is varied to illustrate the effects on results of such a trivial example.
Three calculations with different discretizations are performed:
- fine and structured mesh,
- medium-coarse and structured mesh,
- coarse and unstructured mesh. Only the default settings from the program were used.
1.5 Analytical vs. Numerical
We turn to the analysis of the system. In the current task, we are interested in the velocity field of the fluid and the associated pressure distribution. Additional questions such as temperature distribution or critical rotation speeds can be considered as supplementary analyses. Depending on the analysis, these are significantly more complex and must be treated specifically.
Before proceeding to the solution, considerations must be made regarding the expectations of the system. Two counter-rotating cylinders with fluid adhering to the surfaces are present. It is therefore expected that the fluid velocities at these surfaces correspond to the tangential velocities of the cylinders, i.e., $$u_{\varphi,i}=R_i\Omega_i,\ i=[1,2]$$. Since the cylinders rotate in opposite directions, there will be a change of direction in the velocity distribution.
1.5.1 Velocity
Equation (3) analytically describes the velocity field of the fluid. After substituting the values, the field can be represented graphically. Since we are interested not only in the analytical closed-form solution but also in the numerical approximate solution, we display these in various mesh qualities. For better representation of the results, the velocity field is not shown in its individual components but as a magnitude. While this eliminates any information about the velocity direction at first glance, it proves useful in the subsequent representation.
The analytical solution shows non-zero velocities at the surfaces of the cylinders and a region in between where the velocity actually assumes the value zero. At this point, the change of direction of the velocity occurs. The expectations are thus fulfilled.
In the first superficial examination, the images differ only slightly from each other. What can be recognized without difficulty is that the results with the default mesh are merely less smooth but generally represent the analytical solution. Unfortunately, the numerically generated “colored images” alone are not always meaningful, and numerical values are needed at this point to establish a serious comparison between analytical and numerical results. To analyze the results in more detail, a path is laid through the fluid that represents the physical properties in the fluid. The L2-norm of the velocity is plotted graphically over the radius (from the start of the path at R1 to the end of the path at R2):
The tangential velocities of the fluid should, as expected, correspond at the beginning and end of the plot to the tangential velocities of the cylinders at R1 and R2:
Our expectations are completely fulfilled.
When turning to the numerical results of the CFD analysis, it can be recognized more clearly how large the simulation errors are with poor meshes. Only the good and structured mesh reproduces the analytical results of the calculation with sufficient accuracy.
1.5.2 Pressure
As with the velocity distribution, the velocity field is now determined. Analytically, equation (4) describes the corresponding pressure distribution in the fluid. The distributions are graphically represented below.
The deviation of the numerical results from the analytical solution is also clearly visible in these figures. The influence of the mesh is much greater than in the velocity field.
Here too, the pressure is evaluated over the same path as in the velocity:
Here too, mesh quality plays an enormously important role in the results, and the error with poor meshes is greater than can be seen in the velocity field.
Even with seemingly simple problems, it is surprisingly easy to fall into the trap of errors and inaccuracies, and the resulting calculations can be of questionable quality. The validation of numerical models is therefore of critical importance. Without solid validation, whether through analytical methods or experimental data, there is a risk of delivering worthless or incorrect results.
Modern Computational Fluid Dynamics (CFD) programs in particular often suggest an almost effortless solution to flow problems. With their default settings, they relieve the user of a considerable part of the work. However, there is a potential danger here: the generated result is quickly available but not necessarily correct. Often, the simulation results deviate significantly from reality.
This phenomenon becomes particularly evident in the present simple studies in a two-dimensional space with incompressible, steady-state, and laminar flow conditions. Although the numerical solution appears acceptable at first glance, errors can occur even under seemingly trivial conditions that severely impair the result. This underscores the necessity to carefully validate numerical models and understand the limits of their applicability.
If you do not consider yourself a proficient user of numerical methods or if analytical methods seem too complicated and abstract to you, it is advisable to consult specialists. These specialists have the necessary expertise and experience to solve complex problems from industry and research and achieve accurate results.
$$v_Θ (r)=A⋅r+B/r$$