Proceedings of the 8th International Conference on Heat Transfer and Fluid Flow (HTFF'21)
Prague, Czech Republic Virtual Conference – August, 2021
Paper No. HTFF 126
DOI: 10.11159/htff21.126

Simulations of Low Mach number Flows and Heat Transfer
in Rotating Systems
Ahmed Hodaib*, Isabelle Raspo*, Stéphane Viazzo*, Anthony Randriamampianina*
Aix Marseille Univ, CNRS, Centrale Marseille, M2P2, Marseille, France*

Abstract – Low Mach number flows can be observed in internal air systems of turbomachinery aeroengines, in particular the secondary
air flow extracted from the compressor stage for cooling purpose of the turbine disks. A parallelized domain decomposition technique isproposed for the simulation of three-dimensional flows and heat transfer in a T-shape model rotating cavity, relevant to study this flow.
The Navier-Stokes equations, coupled with energy and perfect gas equations, are solved in the framework of the Low Mach Numberapproximation, allowing a reduction of computational costs by filtering the high-speed sound waves while keeping a good accuracy byconsidering the compressibility effects. A higher-order solver based on a fourth-order compact spatial discretization associ- ated with aparallelized Fourier method is implemented on a staggered grid. A semi-implicit second order scheme is introduced for time integration.
The resulting computational code is parallelized using a hybrid MPI/OpenMP approach. Simulations were performed within an idealized
HP compressor inter-disk cavity, while considering conductive exchanges between the walls and the fluid.

Keywords: Rotating flows, Buoyancy-induced flows, Low Mach number approximation, Domain decomposition, Higherorder compact schemes, Parallel computing.

1. Introduction

Low Mach number flows are involved in many rotating systems. In aircraft engines, we can cite the secondary air flow
passing through the high-pressure HP compressor cavities for cooling purpose of the turbine disks. A better understandingof this complex buoyancy-induced flow is essential to determine the thermal stresses, the radial growth of the blades and thetemperature rise of the air used for cooling, in order to improve the engine performance. The induced convective flow is not
only unsteady and three-dimensional, but also unstable. Moreover, the problem is strongly conjugate: the flow affects thetemperature of the disks, and vice versa [1]. Numerical simulation of Low Mach number flows is not straightforward because
algorithms commonly used for compressible flows either fail to converge or require severe time step restrictions due to thedifference of scales between the characteristic flow velocity and the speed of sound waves. In addition, the solution of thediscrete equations contains pressure fluctuations at the order of the Mach number Ma whereas only pressure fluctuations at
the order of Ma² are present in the continuous equations, which can induce a loss of accuracy [2]. These problems can beovercome using the Low Mach Number (LMN) approximation proposed by Paolucci [3]. Indeed, this asymptotic
development with respect to Ma² enables the filtering of acoustic waves and, consequently, removes the induced severe time
step restrictions. The present paper proposes a higher-order solver based on this approximation combined with a pressurebased method for the simulation of rotating flows subject to high temperature difference. The approach is based on thedecoupled algorithm developed by [4] and improved here by solving a variable-coefficient Poisson equation for the pressure
to enforce the velocity divergence constraint, as recommended by Nicoud [5] and Motheau and Abraham [6]. To handle
semi-complex geometries, as considered in the present work (Fig. 1), a parallelized domain decomposition technique basedon patching-collocation technique is implemented [7]. This new Low Mach number solver is applied to the study ofbuoyancy-induced flows within an idealized High-Pressure compressor inter-disk cavity. Results are validated with respect
to available experimental data [8].

2. Mathematical ModelWe consider 3D non-isothermal flows of air in a T-shape cavity composed of three co-axial cylinders of respective radiirs, a and b (see Fig. 1). This geometry is a simplified model of a HP compressor inter-disk cavity. The whole cavity rotatesaround the horizontal z-axis with an angular velocity Ω, while the inner shaft at radius rs is kept stationary, in accordancewith the experimental setup used by Atkins and Kanjirakkad [8]. An axial throughflow of cooling air is injected between thetwo inner cylinders at radii rs and a, while the outer cylinder at radius b is maintained at a hot temperature. The governing
equations are the compressible Navier-Stokes equations coupled with energy and perfect gas equations in the Low MachHTFF 126-1

Number approximation [3]. Thus, the pressure is split into two parts. The leading order pressure, called thermodynamic
pressure Pth, only depends on time and appears in the energy equation and in the equation of state. In the case of opensystems, it also remains constant in time. The second order pressure, called dynamic pressure Pdyn, is only involved in themomentum equations. Therefore, the Low Mach number equations in the rotating reference frame are:𝜕𝜕𝜕𝜕+ 𝜵𝜵. (𝜌𝜌𝑽𝑽) = 0(1)𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝑽𝑽11𝜌𝜌+ 𝜵𝜵. (𝜌𝜌𝑽𝑽)𝑽𝑽 + 2𝜌𝜌𝒆𝒆𝑧𝑧 × 𝑽𝑽 = −𝜵𝜵𝑃𝑃𝑑𝑑𝑑𝑑𝑑𝑑 +𝚫𝚫𝑽𝑽 +𝜵𝜵(𝜵𝜵. 𝑽𝑽) + 𝜌𝜌𝜌𝜌𝒆𝒆𝑟𝑟 − 2 (sin(𝜃𝜃) 𝒆𝒆𝑟𝑟 + cos(𝜃𝜃) 𝒆𝒆𝜃𝜃 )𝜕𝜕𝜕𝜕𝑅𝑅𝑒𝑒𝜙𝜙
3𝑅𝑅𝑒𝑒𝜙𝜙𝐹𝐹𝑟𝑟
𝜌𝜌𝛾𝛾𝜕𝜕𝑇𝑇+ 𝜌𝜌𝑽𝑽. 𝜵𝜵𝑇𝑇 = −(𝛾𝛾 − 1)𝑃𝑃𝑡𝑡ℎ 𝜵𝜵. 𝑽𝑽 +𝜵𝜵. (𝜆𝜆𝜵𝜵𝑇𝑇)
𝑃𝑃𝑃𝑃𝑃𝑃𝑒𝑒𝜙𝜙𝜕𝜕𝜕𝜕𝑃𝑃𝑡𝑡ℎ = 𝜌𝜌𝜌𝜌

(2)

(3)

(4)

where Pr is the Prandtl number, Reϕ = ρ0Ωb²/η0 is the rotational Reynolds number, Fr = Ωb/�𝑔𝑔𝑔𝑔 is the Froude numberand γ =1.4 is the specific heat ratio. The above dimensionless equations were obtained with b, Ωb and Ω-1 as characteristic
length, velocity and time, respectively, while the other variables are scaled by their initial value.
The velocity divergence can be directly calculated from thermodynamic variables by using continuity and energy
equations combined with the total differential of the equation of state:
𝛾𝛾𝑑𝑑𝑃𝑃− 𝑡𝑡ℎ + 𝑃𝑃𝑃𝑃𝑃𝑃𝑒𝑒 𝜵𝜵. (𝜆𝜆𝜵𝜵𝑇𝑇)
𝑑𝑑𝑑𝑑𝜙𝜙(5)𝜵𝜵. 𝑽𝑽 =𝛾𝛾𝑃𝑃𝑡𝑡ℎFinally, in order to establish a more accurate thermal balance, the present numerical investigation takes into accountconductive exchanges between the solid parts (made of Ti-6Al-4V alloy [8]) and the fluid by solving Fourier heat conduction
within the solid while imposing continuity of heat flux at the interface solid-fluid.

3. Numerical Method

The time discretization of Eqs. (1-3) is done through a semi-implicit second order scheme based on a backward Euler
formula for the time derivative and a second-order Adams-Bashforth evaluation for the non-linear terms, the Coriolis terms
and the crossed derivatives in the diffusive terms resulting from the use of cylindrical coordinates. For space approximation,
we implemented a fourth-order finite difference compact scheme [9] in the radial and axial directions associated with aFourier method in the azimuthal direction and defined on a staggered grid. The mesh is refined near boundaries by applying
a hyperbolic tangent stretching function in order to allow an accurate description of the thin boundary layers developing
along the rotating walls.
The discretized governing equations are solved by a two-step algorithm adapted from the one proposed by [4]. The firststep is dedicated to the computation of the thermodynamic variables (temperature and density) through an iterative process.
The energy equation is first solved as a variable-coefficient Helmholtz equation for temperature 𝑇𝑇 𝑛𝑛+1 . The density andthermal conductivity are split into a constant part and a time-dependent part in order to obtain a constant-coefficient equationthat can be solved using full diagonalization technique [9]. Then the density ρn+1 is calculated from the up-to-datetemperature using the perfect gas equation. When convergence is achieved, the velocity divergence (∇.V)n+1 is computedfrom Eq. (5). In the second step, the variable-density Navier-Stokes equations are solved by a projection-type method derived
from the one proposed by Chorin ([10]) and Teman ([11]) for incompressible flows. A predicted velocity field 𝑉𝑉� 𝑛𝑛+1 is firstcomputed by only taking into account the convective and diffusive terms in the momentum equations, which leads to avariable-coefficient Helmholtz equation, solved with a similar iterative process to the one implemented for the temperatureequation. This predicted velocity is then corrected to satisfy the velocity divergence constraint. This constraint is enforced
through a variable-coefficient Poisson problem for the dynamic pressure, solved iteratively as in [6].
The main drawback of compact schemes stems from their restriction to simple rectangular domains. In order to allow
their use in more complex geometries, we implemented the domain decomposition technique developed by [7]. It is aHTFF 126-2

nonoverlapping patching-collocation method based on an iterative algorithm to enforce the continuity conditions, namelycontinuity of the variables and of their first order normal derivative, at the interface between two adjacent subdomains. The
multidomain approach was implemented using MPI parallelization. In addition, the Fast Fourier Transforms FFT, introducedin the azimuthal direction, were parallelized using the OpenMP library.

4. Results

The configuration of interest is recognized to be one of the most challenging rotating systems for computational fluiddynamics. From the fundamental viewpoint, motion can develop under the concurrent action of three generating mechanisms:
the forced throughflow, the thermal gradient and the rotation. Moreover, thermal gradients may simultaneously arise fromconduction, convection and radiation. Each main mechanism can alternately drive the fluid flow and the heat transfer in very
different ways, which still makes difficult a complete understanding of the complex phenomena under consideration. Owenand Long [1] mentioned that when the central flow is hotter than the shroud, as during acceleration or deceleration of engines,
the flow can be stably stratified. On the other hand, when the shroud is hotter than the axial throughflow, as during cruise
condition, it gives rise to buoyancy-induced flows. However, buoyancy-induced flow, characterized by the occurrence ofcyclonic-anticyclonic vortices resulting from rotation and radial thermal gradient, can coexist with a recirculating vortex
generated by the axial throughflow, depending on the values of parameters.

Fig. 1: Geometrical configuration

Fig. 2: Velocity vector field in the mid-cavity meridional plane(a)(b)

Fig. 3: (a) Normalized temperature in the mid-cavity meridional plane (b) Radial profile of azimuthally-averaged local Nusselt numberHTFF 126-3

We have imposed a large temperature difference between the inlet and the shroud defined by a value of the buoyancy
parameter βΔT=0.29 - outside the domain of validity of the Boussinesq approximation - and a value of Grashof number: Gr=8.9×108 (= (1-a/b)3Reϕ2βΔT). We have considered a value of Rossby number Ro=0.578 (=Wz/Ωa) - the ratio of inertial toCoriolis forces in rotating flows -, where Wz is the mean axial velocity component of the throughflow, and for the axial androtational Reynolds numbers Rez = 3200 (= 2ρ0(a-rs)Wz/η0) and Reϕ = 105. This value of Ro was used by [8] but withhigher values of Rez, Reϕ and Gr.
We have displayed in figure 2 an instantaneous velocity field at the inter-disk mid-plane. The cavity is rotating
anticlockwise. Towards the upper part, we observe the formation of a pair of cyclonic (right) and anticyclonic (left) vortices
generated by a radial inflow coming from the Ekman layer, consistent with experimental evidence [12]. In Fig. 3(a), thecorresponding normalized temperature field is presented. The hot temperature from the shroud has reached the exit of thecavity, as seen in the upper part accordingly to the flow structure, while cold fluid remains in the lower part. A temperaturerise of 9% at the outflow was calculated. In Fig. 3(b), the radial variation of the azimuthally-averaged Nusselt number at thedownstream disk shows the characteristic increase with the radius, with a slight variation towards the hot shroud resulting
from the presence of a thermal plume induced by the centrifugal buoyancy (with a value of Fr=3.92, giving a value ofcentrifugal acceleration Ω²b about 15 times the gravity). This thermal plume is noticeable on the temperature field in a
meridional plane taken in a horizontal section perpendicular to the inter-disk plane considered above. It is worth mentioning
from inspection of Eq. (2) that the centrifugal acceleration is opposed to the gravity in the upper part, unlike in the lowerpart.

5. Conclusion

The purpose of this work was first to check the feasibility of the LMN approach in this specific complex configuration,
before extending our simulations to representative conditions as reported in the experimental studies of Atkins andKanjirakkad [8]. For the values of considered governing parameters, the numerical code was able to capture the features offlow structure and heat transfer in agreement with the behaviors reported in the literature. Computations at higher values ofparameters representative of operating conditions are in progress.

Acknowledgements

This work was granted access to the HPC resources under the allocation 2020-A082A01444 made by GENCI. Theauthors are grateful to the team of IDRIS France for their help on the advanced optimization of the numerical tool.

References[1] Owen J.M., Long C.A., Review of buoyancy-induced flow in rotating cavities, J. Turbomach., 137, 111001, 2015.
[2] Guillard H., Viozat C., On the behaviour of upwind schemes in the low Mach number limit, Comput. Fluids, 28, 63–86,
1999.
[3] Paolucci S., On the filtering of sound from the Navier-Stokes equations, Technical report, Sandia National Lab., USA,
SAND-82-8257, 1982.
[4] Raspo I., Ouazzani J., Un algorithme faible nombre de Mach pour la simulation des écoulements de fluides supercritiques
par des méthodes spectrales, In: Actes du XIXème Congrès de Mécanique, Ed. C Rey, P. Bontoux, A. Chrysochoos,
ISSN 2103-6225, 2009.
[5] Nicoud F., Conservative high-order finite diﬀerence schemes for low-Mach number flows, J. Comput. Phys., 158, 71-97,
2000.
[6] Motheau E., Abraham J., A high-order numerical algorithm for DNS of low-Mach-number reactive flows with detailed
chemistry and quasi-spectral accuracy, J. Comput. Phys., 313, 430-454, 2016.
[7] Louchart O., Randriamampianina A., A spectral iterative domain decomposition technique for the incompressible
Navier–Stokes equations, Appl. Numer. Math., 33, 233-240, 2000.
[8] Atkins N.R., Kanjirakkad V., Flow in a rotating cavity with axial throughflow at engine representative conditions, In
Proceedings of the ASME Turbo Expo 2014: Turbine Technical Conference and Exposition. Volume 5C: Heat Transfer.
Düsseldorf, Germany. June 16–20, 2014. V05CT16A041. ASME.
[9] Abide S., Viazzo S., A 2D compact fourth-order projection decomposition method, J. Comput. Phys., 206, 252-276,2005.
HTFF 125-4

[10] Chorin J., Numerical solution of the Navier–Stokes equations. Math Comput., 22, 745-762, 1968.
[11] Temam R., Sur l’approximation de la solution des équations de Navier–Stokes par la méthode des pas fractionnaires (I),Arch. Rational Mech. Anal., 32, 135-153, 1969.
[12] Farthing, P. R., Long, C. A., Owen, J. M., and Pincombe, J. R., Rotating Cavity With Axial Throughflow of Cooling
Air: Flow Structure, ASME J. Turbomach., 114(1), 237-246, 1992.

HTFF 125-5