Nomenclature

$\omega$

vorticity (1/s)

n

normal vector

$\rho$

density (kg/m³)

$\mu$

dynamics viscosity coefficient (Pa·s)

T

total pressure ratio

P ^*

total pressure (Pa)

T ^*

total temperature (K)

$\dot m$

rate of flow (kg/s)

η ^*

efficiency

π

total pressure ratio

BVF

boundary vorticity flux (m/s ²)

BEF

boundary enstrophy flux (kg/s³)

$\omega_{\theta}$

bircumferential vorticity (1/s)

2.0 Vorticity dynamics theories

BVF is considered the root of multi-scale flow separation in the passage [Reference Lu, Li and Zuo31–Reference Peng, Ouyang and Wu33]. The essence of BVF diagnosis method is to express the forces and moments acting on the blade as the integral form of the first or second moment of BVF, by observing the distribution characteristics of BVF on the blade surface, the fluid’s positive or negative contribution regions to the blade moment can be directly determined.

The concept of BVF was first proposed by Lighthill [Reference Lighthill34]; this parameter represents the rate at which the vorticity diffuses into the fluid through viscous diffusion, and it describes the vorticity flux diffused through a unit area in unit time. It is an indispensable basic tool for studying the interaction between vorticity and wall. The equations are as follows:

(1) \begin{align}\omega = \nabla \times u\end{align}

(2) \begin{align}\sigma = {{\partial (\mu \omega )} \over {\partial n}}\end{align}

Where $\omega$ is the vorticity, $\mu$ is the dynamics viscosity coefficient, n is the unit normal vector for wall. For three-dimensional viscous compressible flow, BVF consists of three parts:

(3) \begin{align}\sigma = {\sigma _a} + {\sigma _p} + {\sigma _\tau }\end{align}

It refers to $\sigma_{a}$ produced by inertia force, $\sigma_{p}$ produced by normal stress and $\sigma_{\tau}$ produced by shear stress on the surface of the fluid, respectively.

(4) \begin{align}{\sigma _a} = n \times a\end{align}

(5) \begin{align}{\sigma _p} = \frac{1}{\rho }n \times \nabla p\end{align}

(6) \begin{align}{\sigma _\tau } = \frac{1}{\rho }(n \times \nabla ) \times (\mu \omega )\end{align}

Where p is the pressure, $\rho$ is the density.

The flow in the high-load counter-rotating compressor passage is dominated by high Reynolds number flow, and the order of magnitude of $\sigma_{\tau}$ is much smaller than that of $\sigma_{p}$ , which can be basically ignored in the flow analysis. The blade is in a uniform rotation state, and there is wall no slip, so BVF only needs to consider $\sigma_{p}$ generated by normal stress. For the counter-rotating compressor, the normal stress is the pressure. The Equation (7) is as follows:

(7) \begin{align}\sigma = {\sigma _p} = \left| \begin{array}{c@{\quad}c@{\quad}c} x & y & z \\ i & j & k \\ {{{\partial p} \mathord{\left/ {\vphantom {{\partial p} {\partial x}}} \right. } {\partial x}}} & {{{\partial p} \mathord{\left/ {\vphantom {{\partial p} {\partial y}}} \right. } {\partial y}}} & {{{\partial p} \mathord{\left/ {\vphantom {{\partial p} {\partial z}}} \right. } {\partial z}}} \end{array} \right|\end{align}

The axial moment generated by the pressure gradient on the blade surface can be expressed as follows:

(8) \begin{align}{M_Z} = - {\int_S {(r \times pn)} _z}dS = - \frac{1}{2}\int_S {\rho (r \cdot r)} {\sigma _p}_zdS + \frac{1}{2}\oint_{\partial S} {\rho (r \cdot r)} dz\end{align}

Where s is the surface area of blade, ∂S is the boundary of blade surface, $\sigma_{pz}$ is the axial component of $\sigma_{p}$ .

It can be seen from Equation (8) that the first term on the right is the surface integral of $\sigma_{pz}$ on the blade, and the second term on the right is the line integral of the pressure moment on the boundary. The $\sigma_{pz}$ is shown in Equation (9), which is mainly related to the radial and circumferential pressure gradients and is the main reason for the formation of the passage vortex and radial secondary flow [Reference Wu and Jiang35]. The axial moment is closely related to the net power of the compressor, Equation (10) gives the relationship between the net power and $\sigma_{pz}$ on the blade surface, and it can be analysed that $\sigma_{pz}$ positive peak will be detrimental to the net power. Therefore, BVF diagnosis only needs to focus on the distribution characteristics of $\sigma_{pz}$ on the blade surface. For the convenience of writing, the BVF diagnostic method mentioned below diagnoses the distribution characteristics of $\sigma_{pz}$ on the blade surface.

(9) \begin{align}{\sigma _p}_z = \frac{1}{\rho }\left( {\frac{{\partial p}}{{\partial y}}i - \frac{{\partial p}}{{\partial x}}j} \right)\end{align}

(10) \begin{align}{L_u} = - \Omega {M_z} = \Omega \left( { - \frac{1}{2}\int_S {\rho (r \cdot r)} {\sigma _p}_zdS + \frac{1}{2}\oint_{\partial S} {\rho (r \cdot r)} dz} \right)\end{align}

The vector lines of vorticity dynamics parameters can also predict the three-dimensional boundary layer separation phenomenon. When skin-friction vector $\tau$ -lines on the blade surface converge, the BVF $\sigma_{p}$ -lines turn basically along the direction of the skin-friction vector $\tau$ -lines, that is, parallel to the skin-friction vector $\tau$ -lines, and there is a positive peak of $\sigma_{pz}$ , the flow in the nearby area will be separated. In addition, when skin-friction vector $\tau$ -lines on the blade surface converge, and the curvature of the vorticity vector $\omega$ -lines reach the local maximum value at the same time, it can also be determined that there is flow separation, which is the position that needs to be focused on optimised [Reference Zhang, Cai and Li36, Reference Wu, Ma and Zhou37]. The expression of the above theory is as follows:

(11) \begin{align}\frac{{{\sigma _p} \times \tau }}{{\left| {{\sigma _p}} \right|\left| \tau \right|}} = b({{\mathop{\rm Re}\nolimits} ^{ - 1/8}})\end{align}

(12) \begin{align}\tau = \mu \omega \times n\end{align}

Where b(Re^-1/8) is the value of Re^−1/8 at a point on the blade surface.

Equation (11) shows that the angle between BVF $\sigma_{p}$ -lines and the skin-friction vector $\tau$ -lines is inversely proportional to Reynolds number, the smaller the Reynolds number, the larger the angle.

The surface integral of BVF on the non-rotating solid wall is always zero, which makes it impossible to use BVF diagnostic method to study how much vorticity diffuses into the fluid on a closed surface. Enstrophy is another important concept in vorticity dynamics theory, which is defined as a scalar 0.5 $\omega$ ². Its changes involve dynamics properties and kinematic properties, and it can be used to judge whether the flow field is swirling and the degree of swirling. The square of the vorticity in enstrophy can avoid that $\omega$ in the closed vortex tube will offset each other when integrating. Based on the Navier-Stokes energy equation, it is analysed that enstrophy can reflect the energy dissipation characteristics of the fluid. Similar to the mathematical relationship of BVF, the normal gradient of enstrophy on the blade surface is the BEF, the co-analysis of BVF and BEF can comprehensively study the physical mechanism of vorticity entering the fluid from the boundary.

The definition formula of BEF is as follows:

(13) \begin{align}\eta * \equiv \nu \frac{\partial }{{\partial n}}\left( {\frac{1}{2}{\omega ^2}} \right) = \frac{1}{\rho }\omega \cdot \sigma \end{align}

It can be seen that BEF and the projections of BVF along the vorticity direction are closely related. It can be concluded from the above Equation (13) that in some regions of wall, when $\eta$ > 0, the direction of BVF and vorticity are the same. In this case, the shedding vorticity from the wall will increase the vorticity near the wall, and this area of wall is the vorticity source. When $\eta$ <0, the direction of BVF and vorticity are the opposite. In this case, the shedding vorticity from the wall will offset the vorticity near the wall, and this area of wall is the vorticity convergence.

Wu et al. [Reference Wu, Ma and Zhou37] summarised the distribution laws of different parameters before and after flow separation in a two-dimensional plate boundary layer, as shown in Fig. 1. For the fluid flow with an adverse pressure gradient in the boundary layer, the sign of the BVF is opposite to the vorticity near the wall, the shedding vorticity from the wall will offset the original vorticity, and the continuous development of this situation will cause $\partial$ u/ $\partial$ y, hence flow separation occurs according to Prandtl’s separation theory. It can be seen that the appearance of BVF positive peak is the sufficient condition for flow separation. In addition, when BEF goes from less than zero to greater than zero, airflow also starts to separate.

Figure 1.

Parameter distribution of boundary layer during transition from attachment to separation. (a) velocity, (b) vorticity, (c) BVF and (d) BEF [Reference Wu, Ma and Zhou37].

There is a direct mathematical relationship between the aerodynamic performance of the compressor and the circumferential vorticity, and the distribution of the circumferential vorticity directly affects the flow state in the compressor passage [Reference Wu18]. Let S ₁ , and S ₂ denote inlet plane and outlet plane for compressor, respectively. The relationship between total pressure ratio and circumferential vorticity is as follows:

(14) \begin{align}\pi = {{p_2^*} \over {p_1^*}} = {{{{\int_{{{\rm{S}}_{\rm{2}}}} {\rho {u_z}{p^*}dS} } \over {\int_{{{\rm{S}}_{\rm{2}}}} {\rho {u_z}dS} }}} \over {{{\int_{{{\rm{S}}_{\rm{1}}}} {\rho {u_z}{p^*}dS} } \over {\int_{{{\rm{S}}_{\rm{1}}}} {\rho {u_z}dS} }}}} = {{{{\int_{{{\rm{S}}_{\rm{2}}}} {\rho {u_z}{p^*}dS} } \over {{{\dot m}_2}}}} \over {{{\int_{{{\rm{S}}_{\rm{1}}}} {\rho {u_z}{p^*}dS} } \over {{{\dot m}_1}}}}} = {{\int_{{{\rm{S}}_{\rm{2}}}} {\rho {u_z}{p^*}dS} } \over {\int_{{{\rm{S}}_{\rm{1}}}} {\rho {u_z}{p^*}dS} }} = {{\int_{{{\rm{S}}_{\rm{2}}}} {\rho {u_z}{p^*}dS} } \over {{C_1}}}{\rm{ = }}{{ - \int_{{{\rm{S}}_{\rm{2}}}} {r\rho {p^*}{{\partial {u_z}} \over {\partial r}}dS} } \over {2{C_1}}} \approx {{\int_{{{\rm{S}}_{\rm{2}}}} {r\rho {p^*}{\omega _\theta }dS} } \over {{C_2}}}\end{align}

(15) \begin{align}{C_2} = 2{C_1} = \int_{{{\rm{S}}_{\rm{1}}}} {\rho {u_z}{p^*}dS} \approx - 0.5\int_{{{\rm{S}}_{\rm{1}}}} {r\rho {p^*}\frac{{\partial {u_z}}}{{\partial r}}dS} = - \pi \int_{{r_1}}^{{r_2}} {{r^2}\rho {p^*}\frac{{\partial {u_z}}}{{\partial r}}dr} \approx \pi \int_{{r_1}}^{{r_2}} {{r^2}\rho {p^*}{\omega _\theta }dr} \end{align}

where, C ₁ and C ₂ are constants at inlet plane of compressor, P ^* is the total pressure, $\omega_{\theta}$ is circumferential vorticity, u _z is the axial component of velocity, r is the radial position.

Using a similar method to the total pressure ratio, the relation of total temperature ratio by circumferential vorticity expression is written as follows:

(16) \begin{align}T = {{T_2^*} \over {T_1^*}} = {{{{\int_{{{\rm{S}}_{\rm{2}}}} {\rho {u_z}{T^*}dS} } \over {\int_{{{\rm{S}}_{\rm{2}}}} {\rho {u_z}dS} }}} \over {{{\int_{{{\rm{S}}_{\rm{1}}}} {\rho {u_z}{T^*}dS} } \over {\int_{{{\rm{S}}_{\rm{1}}}} {\rho {u_z}dS} }}}} = {{{{\int_{{{\rm{S}}_{\rm{2}}}} {\rho {u_z}{T^*}dS} } \over {{{\dot m}_2}}}} \over {{{\int_{{{\rm{S}}_{\rm{1}}}} {\rho {u_z}{T^*}dS} } \over {{{\dot m}_1}}}}} = {{\int_{{{\rm{S}}_{\rm{2}}}} {\rho {u_z}{T^*}dS} } \over {\int_{{{\rm{S}}_{\rm{1}}}} {\rho {u_z}{T^*}dS} }} = {{\int_{{{\rm{S}}_{\rm{2}}}} {\rho {u_z}{T^*}dS} } \over {{C_3}}}{\rm{ = }}{{ - \int_{{{\rm{S}}_{\rm{2}}}} {r\rho {T^*}{{\partial {u_z}} \over {\partial r}}dS} } \over {2{C_3}}} \approx {{\int_{{{\rm{S}}_{\rm{2}}}} {r\rho {T^*}{\omega _\theta }dS} } \over {{C_4}}}\end{align}

(17) \begin{align}{C_{\rm{4}}} = 2{C_{\rm{3}}} = \int_{{{\rm{S}}_{\rm{1}}}} {\rho {u_z}{T^*}dS} \approx - 0.5\int_{{{\rm{S}}_{\rm{1}}}} {r\rho {T^*}{{\partial {u_z}} \over {\partial r}}dS} = - \pi \int_{{r_1}}^{{r_2}} {{r^2}\rho {T^*}{{\partial {u_z}} \over {\partial r}}dr} \approx \pi \int_{{r_1}}^{{r_2}} {{r^2}\rho {T^*}{\omega _\theta }dr} \end{align}

Where C ₃ and C ₄ are constants at inlet plane of compressor. T ^* is the total temperature.

The compressor efficiency equation based on Equations (14–17) is as follows:

(18) \begin{align}\eta = \frac{{{{\left( {\frac{{p_2^*}}{{p_1^*}}} \right)}^{\frac{{\gamma - 1}}{\gamma }}} - 1}}{{\frac{{T_2^*}}{{T_1^*}} - 1}} = \frac{{\pi _{tt}^{\frac{{\gamma - 1}}{\gamma }} - 1}}{{{T_{tt}} - 1}} \approx \frac{{{{\left( {\frac{{\int_{{{\rm{S}}_{\rm{2}}}} {r\rho {p^*}{\omega _\theta }dS} }}{{{C_2}}}} \right)}^{\frac{{\gamma - 1}}{\gamma }}} - 1}}{{\frac{{\int_{{{\rm{S}}_{\rm{2}}}} {r\rho {T^*}{\omega _\theta }dS} }}{{{C_{\rm{4}}}}} - 1}}\end{align}

3.0 Calculation model and aerodynamic optimisation design platform

The object of diagnosis and optimisation design is a high-load counter-rotating compressor with a 1/2 + 1 aerodynamic configuration, and its geometric model is shown in Fig. 2. There are three rows of blades from the inlet to the outlet, which are low-pressure rotor (R1), high-pressure rotor (R2), and stator in turn. The three-dimensional fluid simulation software NUMECA is used for steady calculation. Among them, the Spalart-Allmaras model is selected as the turbulence model, which can accurately simulate the shock wave and separation phenomenon in the passage, and this turbulence model has been used by some scholars [Reference Wu, Ma and Zhou37–Reference Li, Wu and Guo40]. At the inlet, the total temperature is 288.15 K and the total pressure is 101325 Pa. At the outlet, the static pressure at the design point is 274 kPa. The numerical method verification are described in (Reference Yan, Chen, Fang and Yan41).

Figure 2.

Geometric model.

The geometric mesh division is completed in the NUMECA/AutoGrid5 module. The schematic diagram of meridional plane mesh is shown in Fig. 3. The O₄H grid topology is selected. In order to accurately solve the internal flow of the boundary layer, the surface of the blade profile adopts the structured grid, and the height of the first layer of the grid is 0.005 mm. Due to the complex shock wave problem in the high-Mach-number compressor channel, the degree to which the number of grids affects the simulation results must be investigated. The results of the global parameter grid independence tests are shown in Table 1. When the grid number exceeds 2.4 million, the relative changes between the mass flow rate, total pressure ratio and adiabatic efficiency are all less than 0.1% at the design point, and the simulation results tend to be stable. In addition, Fig. 4 shows the isentropic Mach number of the 50% span section of the R1 blade under different cases, and the results show that under different grid numbers, its value changes relatively little. Considering the computational cost, the case 4 grid parameter setting is adopted in all subsequent examples in this paper.

Figure 3.

The schematic diagram of the meridional plane mesh.

The aerodynamic optimisation design platform developed by the research group uses the numerical optimisation method based on an artificial neural network and genetic algorithm, and it is finally integrated on the NUMECA/Design3D module by coupling computational fluid dynamics technology and a self-edited Python script. The logical relationship is shown in Fig. 5. The user defined module can complete the functions of blade parametric modeling, grid division, physical model selection and boundary condition setting. The Python script in the multi-stage blade optimisation script module is the highlight of the highly automated platform, which first determines the blade geometric parameters that need to be modified and defines the geometric parameters as optimisation variables that the platform can recognise, and then completes the link of the geometry functions. Second, the method determines the range of optimisation variables, which are used to limit the parameter search space of the optimisation algorithm, and then interacts with the user defined module and the numerical optimisation module to complete the work of generating geometric files, perform grid generation and assembly of the meshes. The resulting grid is invoked by the solver of the numerical optimisation module, and the objective function is automatically obtained. For the effective samples generated by the numerical optimisation module, artificial neural networks can be used to find the relationship between the objective function and the optimisation variables, and then the optimal value of the objective function can be found through the genetic algorithm to obtain the optimal blade geometry.

Table 1.

Grid independence verification

The parameterisation of camber line is shown in Fig. 6. To keep the leading-edge position of profile unchanged, the camber line is fitted with a high-order Bezier curve controlled by five discrete points, which can ensure a smooth blade profile and avoid large aerodynamic loss. Each blade has four parametric cross-sections equidistantly distributed along the spanwise direction. They are represented by the red dotted line in Fig. 6. The specific optimisation method is to add a perturbation amount to each discrete point of the camber line of the R1 and R2 blade designated sections, respectively, and optimise the three-dimensional blade channel geometry by controlling the shape of the camber line. The purpose of optimisation is to improve the compressor efficiency on the premise of ensuring the total pressure ratio. The focus of this paper is to study the application of vorticity dynamics theory in the flow field diagnosis and aerodynamic performance optimisation design of counter-rotating compressor. Therefore, the specific content of optimisation design will not be described too much. The optimisation process and results are described in (Reference Yan, Chen, Fang and Yan41).

Figure 4.

Isentropic Mach number.

Figure 5.

Aerodynamic optimisation design platform.