asymmetry as shown in figure 1. Obviously, to control sound and vibration levels, the number of impeller blades and the number of diffuser vanes have no common divisor. Therefore, for an unsteady solution, the computational domain must cover all the fluid domain. In order to better approximate the near-wall viscous characteristics, the *k* − *ω* SST (Shear Stress Transport) model was adopted. The model treatment close to the wall combines a correction for high and low Reynolds number to predict separation on smooth surfaces, see Menter (1993) paper. Typically, this model gives a realistic estimation of the generation of the turbulent kinetic energy at the stagnation points. The SST model performance has been studied in a large number of cases. The model was rated the most accurate for aerodynamic applications in the NASA Technical Memorandum written by Bardina et al. (1997). As the geometry of the fan is complex, a hybrid mesh is used, composed of tetrahedral elements for the impeller, the diffuser and return channel volumes, and composed of hexahedral elements for upstream and downstream volumes of fluid. A previous study of Khelladi et al. (2005a) shows that a grid of 4.4 × 106 meshes is considered to be sufficiently reliable to make the numerical modeling results independent of the mesh size.

#### **2.3 Measurement methodology**

Measurements were carried out on a test bench (see figure 2) equipped with an airtight box (0.6 x 0.6 x 0.6 m), placed upstream of the centrifugal fan. The flow rate from 12 to 60 l/s is controled by changing the diameter of a diaphragm orifice. The unsteady aerodynamic pressure was measured at various positions within the impeller-diffuser unit. Kulite type dynamic sensors, with a diameter of 1.6 mm and a band-width of 125 kHz were used. They are placed at the impeller inlet (figure 2: A), the impeller-diffuser interface (figure 2: B), the diffuser centreline (figure 16:) and the return channel outlet (figure 2: C). This allows the measurement of a static pressure up to 140 mbar and a fluctuating component up to 194 dBA. These aerodynamic data are transmitted to a digital oscilloscope (Gould Nicolet: Sigma 90) with 8 simultaneous channels whose band-width is 25 MHz and has a resolution of 12 bits. Figure 3 presents a schematic of the measuring equipment.

Fig. 2. Diagram of the test bench

#### **2.3.1 Overall performances and flow field description**

4 Will-be-set-by-IN-TECH

asymmetry as shown in figure 1. Obviously, to control sound and vibration levels, the number of impeller blades and the number of diffuser vanes have no common divisor. Therefore, for an unsteady solution, the computational domain must cover all the fluid domain. In order to better approximate the near-wall viscous characteristics, the *k* − *ω* SST (Shear Stress Transport) model was adopted. The model treatment close to the wall combines a correction for high and low Reynolds number to predict separation on smooth surfaces, see Menter (1993) paper. Typically, this model gives a realistic estimation of the generation of the turbulent kinetic energy at the stagnation points. The SST model performance has been studied in a large number of cases. The model was rated the most accurate for aerodynamic applications in the NASA Technical Memorandum written by Bardina et al. (1997). As the geometry of the fan is complex, a hybrid mesh is used, composed of tetrahedral elements for the impeller, the diffuser and return channel volumes, and composed of hexahedral elements for upstream and downstream volumes of fluid. A previous study of Khelladi et al. (2005a) shows that a grid of 4.4 × 106 meshes is considered to be sufficiently reliable to make the numerical modeling

Measurements were carried out on a test bench (see figure 2) equipped with an airtight box (0.6 x 0.6 x 0.6 m), placed upstream of the centrifugal fan. The flow rate from 12 to 60 l/s is controled by changing the diameter of a diaphragm orifice. The unsteady aerodynamic pressure was measured at various positions within the impeller-diffuser unit. Kulite type dynamic sensors, with a diameter of 1.6 mm and a band-width of 125 kHz were used. They are placed at the impeller inlet (figure 2: A), the impeller-diffuser interface (figure 2: B), the diffuser centreline (figure 16:) and the return channel outlet (figure 2: C). This allows the measurement of a static pressure up to 140 mbar and a fluctuating component up to 194 dBA. These aerodynamic data are transmitted to a digital oscilloscope (Gould Nicolet: Sigma 90) with 8 simultaneous channels whose band-width is 25 MHz and has a resolution of 12 bits.

results independent of the mesh size.

Figure 3 presents a schematic of the measuring equipment.

**2.3 Measurement methodology**

Fig. 2. Diagram of the test bench

Experimental measurements were conducted in order to get the overall and local air-flow characteristics. Two numerical simulations were carried out, the first one does not take into account the cavity between the impeller shroud and the casing, while the second one does. The simulation results relating to the pressure fluctuations in various parts of the fan are compared with the experimental results. For all pressure measurements the uncertainty is ±2.5 mbar.

Figure 3 represents the evolution of the pressure at the inlet of the impeller and the outlet of the return channel according to the flow rate mesured with an uncertainty of ±0.4 l/s. At the outlet of the return channel measurements and numerical results are in good agreement. A phenomenon of blocking produce the increase in static pressure with the flow rate. At the impeller inlet, the calculated pressure curve corresponding to a case without cavity between the impeller and the casing, predicts obviously more depression compared to tests and case with cavity, particularly with flow rates from 21 to 40 l/s. The pressure rise difference between cases with and without gap is due to the presence of strong bakflow vortex at the impeller entrance which leads to the reduction in total and static pressure rise (pressure difference between the inlet and the outlet), see Khelladi et al. (2005a). The effect of vortices decreases at partial flows. The pressure curve calculated taking into account the cavity is in good agreement with the experimental results particularly with flow rates greater or in the vincinity of the nominal flow-rate of 35 l/s. When the fan works at partial flow rates the solution is more difficult to achieve accurately, specially because of the presence of separated flows.

At the nominal flow-rate, the difference of static pressure between experimental result, the case without cavity (Case I) and with cavity (Case II) can be respectively assessed to be of 40% and 3% (experimental results are used as reference). Therefore, it can be considered that the numerical model with the cavity is in adequate accuracy with the measurement and that the effect of the cavity on the overall performance can not be ignored.

The modeling of the cavity, specially the meshing procedure is time and resources consuming. Therefore, when the computation aim the determination of global performance of the fan, volumetric losses due to the cavity and more precisely due to the clearance can be separately modeled using the results of the numerical case without cavity.

The first step is to evaluate the pressure *pf* prevailing at the ring seals in view to obtain the volumetric losses and the mass flow rate shift due to the backflow through the cavity. The model is based on simple physical considerations traducing radial balance of the force acting on fluid particles in the volume between the rotating shroud and the fix external casing. For simplification, each part is modeled by a disc. Then, let consider the elementary volume of a ring of fluid of internal radius r, external radius(r+dr) and of lateral length b. This volume is submitted to pressure and centrifugal forces respectively 2*π* · *r* · *b* · *p* for inward force, 2*π* · *b* · (*p* + *dp*) for outward (2*π* · *dr* · *b*(*p* + *dp*) is neglectible) and 2*π* · *r* · *dr* · *b* · *ρ* · *ω*2 · *r* for the centrifugal force. Equilibrium involve force balance to be zero in radial direction. The equilibrium equation is given by,

$$2\pi \cdot r \cdot p \cdot b + \rho \cdot 2\pi \cdot r \cdot dr \cdot \omega^2 \cdot r = (p + dp)2\pi \cdot r \cdot b \tag{1}$$

and after simplification,

$$dp = \rho \omega^2 \cdot r \cdot dr \tag{2}$$

Figure 4 represents the difference of pressure between experimental measurement, the results of full computation and the results of computation without cavity modified with the model.

Aerodynamic and Aeroacoustic Study of a High Rotational Speed Centrifugal Fan 133

Figure 5 represents the pressure fluctuation versus time, at point p, as obtained by calculations and tests at the operating point. The uncertainty in time measurement was ±2.5 × 10−6s . Notice that the measurement point is located at the midline of two diffuser blades on the impeller-diffuser interface, so that the blade to blade interaction between the impeller and the diffuser and the separation at the blades leading edge of the diffuser do not disturb the

Fig. 5. Variation of the pressure at a point on the impeller-diffuser interface versus time - flow

On the figure, one distinguishes the impeller blade passage as characterized by pressure peaks. The pressure decreases quickly just after the blade passage and increases gradually with the approach of the following blade. The peak to peak pressure gradient due to the blade impeller passing is very important (30 mbar) compared to (it will be show later) the pressure rise in the entire diffuser. A significant difference between theoretical results with and without cavity is noted at low pressure peaks. The comparison between tests and

Fig. 4. Pressure difference between measurement and predictions

measurement.

rate = 35 l/s.

The integration of the relation from the exit of the impeller to the clearance located respectively at radius *R*2 and *Rf* provides the expression:

$$p\_f = p\_2 - \frac{\rho}{8} \omega^2 \left(R\_2^2 - R\_f^2\right) \tag{3}$$

where *p*2 is the static pressure that prevails at this location on the theoretical machine. *pf* will be used to estimate the head loss due to the leakage: Δ*pf* = *pf* − *p*1. The velocity of the flow between the static and rotating disc vary from zero to *ωRi*. The effective velocity of the flow in the clearance is modeled by *kωRi*. The conventional assumption is made that the flow between the static and the rotating surface has almost the characteristics of a Couette flow and *k* ≈ 0.5. This proposal has been verified by a simple CFD calculation in an equivalent machine.

The leakage rate *qvf* through the gap is calculated by,

$$
\sigma v\_f = \mathbb{C}\_f \mathbb{S}\_j = \mathbb{C}\_f 2\pi R\_f j \tag{4}
$$

where *Cf* is the actual speed of the fluid in the gap.

The velocity *Cf* is calculated from the expression (5). The first term on the right end side is the loss due to the contraction and expansion through the gap. Factor 1.5 is empirical. The second term on the right is the loss by friction on the length of the ring where *j* is the hydraulic diameter. *L* is the axial length of the clearance and *λ* the friction loss coefficient for a smooth turbulent flow.

$$\frac{\Delta p\_f}{\rho} = 1.5 \frac{\mathbb{C}\_f^2}{2} + \lambda \frac{L}{2j} \frac{\mathbb{C}\_f^2}{2} \tag{5}$$

with,

$$
\lambda = \frac{0.316}{\Re e^{1/4}} \text{ and, } \Re e = \frac{\mathbb{C}\_f 2j}{\nu} \tag{6}
$$

The quantity Δ*pf* is subtracted from the charge of the fan computed without clearance and the flow rate is shifted of *qvf* . The final results are confronted figure (3) where it can be seen that the results are in good agreement with experiment on a wider range that the full simulation.

Fig. 3. Comparison of overall performance between calculation and measurement

6 Will-be-set-by-IN-TECH

The integration of the relation from the exit of the impeller to the clearance located respectively

(3)

8 *ω*2 *R*2 2 − *R*2 *f*

where *p*2 is the static pressure that prevails at this location on the theoretical machine. *pf* will be used to estimate the head loss due to the leakage: Δ*pf* = *pf* − *p*1. The velocity of the flow between the static and rotating disc vary from zero to *ωRi*. The effective velocity of the flow in the clearance is modeled by *kωRi*. The conventional assumption is made that the flow between the static and the rotating surface has almost the characteristics of a Couette flow and *k* ≈ 0.5. This proposal has been verified by a simple CFD calculation in an equivalent

The velocity *Cf* is calculated from the expression (5). The first term on the right end side is the loss due to the contraction and expansion through the gap. Factor 1.5 is empirical. The second term on the right is the loss by friction on the length of the ring where *j* is the hydraulic diameter. *L* is the axial length of the clearance and *λ* the friction loss coefficient for a smooth

> *C*2 *f* 2 + *λ L* 2*j*

�*e*1/4 and, �*e* = *Cf* 2*j*

The quantity Δ*pf* is subtracted from the charge of the fan computed without clearance and the flow rate is shifted of *qvf* . The final results are confronted figure (3) where it can be seen that the results are in good agreement with experiment on a wider range that the full simulation.

*C*2 *f*

Δ*pf*

*λ* = 0.316

*ρ* = 1.5

Fig. 3. Comparison of overall performance between calculation and measurement

*qvf* = *Cf Sj* = *Cf* 2*πRf j* (4)

2 (5)

*ν* (6)

*pf* = *p*2 − *ρ*

at radius *R*2 and *Rf* provides the expression:

The leakage rate *qvf* through the gap is calculated by,

where *Cf* is the actual speed of the fluid in the gap.

machine.

turbulent flow.

with,

Figure 4 represents the difference of pressure between experimental measurement, the results of full computation and the results of computation without cavity modified with the model.

Fig. 4. Pressure difference between measurement and predictions

Figure 5 represents the pressure fluctuation versus time, at point p, as obtained by calculations and tests at the operating point. The uncertainty in time measurement was ±2.5 × 10−6s . Notice that the measurement point is located at the midline of two diffuser blades on the impeller-diffuser interface, so that the blade to blade interaction between the impeller and the diffuser and the separation at the blades leading edge of the diffuser do not disturb the measurement.

Fig. 5. Variation of the pressure at a point on the impeller-diffuser interface versus time - flow rate = 35 l/s.

On the figure, one distinguishes the impeller blade passage as characterized by pressure peaks. The pressure decreases quickly just after the blade passage and increases gradually with the approach of the following blade. The peak to peak pressure gradient due to the blade impeller passing is very important (30 mbar) compared to (it will be show later) the pressure rise in the entire diffuser. A significant difference between theoretical results with and without cavity is noted at low pressure peaks. The comparison between tests and

the leading edge, in turn, it produce an acceleration regime in the opposite side. Naturally, this kind of flow increases both the static loading of the diffuser and the force fluctuation due to impeller/diffuser interaction. Not only the through flow capability but also the overall

Aerodynamic and Aeroacoustic Study of a High Rotational Speed Centrifugal Fan 135

Fig. 7. Comparison of contour of static pressure at a blade-to-blade section of mid radial span

The comparison of the static pressure at the section plotted across the rotating axis for the case with/without cavity is presented in figure 8, figure 8A) and figure 8B) are for Case I and Case II, respectively. The static pressure through the fan inlet and outlet can be roughly figured out for both cases, for Case I it is about 1.3 times bigger than that of Case II. The static pressure is gradually increased from impeller inlet to outlet, inducing a returned flow inside the cavity from the outlet side of impeller back to the inlet and generating a re-circulation flow pattern, see Khelladi et al. (2005a). As can be seen in figure 8B), the static pressure inside the cavity deacrease slowly from the impeller outlet to the gap. At the gap the static pressure of cavity is bigger than that of impeller, which produces an unavoidable source of flow re-circulation

Fig. 8. Comparison of contour of static pressure at a cross-rotating-axis section for the case,

The comparison of the total pressure at the cross-rotating-axis section is given in figure 9. The maximum total pressure rise for Case I is roughly 2 times the one of Case II. For Case I, the total pressure is uniformly increased through impeller and diffuser but high distortions are found at the cross-over bend and return channel. For Case II, it exists a non-uniformity at impeller inlet, cross-over bend and return channel. The almost constant total pressure inside the cavity

for the case, A) without tip clearance and B) with tip clearance

distortion ahead of the impeller as well as loss generation.

A) without tip clearance and B) with tip clearance

performance will be degraded due to the flow blockage existing in the diffuser.

numerical simulations validates this result. The calculation result taking into account the cavity is closer to the tests than that without. The general shape of the two calculation results are quite similar, the same curvatures are observed at the same time positions. So the axial cavity only affects the peak to peak pressure gradient.

Figure 16 represents the variation of the static pressure along the centerline of the diffuser as measured by flush mounted sensors that do not disturb the flow. The uncertainty in the curvilinear coordinate distance (*s*) is ±0.1 mm.

Fig. 6. Evolution of the static pressure in the curvilinear direction along the axial clearance of the diffuser - flow rate = 35 l/s.

The first three points (curvilinear coordinate = 0 to 0.008 m) of measurement in figure 6 are located in a recirculation zone, according to Khelladi et al. (2005a) and the difference between pressure distributions could be seen as a question of static pressure measurment. However, the results of the simulation where the gap is taken into account, is to bee seen mainly in a shift of the pressure rise characteristic. The internal mass flow through the impeller is higher than the mass flow measured at the inlet of the stage due to the leakage flow. This means a given pressure rise across the impeller is attributed to a lower mass flow (internal mass flow trough the impeller minus leakage mass flow). Just after this zone (curvilinear coordinate = 0.008 to 0.01 m) the pressure falls quickly at the leading edge of the diffuser to increase until reaching its maximum value on the outlet side of the diffuser (curvilinear coordinate = 0.01 to 0.045 m). The diffuser transforms a part of the kinetic energy to pressure energy, it increases the pressure from -100 to 30 mbar.

A blade-to-blade section was selected. This sectional view of the static pressure at the section of blade-to-blade at middle radial span (5 mm away from the hub)is shown in figure 7A) and figure 7B) for the Case I and II, respectively. The static pressure difference between the impeller inlet and diffuser outlet can be roughly assessed. The difference for Case I is about twice of the Case II. This implies that the existence of the leakage flow reduces the capability of static pressure rise. The inflow condition to the downstream component(i.e., the diffuser) depends on the flow field at impeller exit. Referring to both cases, it can be found that the flow around one of the diffuser vanes is characterized by the fact that the incidence angle is too positive introducing a strong stagnation point at the suction side near 8 Will-be-set-by-IN-TECH

numerical simulations validates this result. The calculation result taking into account the cavity is closer to the tests than that without. The general shape of the two calculation results are quite similar, the same curvatures are observed at the same time positions. So the axial

Figure 16 represents the variation of the static pressure along the centerline of the diffuser as measured by flush mounted sensors that do not disturb the flow. The uncertainty in the

Fig. 6. Evolution of the static pressure in the curvilinear direction along the axial clearance of

The first three points (curvilinear coordinate = 0 to 0.008 m) of measurement in figure 6 are located in a recirculation zone, according to Khelladi et al. (2005a) and the difference between pressure distributions could be seen as a question of static pressure measurment. However, the results of the simulation where the gap is taken into account, is to bee seen mainly in a shift of the pressure rise characteristic. The internal mass flow through the impeller is higher than the mass flow measured at the inlet of the stage due to the leakage flow. This means a given pressure rise across the impeller is attributed to a lower mass flow (internal mass flow trough the impeller minus leakage mass flow). Just after this zone (curvilinear coordinate = 0.008 to 0.01 m) the pressure falls quickly at the leading edge of the diffuser to increase until reaching its maximum value on the outlet side of the diffuser (curvilinear coordinate = 0.01 to 0.045 m). The diffuser transforms a part of the kinetic energy to pressure energy, it increases

A blade-to-blade section was selected. This sectional view of the static pressure at the section of blade-to-blade at middle radial span (5 mm away from the hub)is shown in figure 7A) and figure 7B) for the Case I and II, respectively. The static pressure difference between the impeller inlet and diffuser outlet can be roughly assessed. The difference for Case I is about twice of the Case II. This implies that the existence of the leakage flow reduces the capability of static pressure rise. The inflow condition to the downstream component(i.e., the diffuser) depends on the flow field at impeller exit. Referring to both cases, it can be found that the flow around one of the diffuser vanes is characterized by the fact that the incidence angle is too positive introducing a strong stagnation point at the suction side near

cavity only affects the peak to peak pressure gradient.

curvilinear coordinate distance (*s*) is ±0.1 mm.

the diffuser - flow rate = 35 l/s.

the pressure from -100 to 30 mbar.

the leading edge, in turn, it produce an acceleration regime in the opposite side. Naturally, this kind of flow increases both the static loading of the diffuser and the force fluctuation due to impeller/diffuser interaction. Not only the through flow capability but also the overall performance will be degraded due to the flow blockage existing in the diffuser.

Fig. 7. Comparison of contour of static pressure at a blade-to-blade section of mid radial span for the case, A) without tip clearance and B) with tip clearance

The comparison of the static pressure at the section plotted across the rotating axis for the case with/without cavity is presented in figure 8, figure 8A) and figure 8B) are for Case I and Case II, respectively. The static pressure through the fan inlet and outlet can be roughly figured out for both cases, for Case I it is about 1.3 times bigger than that of Case II. The static pressure is gradually increased from impeller inlet to outlet, inducing a returned flow inside the cavity from the outlet side of impeller back to the inlet and generating a re-circulation flow pattern, see Khelladi et al. (2005a). As can be seen in figure 8B), the static pressure inside the cavity deacrease slowly from the impeller outlet to the gap. At the gap the static pressure of cavity is bigger than that of impeller, which produces an unavoidable source of flow re-circulation distortion ahead of the impeller as well as loss generation.

Fig. 8. Comparison of contour of static pressure at a cross-rotating-axis section for the case, A) without tip clearance and B) with tip clearance

The comparison of the total pressure at the cross-rotating-axis section is given in figure 9. The maximum total pressure rise for Case I is roughly 2 times the one of Case II. For Case I, the total pressure is uniformly increased through impeller and diffuser but high distortions are found at the cross-over bend and return channel. For Case II, it exists a non-uniformity at impeller inlet, cross-over bend and return channel. The almost constant total pressure inside the cavity

Fig. 11. Schematic of the calculation domain with the flow components of the centrifugal fan

Aerodynamic and Aeroacoustic Study of a High Rotational Speed Centrifugal Fan 137

Fig. 12. Distribution of normalized vorticity amplitude along the centrifugal fan (Case I)

*Cp* = *Ps*2 − *Ps*1 *Pt*1 − *Ps*1

Where subscripts '*s*' denotes static pressure, '*t*' the total pressure, '2' the exit section of the component, and '1' the inlet section of the component. The static pressure recovery coefficient along the stationary parts of Case I is shown in figure 13. In this case, the maximum pressure recovery, around 0.53, occurs at the cross-over bend. Both *A* and *D* regions have the same

> *ω* = *Pt*2 − *Pt*1 *Pt*1 − *Ps*1

In figure (14) one can see that the main part of losses occurs in the return channel and the cross over bend. This is due to the flow separation resulting from the strong meridional curvature and by the 90 degree change of direction of the flow. The comparison between cases I and II is listed in Table 3. We can observe that for Case I the maximum loss occurs at the cross-over

(7)

(8)

The static pressure rise coefficient for the stationary compartment can be written as

(half view)

pressure recovery, around 0.15.

The loss coefficient for the stationary parts is given by

is found to generate a stagnation point near the impeller inlet. As mentioned-above, the cross-over bend and return channel are two dominant components with high loss generation.

Fig. 9. Comparison of contour of total pressure at a cross-rotating-axis section for the case, A) without tip clearance and B) with tip clearance

Figure 10 shows the solution of entropy for both cases. By this one can be quantitatively understood the loss distribution. For Case I, as shown in figure 10A), very high entropy gradients occur through impeller exit, cross-over bend and return channel. This can also be found in figure 10B) for Case II, typically, a big jump appears at the impeller inlet close to the gap where a strong stagnation point is produced. Due to the high rotating speed, the absolute velocity at impeller exit is so high that a strong swirling flow occurs at the cross-over bend. Due to the narrow frontal area at the cross-over bend and return channel, the strong flow blockage and vortex is accumulated.

Fig. 10. Comparison of contour of entropy at a cross-rotating-axis section for the case, A) without tip clearance and B) with tip clearance

#### **2.3.2 Loss generation**

In order to easily show (quantify) the loss generation for each flow component, a half view of the flow components into the centrifugal fan is schematically plotted as shown in figure 11. The impeller inlet is indicated by *A* in the figure and the diffuser by *B*. The cross-over bend is indicated by *C*. *D* region is the return channel and the downstream duct is indicated by *E*.

The vorticity amplitude for each component is normalized by that of the impeller outlet. As shown in figure 12 (Case I), the maximum vorticity amplitude occurs at the impeller due to its high rotating speed. *A*, *B* and *C* regions have equivalent ones with about half of the one of the impeller implying a strong swirling flow and is probably the cause of flow blockages and total pressure drop.

10 Will-be-set-by-IN-TECH

is found to generate a stagnation point near the impeller inlet. As mentioned-above, the cross-over bend and return channel are two dominant components with high loss generation.

Fig. 9. Comparison of contour of total pressure at a cross-rotating-axis section for the case, A)

Figure 10 shows the solution of entropy for both cases. By this one can be quantitatively understood the loss distribution. For Case I, as shown in figure 10A), very high entropy gradients occur through impeller exit, cross-over bend and return channel. This can also be found in figure 10B) for Case II, typically, a big jump appears at the impeller inlet close to the gap where a strong stagnation point is produced. Due to the high rotating speed, the absolute velocity at impeller exit is so high that a strong swirling flow occurs at the cross-over bend. Due to the narrow frontal area at the cross-over bend and return channel, the strong flow

Fig. 10. Comparison of contour of entropy at a cross-rotating-axis section for the case, A)

In order to easily show (quantify) the loss generation for each flow component, a half view of the flow components into the centrifugal fan is schematically plotted as shown in figure 11. The impeller inlet is indicated by *A* in the figure and the diffuser by *B*. The cross-over bend is indicated by *C*. *D* region is the return channel and the downstream duct is indicated by *E*. The vorticity amplitude for each component is normalized by that of the impeller outlet. As shown in figure 12 (Case I), the maximum vorticity amplitude occurs at the impeller due to its high rotating speed. *A*, *B* and *C* regions have equivalent ones with about half of the one of the impeller implying a strong swirling flow and is probably the cause of flow blockages and

without tip clearance and B) with tip clearance

without tip clearance and B) with tip clearance

**2.3.2 Loss generation**

total pressure drop.

blockage and vortex is accumulated.

Fig. 11. Schematic of the calculation domain with the flow components of the centrifugal fan (half view)

Fig. 12. Distribution of normalized vorticity amplitude along the centrifugal fan (Case I)

The static pressure rise coefficient for the stationary compartment can be written as

$$\mathcal{C}\_p = \frac{P\_{s2} - P\_{s1}}{P\_{t1} - P\_{s1}} \tag{7}$$

Where subscripts '*s*' denotes static pressure, '*t*' the total pressure, '2' the exit section of the component, and '1' the inlet section of the component. The static pressure recovery coefficient along the stationary parts of Case I is shown in figure 13. In this case, the maximum pressure recovery, around 0.53, occurs at the cross-over bend. Both *A* and *D* regions have the same pressure recovery, around 0.15.

The loss coefficient for the stationary parts is given by

$$
\omega = \frac{P\_{t2} - P\_{t1}}{P\_{t1} - P\_{s1}} \tag{8}
$$

In figure (14) one can see that the main part of losses occurs in the return channel and the cross over bend. This is due to the flow separation resulting from the strong meridional curvature and by the 90 degree change of direction of the flow. The comparison between cases I and II is listed in Table 3. We can observe that for Case I the maximum loss occurs at the cross-over

performances for both cases are listed in Table 4. It is found that the total efficiency has 1% drop when the cavity is taken into account, meanwhile, both the total-pressure ratio and

Aerodynamic and Aeroacoustic Study of a High Rotational Speed Centrifugal Fan 139

*η* (total-total efficiency) 0.87 0.86 *π* (total-total pressure ratio) 1.15 1.12 *τ* (total temperature ratio) 1.065 1.057

In this section, a rotational speed centrifugal fan with compact dimensions was studied numerically and experimentally. The computational models with/without cavity (between the shrouded impeller and the fixed outer casing) were set up for unsteady solutions. A volumetric loss model was coupled to the case without cavity and compared to the later configurations. The effects of the cavity on the calculated flow field and the loss generation were qualitatively analyzed. The comparison of fan performances between the prediction and measurement showed a favorable tendency. The result obtained by the volumetric loss model, which is easy to setup, coupled to the configuration without cavity gives as good result as the computational model with cavity. This result can be very interesting knowing that the model with cavity is more consumptive in grid size and computational effort than the case without cavity. In terms of the flow restriction obtained for the centrifugal fan, a new empirical constant accounting for the effects of a cavity on the flow restriction of the centrifugal fan system was approximated. In presence of the cavity, due to the difference of static pressure, a leakage flow due the clearance appears and is injected in the main stream at the impeller inlet. Subsequently, a stagnation point appears and changes the incidence angle

**3.1 Computational aeroacoustics for compelex geometries – application to subsonic**

The sound generation by a flow and its propagation are a matter of aerodynamics. Indeed, the conservation equations of mass and momentum govern both the flow dynamics and the resulting acoustic phenomena. However, the features of the aerodynamic flow and the sound are different. The first is convective and/or diffusive and the second is propagative with very low attenuation due to viscosity. On the other hand, aeroacoustic problems present a wider

Aeroacoustic noise optimization is the main topic of many widespread research studies of industrial interest, see Khelladi et al. (2008); Maaloum et al. (2004) works. In fact, the noise level emitted by a device could determine the success or failure of a new prototype. On the other hand turbomachines are widely found in industrial applications. In these devices the

The prediction of aerodynamic noise benefits from recent developments in numerical methods and computer science. However, despite the knowledge accumulated over the past few

Table 4. Comparison of overall performance between Case I and Case II

**Case I Case II**

temperature ratio drop too.

**2.4 Concluding remarks**

and then the flow performances as well.

range of wave-lengths than those of aerodynamic ones.

level of sound generated is a very important parameter of design.

**3. Aeroacoustic study**

**turbomachines**

Fig. 13. Static pressure recovery for stationary sections (Case I)

bend while it moves downstream to the return channel for Case II. This is due to the severe flow blockage occurring in Case II.

Fig. 14. Distribution of total pressure loss along the centrifugal fan (Case I)

Table 3. Comparison of loss distribution into the centrifugal fan components

The isentropic efficiency of the impeller is defined as

$$\eta = \frac{\left(\frac{P\_{l2}}{P\_{l1}}\right)^{\frac{\gamma - 1}{\gamma}} - 1}{\frac{T\_{l2}}{T\_{l1}} - 1} \tag{9}$$

Considering the viscous effect to the total pressure rise, the total pressure coefficient can be expressed as

$$
\pi = \pi^{\frac{\gamma - 1}{\gamma}} \exp\left(\frac{-\Delta s}{R}\right) \tag{10}
$$

where *τ* = 1 + *U*2*Cθ*2 *CpTt*1 , *U*2 is the rotating velocity at the impeller exit, *Cθ*2 is the circumferential component of the absolute velocity at the impeller exit, *Tt*1 is the total temperature at the impeller inlet. Δ*s* is the entropy increment and *R* is the gas constant. The calculated overall performances for both cases are listed in Table 4. It is found that the total efficiency has 1% drop when the cavity is taken into account, meanwhile, both the total-pressure ratio and temperature ratio drop too.

Table 4. Comparison of overall performance between Case I and Case II

#### **2.4 Concluding remarks**

12 Will-be-set-by-IN-TECH

bend while it moves downstream to the return channel for Case II. This is due to the severe

**Section Case I Case II A** 0.065 0.50 **Impeller** 0.13 0.14 **B** 0.035 0.79 **C** 0.93 0.79 **D** 0.57 0.97

Fig. 13. Static pressure recovery for stationary sections (Case I)

Fig. 14. Distribution of total pressure loss along the centrifugal fan (Case I)

Table 3. Comparison of loss distribution into the centrifugal fan components

*η* =

*π* = *τ*

 *Pt*2 *Pt*1 *γ*−1 *γ* − 1

*γ*−1 *γ* exp

*Tt*2

Considering the viscous effect to the total pressure rise, the total pressure coefficient can be

component of the absolute velocity at the impeller exit, *Tt*1 is the total temperature at the impeller inlet. Δ*s* is the entropy increment and *R* is the gas constant. The calculated overall

−Δ*s R*

*CpTt*1 , *U*2 is the rotating velocity at the impeller exit, *Cθ*2 is the circumferential

*Tt*1 − 1 (9)

(10)

The isentropic efficiency of the impeller is defined as

expressed as

where *τ* = 1 + *U*2*Cθ*2

flow blockage occurring in Case II.

In this section, a rotational speed centrifugal fan with compact dimensions was studied numerically and experimentally. The computational models with/without cavity (between the shrouded impeller and the fixed outer casing) were set up for unsteady solutions. A volumetric loss model was coupled to the case without cavity and compared to the later configurations. The effects of the cavity on the calculated flow field and the loss generation were qualitatively analyzed. The comparison of fan performances between the prediction and measurement showed a favorable tendency. The result obtained by the volumetric loss model, which is easy to setup, coupled to the configuration without cavity gives as good result as the computational model with cavity. This result can be very interesting knowing that the model with cavity is more consumptive in grid size and computational effort than the case without cavity. In terms of the flow restriction obtained for the centrifugal fan, a new empirical constant accounting for the effects of a cavity on the flow restriction of the centrifugal fan system was approximated. In presence of the cavity, due to the difference of static pressure, a leakage flow due the clearance appears and is injected in the main stream at the impeller inlet. Subsequently, a stagnation point appears and changes the incidence angle and then the flow performances as well.

#### **3. Aeroacoustic study**

#### **3.1 Computational aeroacoustics for compelex geometries – application to subsonic turbomachines**

The sound generation by a flow and its propagation are a matter of aerodynamics. Indeed, the conservation equations of mass and momentum govern both the flow dynamics and the resulting acoustic phenomena. However, the features of the aerodynamic flow and the sound are different. The first is convective and/or diffusive and the second is propagative with very low attenuation due to viscosity. On the other hand, aeroacoustic problems present a wider range of wave-lengths than those of aerodynamic ones.

Aeroacoustic noise optimization is the main topic of many widespread research studies of industrial interest, see Khelladi et al. (2008); Maaloum et al. (2004) works. In fact, the noise level emitted by a device could determine the success or failure of a new prototype. On the other hand turbomachines are widely found in industrial applications. In these devices the level of sound generated is a very important parameter of design.

The prediction of aerodynamic noise benefits from recent developments in numerical methods and computer science. However, despite the knowledge accumulated over the past few

as fan blades, helicopter rotor, compressors or turbines, etc. Thus, in these latter kinds of

Aerodynamic and Aeroacoustic Study of a High Rotational Speed Centrifugal Fan 141

The first attempt to formulate a theory about the acoustics of propellers was conducted by Lynam & Webb. (1919). They showed that the rotation of the blades of a propeller causes a periodic modulation of the fluid flow and associated acoustic disturbance. Another approach, initiated by Bryan (1920), is to study the propagation of a source point in uniform motion.

Gutin (1948) was the first to establish a theoretical formalism of a steady noise source through linear acoustics. He showed that steady aerodynamic forces correspond to dipole source distribution on the disc of a propeller. This model proves to be incomplete because, in reality, the noise emitted by rotating blades extends rather high frequencies. The sound at high

Advances in the prediction of noise from the airflow, are based on of Lighthill (1952; 1954) investigations. In his analogy, the generated noise is mathematically reduced to the study of wave propagation in a medium at rest, in which the effect of the flow is replaced by a distribution of sources. The pressure is therefore regarded as characterizing a sound field of small amplitude carried by a fluid, whose properties are uniform throughout the area at rest. The major intake of Lighthill is to include nonlinear terms expressing the noise generation by

Curle (1955) extended the Lighthill's analogy to include solid boundaries by treating them as distributions of surface loads. Subsequently, Ffowcs Williams & Hawkings (1969) (FW&H) have extended this approach by taking into account the motion of solid surfaces in the flow.

The concept of far-field, relative to the effects of the flow compressibility, concerns the propagation of acoustic waves produced by a pressure change in the propagation medium. The occurred disturbance propagates gradually by molecular excitement to the observer

Unlike the far-field, near-field includes the sound due to the fluid compressibility and another component called the aerodynamic disturbance field or pseudo-sound. It consists of all pressure fluctuations governed primarily by the incompressibility directly related to the flow.

In the common formulation of aeroacoustic analogy solution, the noise is radiated in free and far field. As strong hypotheses: reflections, diffractions, scattering as well as the confinement effects are not taken into account. These hypotheses make very easy the use aeroacoustic analogy for noise prediction of open rotors or free jets for example, but they are also its weaknesses in case of confining. In paper of Khelladi et al. (2008) it was shown that using the FW&H formulation to model the noise generated by a centrifugal fan does not match measurements because of the presence of a casing. Taking into account the sound attenuation of the casing to correct the directivity has not really solved the problem. For this kind of problems, it was therefore concluded that the aeroacoustic analogy does not obtain accurate

applications, we need to deal with flow through complex geometries.

This had as main feature the introduction of the concept of delayed time.

frequencies is a consequence of the unsteadiness of aerodynamic loads.

turbulent flow.

**Far and near-fields**

results.

located far from the source.

These fluctuations are local and are not propagative.

**3.2 Limits of aeroacoustic analogy and alternative approaches**

decades on the mechanisms of noise generation on complex systems as for example air delivery systems, the prediction of such a flow field and the resulting acoustic pressure, by numerical methods is still difficult. This is due to our inability to model the turbulent viscous flow with enough accuracy on complex geometries and to the complicated nature of flow through turbomachines. Until now, there is still no consensus about the aeroacoustic approach to adopt, and actually, it depends on the application. In the following we present a succinct description of the most commonly used approaches.

Previously to our exposition, we recall the concepts of far and near-fields. The concept of far-field, relative to the effects of the flow compressibility, concerns the propagation of acoustic waves produced by a pressure change in the propagation medium. The occurred disturbance propagates gradually by molecular excitement to the observer far from the source. Unlike the far-field, near-field includes the sound due to the fluid compressibility and another component called the aerodynamic disturbance field or pseudo-sound. It consists of all pressure fluctuations governed primarily by the incompressibility directly related to the flow. These fluctuations are local and are not propagative.

In Computational Aeroacoustics (CAA), two computational approaches are possible:

#### **Direct approach**

This approach consists of adjusting the aerodynamic numerical modeling to the acoustics requirements. In other words, it is needed to use numerical schemes adapted to the acoustic propagation, providing low-dissipation and low-dispersion. However, the complexity of implementing these schemes and the far field constraint, where the grid must extend over very large distances, greatly increases the computational costs and makes using this approach very difficult for complex geometries.

#### **Hybrid approach**

This approach can be divided into two types of modelling.

#### *Aeroacoustics of complex geometries*

Aeroacoustics is a science dealing with the sound generated either by the flow itself, as free jet turbulence or by its interaction with a moving or static surface, rigid or deformable, as fan blades, helicopter rotor, compressors or turbines, etc. Thus, in these latter kinds of applications, we need to deal with flow through complex geometries.

The first attempt to formulate a theory about the acoustics of propellers was conducted by Lynam & Webb. (1919). They showed that the rotation of the blades of a propeller causes a periodic modulation of the fluid flow and associated acoustic disturbance. Another approach, initiated by Bryan (1920), is to study the propagation of a source point in uniform motion. This had as main feature the introduction of the concept of delayed time.

Gutin (1948) was the first to establish a theoretical formalism of a steady noise source through linear acoustics. He showed that steady aerodynamic forces correspond to dipole source distribution on the disc of a propeller. This model proves to be incomplete because, in reality, the noise emitted by rotating blades extends rather high frequencies. The sound at high frequencies is a consequence of the unsteadiness of aerodynamic loads.

Advances in the prediction of noise from the airflow, are based on of Lighthill (1952; 1954) investigations. In his analogy, the generated noise is mathematically reduced to the study of wave propagation in a medium at rest, in which the effect of the flow is replaced by a distribution of sources. The pressure is therefore regarded as characterizing a sound field of small amplitude carried by a fluid, whose properties are uniform throughout the area at rest. The major intake of Lighthill is to include nonlinear terms expressing the noise generation by turbulent flow.

Curle (1955) extended the Lighthill's analogy to include solid boundaries by treating them as distributions of surface loads. Subsequently, Ffowcs Williams & Hawkings (1969) (FW&H) have extended this approach by taking into account the motion of solid surfaces in the flow.

#### **Far and near-fields**

14 Will-be-set-by-IN-TECH

decades on the mechanisms of noise generation on complex systems as for example air delivery systems, the prediction of such a flow field and the resulting acoustic pressure, by numerical methods is still difficult. This is due to our inability to model the turbulent viscous flow with enough accuracy on complex geometries and to the complicated nature of flow through turbomachines. Until now, there is still no consensus about the aeroacoustic approach to adopt, and actually, it depends on the application. In the following we present a succinct

Previously to our exposition, we recall the concepts of far and near-fields. The concept of far-field, relative to the effects of the flow compressibility, concerns the propagation of acoustic waves produced by a pressure change in the propagation medium. The occurred disturbance propagates gradually by molecular excitement to the observer far from the source. Unlike the far-field, near-field includes the sound due to the fluid compressibility and another component called the aerodynamic disturbance field or pseudo-sound. It consists of all pressure fluctuations governed primarily by the incompressibility directly related to the flow.

In Computational Aeroacoustics (CAA), two computational approaches are possible:

This approach consists of adjusting the aerodynamic numerical modeling to the acoustics requirements. In other words, it is needed to use numerical schemes adapted to the acoustic propagation, providing low-dissipation and low-dispersion. However, the complexity of implementing these schemes and the far field constraint, where the grid must extend over very large distances, greatly increases the computational costs and makes using this approach

• The first one is to use the direct approach near disturbances in which the acoustic waves are propagated over a short distance. They are then propagated using an adapted propagation operator, as Kirchhoff's equation, according to Farassat & Myers (1988), for example, to the far field. For adapted wave operator we mean a wave equation or other conservation equation system that permits an acoustic wave to propagate from a given acoustic source. However, the simulation of the flow field requires DNS or LES, and the treatment of boundary conditions must be done with utmost care to ensure an accurate transition

• The second consists of separating aerodynamics and acoustics computations. Thus, acoustic sources are given by aerodynamic calculation and propagated using wave equation (Ffowcs Williams & Hawkings (1969), Kirchhoff,...), linearized Euler equations or other approaches like LPCE presented by Moon & Seo (2006),.... The constraints and the

Aeroacoustics is a science dealing with the sound generated either by the flow itself, as free jet turbulence or by its interaction with a moving or static surface, rigid or deformable,

description of the most commonly used approaches.

These fluctuations are local and are not propagative.

This approach can be divided into two types of modelling.

computation time is considerably reduced compared to DNS.

very difficult for complex geometries.

between near and far fields.

*Aeroacoustics of complex geometries*

**Direct approach**

**Hybrid approach**

The concept of far-field, relative to the effects of the flow compressibility, concerns the propagation of acoustic waves produced by a pressure change in the propagation medium. The occurred disturbance propagates gradually by molecular excitement to the observer located far from the source.

Unlike the far-field, near-field includes the sound due to the fluid compressibility and another component called the aerodynamic disturbance field or pseudo-sound. It consists of all pressure fluctuations governed primarily by the incompressibility directly related to the flow. These fluctuations are local and are not propagative.

#### **3.2 Limits of aeroacoustic analogy and alternative approaches**

In the common formulation of aeroacoustic analogy solution, the noise is radiated in free and far field. As strong hypotheses: reflections, diffractions, scattering as well as the confinement effects are not taken into account. These hypotheses make very easy the use aeroacoustic analogy for noise prediction of open rotors or free jets for example, but they are also its weaknesses in case of confining. In paper of Khelladi et al. (2008) it was shown that using the FW&H formulation to model the noise generated by a centrifugal fan does not match measurements because of the presence of a casing. Taking into account the sound attenuation of the casing to correct the directivity has not really solved the problem. For this kind of problems, it was therefore concluded that the aeroacoustic analogy does not obtain accurate results.

**3.3 Linearised Euler equations**

being *S* a source term and

**3.4 Boundary conditions**

*U* =

⎛

*ρ ρu ρv p*

⎞

⎟⎟⎠

*H* =

(2002) presented a methodology to compute the sources.

enforced and the information leaving the domain.

on the rigid wall boundaries, an external mirror fictitious state *U*∗−.

**3.4.1 Reflecting boundary conditions**

⎜⎜⎝

written in conservative form are the following:

*∂U ∂t* +

⎛

⎜⎜⎝

*F* =

⎛

⎜⎜⎜⎜⎜⎝

*∂F ∂x* + *∂G ∂y*

*ρu*0 + *ρ*0*u p* + *ρ*0*u*0*u ρ*0*u*0*v u*0 *p* + *γp*0*u*

(*ρ*0*u* + *u*0*ρ*) *∂u*0

(*ρ*0*u* + *u*0*ρ*) *∂v*0

⎞

⎟⎟⎠

0

(*γ* − 1) *p*∇ · *ν*0 − (*γ* − 1) *ν* · ∇*p*0

*∂x*

*∂x*

where *ν* = (*u*, *v*), *ν*0 = (*u*0, *v*0), *ρ* is the density, *p* the pressure and *γ* = 1.4. Subscript 0 is referring to mean values. In case of an uniform mean flow, *H* is null. To solve these equations, we need to compute previously the acoustic sources S, by using LES or DNS. Bogey et al.

To solve the above system which is a hyperbolic system, we used FV-MLS method presented and commented by Khelladi et al. (2011). The FV-MLS is a new high order finite volume method for unstructured grids based on moving least squares approximation, it was demonstrated that this method is particularly efficient for this kind of hyperbolic systems. Concerning the numerical developments, we choose to present only the construction of boundary conditions which are one of the most critical points of this kind of problems.

In CAA the treatment of boundary conditions plays a key role according to Colonius et al. (1993), since even small spurious disturbances when the waves leave the domain can distort the acoustic field. In the following we expose our approach to the boundary conditions. For our modeling, the boundary conditions enter in the discretized equations through a proper definition of the numerical flux that can be written as *H*(*U*+,*U*∗−,*n*), where *n* points outward from the domain and *U*∗− is the external state variable. Depending on the boundary type, the construction of *U*∗− accounts for, both, the physical boundary conditions that must be

A perfectly reflecting boundary condition is easily obtained by defining, at each Gauss points

*G* =

+ (*ρ*0*v* + *v*0*ρ*) *∂u*0

+ (*ρ*0*v* + *v*0*ρ*) *∂v*0

⎛

⎜⎜⎝

*∂y*

*∂y*

In many aeroacoustic applications we can assume that problems are linear, see Colonius & Lele (2004). In those cases, it is possible to linearize the Euler equations around a (mean) stationary solution *U*0 = (*ρ*0, *u*0, *v*0, *p*0). Thus, we can write the Linearized Euler Equations

Aerodynamic and Aeroacoustic Study of a High Rotational Speed Centrifugal Fan 143

+ *H* = *S* (11)

⎞

⎟⎟⎠

(12)

(13)

*ρv*0 + *ρ*0*u ρ*0*v*0*u p* + *ρ*0*v*0*v v*0 *p* + *γp*0*v*

⎞

⎟⎟⎟⎟⎟⎠

To take into account confining effects, it is expected that LEE can give satisfactory results. In fact, with LEE one can use the same acoustic sources as FW&H for example, and in addition reflections, diffractions and scattering are naturally taken into account by an adequate choice of boundary conditions.

In fact, figure (15) shows the predicted overall sound level compared to measurements of the centrigugal fan presented above into the aerodynamic section. The difference between the overall sound level and mesurements is presented in figure (16). As shown in this figure, the difference between measurements and calculation is about 10 dB (over a maximum value of 110 dB) at the radial direction around 90◦ and 270◦, the calculation predict more noise at this direction than measurements. In addition to the assumptions stated above, the reason is that in the theoretical model all the surfaces are considered as acoustically transparent, the effect of the acoustic attenuation of the diffuser blades in the radial direction is not taken into account. In the axial direction, the difference do not exceed 5 dB by taking into account the attenuation of the casing. In this direction the fan radiates in free field.

Fig. 15. Numerical and experimental overall noise directivity - *m* = 1, *r*0 = 1 m, *ϕ* = *π*/4.

Fig. 16. Difference between directivities obtained by the measurements and the calculations (Δ*Lp* = *Lpexp* − *Lpnum* ) at *m* = 1, *r*0 = 1 m and *ϕ* = *π*/4.

To take into account confining effects, it is expected that LEE can give satisfactory results. In fact, with LEE not only one can uses the same acoustic sources as FW&H for example, but also reflexions, diffractions and scattering are naturally taken into account by an adequate choice of boundary conditions.

#### **3.3 Linearised Euler equations**

16 Will-be-set-by-IN-TECH

To take into account confining effects, it is expected that LEE can give satisfactory results. In fact, with LEE one can use the same acoustic sources as FW&H for example, and in addition reflections, diffractions and scattering are naturally taken into account by an adequate choice

In fact, figure (15) shows the predicted overall sound level compared to measurements of the centrigugal fan presented above into the aerodynamic section. The difference between the overall sound level and mesurements is presented in figure (16). As shown in this figure, the difference between measurements and calculation is about 10 dB (over a maximum value of 110 dB) at the radial direction around 90◦ and 270◦, the calculation predict more noise at this direction than measurements. In addition to the assumptions stated above, the reason is that in the theoretical model all the surfaces are considered as acoustically transparent, the effect of the acoustic attenuation of the diffuser blades in the radial direction is not taken into account. In the axial direction, the difference do not exceed 5 dB by taking into account the attenuation

Fig. 15. Numerical and experimental overall noise directivity - *m* = 1, *r*0 = 1 m, *ϕ* = *π*/4.

Fig. 16. Difference between directivities obtained by the measurements and the calculations

To take into account confining effects, it is expected that LEE can give satisfactory results. In fact, with LEE not only one can uses the same acoustic sources as FW&H for example, but also reflexions, diffractions and scattering are naturally taken into account by an adequate choice

) at *m* = 1, *r*0 = 1 m and *ϕ* = *π*/4.

of the casing. In this direction the fan radiates in free field.

of boundary conditions.

(Δ*Lp* =

*Lpexp* − *Lpnum*

of boundary conditions.

  In many aeroacoustic applications we can assume that problems are linear, see Colonius & Lele (2004). In those cases, it is possible to linearize the Euler equations around a (mean) stationary solution *U*0 = (*ρ*0, *u*0, *v*0, *p*0).  Subscript 0 is referring to mean values. In case of an uniform mean flow, *H* is null. To solve these equations, we need to compute previously the acoustic sources S, by using LES or DNS. Bogey et al. (2002) presented a methodology to compute the sources.

To solve the above system which is a hyperbolic system, we used FV-MLS method presented and commented by Khelladi et al. (2011). The FV-MLS is a new high order finite volume method for unstructured grids based on moving least squares approximation, it was demonstrated that this method is particularly efficient for this kind of hyperbolic systems. Concerning the numerical developments, we choose to present only the construction of boundary conditions which are one of the most critical points of this kind of problems.

#### **3.4 Boundary conditions**

In CAA the treatment of boundary conditions plays a key role according to Colonius et al. (1993), since even small spurious disturbances when the waves leave the domain can distort the acoustic field. In the following we expose our approach to the boundary conditions. For our modeling, the boundary conditions enter in the discretized equations through a proper definition of the numerical flux that can be written as *H*(*U*+,*U*∗−,*n*), where *n* points outward from the domain and *U*∗− is the external state variable. Depending on the boundary type, the construction of *U*∗− accounts for, both, the physical boundary conditions that must be enforced and the information leaving the domain.

#### **3.4.1 Reflecting boundary conditions**

A perfectly reflecting boundary condition is easily obtained by defining, at each Gauss points on the rigid wall boundaries, an external mirror fictitious state *U*∗−.

the Jacobien matrix of system (11) and |**P**| = **V**−1|**D**|**V**, **D** and **V** are respectively, eigenvalues

Aerodynamic and Aeroacoustic Study of a High Rotational Speed Centrifugal Fan 145

2*c*0 (−*L*1 + *L*2) −1

2 (*L*1 + *L*2 − 2*L*4) *nx*

2 (−*L*1 + *L*2) 1

*xL*4

2 (*L*1 + *L*2) + *n*2

*L*1 = |*V*0 · *n* − *c*0|

*L*3 = *L*4 = |*V*0 · *n*|

This case reproduces the example of Bailly & Juvé (2000). The radiation of a monopole source is computed in a subsonic mean flow, with Mach number *Mx* = 0.5. The source is located at

2

where the angular frequency is *ω* = 2*π*/30 and *t* is the time coordinate. The wave length is *λ* = 30 units, and the computational domain is a square with 200 units for each side.

testing the stability and the behavior of the proposed method for the boundary conditions, an unstructured grid absorbing layer has been added. The absorbing layer is placed from the boundary of the computational domain to *x* = ±300 and *y* = ±300. In figure 17 (top) it is shown the grid used for the resolution of this problem. To build this grid, 800 equally spaced nodes at the circumference of the computational domain are used and 120 nodes at the outer

In addition to the absorbing boundary condition given by equation (17), the shape filter

Two acoustic waves propagate upstream and downstream of the source, and due to the effect of the mean flow, the apparent wavelength is modified and it is different upstream (*λ*1 =

In figure 18 pressure isocontours for different non-dimensional times *t* are shown. The pressure profile along axis *y* = 0 at time *t* = 270 is reproduced in figure 19, and also match

parameters of the absorbing layer are *κx* = *κy* = 8, see Nogueira et al. (2010) paper.

A 5th order mass matrix based FV-MLS solver is used for this example.

(1 − *Mx*) *λ*) and downstream (*λ*2 = (1 + *Mx*) *λ*) of the source.

2 + (*y* − *ys*)

2

�

*ρ*0*c*0/Δ*x*, 0, 0, *ρ*0*c*3

*c*2 0 *L*3 + 1

*ny*

*L*2 = |*V*0 · *n* + *c*0| (19)

sin (*ωt*) × [1, 0, 0, 1]

0/Δ*x* �*T*

2*c*0 (*L*1 + *L*2)

⎞

⎟⎟⎟⎟⎟⎠

*T* (20)

. With the aim of

(18)

2*c*0 (−*L*1 + *L*2)

2*c*0 (−*L*1 + *L*2)

2 (*L*1 + *L*2)

diagonal matrix and eigenvectors matrix of **P**. |**P**| is then given by,

2 (*L*1 + *L*2) + *n*2

2*c*0 (−*L*1 + *L*2) *ny*

2 (−*L*1 + *L*2) *ny c*0

2 (*L*1 + *L*2 − 2*L*4) *n*2

*yL*4 *nxny*

*y*

where,

**3.5 Application**

⎛

*L*3 *nx*

0 *n*2 *x*

0 *nxny*

with, *V*0 = (*u*0, *v*0) and *c*0 the speed of sound.

2 exp �

The source term is made dimensionless with �

−*ln*(2)

(*x* − *xs*)

**3.5.1 Validation - convected monopole**

*xs* = *ys* = 0, and is defined as:

boundaries circumference.

the results in Bailly & Juvé (2000).

*Sp* = 1

0 *nx c*0

⎜⎜⎜⎜⎜⎝

The external state is then expressed as

$$\mathbf{U}^{\*-} = \mathbf{R}\mathbf{U}^{+} \tag{14}$$

where *R* is a transition matrix function of *n*'s components, it reads

$$\mathbf{R} = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 - 2n\_X^2 & -2n\_X n\_y & 0 \\ 0 & -2n\_X n\_y & 1 - 2n\_y^2 & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix} \tag{15}$$

Using this condition, the mass flux computed by the Riemann solver is zero and the non-permeability condition is satisfied.

#### **3.4.2 Absorbing boundary conditions**

Constructing absorbing (non-reflecting) boundary conditions for CAA is pretty delicate because of the high sensitivity of the accuracy to the small spurious wave reflexions at far field boundaries. Approaches based on the characteristics theory are not suited for CAA problems, other approaches, such as Perfectly Matched Layers (PML), presented by Hu (1996), and radial boundary condition of Tam & Webb (1993) are more indicated and widely discussed in the literature for finite differences schemes.

In this work we employ upwinding technique used by Bernacki et al. (2006) with DG to select only outgoing waves at the outer boundaries. Intuitively, it means that whole waves energy is dissipated at boundaries but unfortunately nothing prove that energy is actually dissipated and no spurious wave reflexions persist. To overcome this inconvenient, we join to the above procedure a grid stretching zone, see Nogueira et al. (2010). Grid stretching transfers the energy of the wave into increasingly higher wavenumber modes and the numerical scheme removes this high-frequency content. With this process most of the energy of the wave is dissipated before reaching the boundaries. For a high wavenumber the numerical method introduces more dissipation.

At the grid stretching zone, it is possible to use the MLS method as a filter in unstructured grids. The filtering process is developed by the application of a MLS reconstruction of the variables, i.e:

$$\mathbf{U}(\mathbf{x}) = \sum\_{j=1}^{n\_{x\_I}} \mathbf{U}(\mathbf{x}) N\_j(\mathbf{x}) \tag{16}$$

where, *U* is the reconstructed variable, *U*¯ is the filtered variable and *N* is the MLS shape function. This reconstruction is performed by using a kernel with shape parameters favoring dissipative behavior that the ones used to the approximation of the variables. The value of these parameters determines the range of frequencies to be filtered.

At the outer boundaries, we propose the following explicit numerical flux,

$$\mathbf{H}(\mathbf{U}^{\mathrm{n}}, \mathbf{U}^{\star \mathrm{n}}, \mathfrak{n}) = \frac{1}{2} (\mathbb{P}(\mathbf{U}^{\mathrm{n}}) \cdot \mathfrak{n} + |\mathbb{P}| \mathbf{U}^{\mathrm{n}-1}) \tag{17}$$

with,

*U*∗*n* is the fictitious state corresponding to the absorbing side ensuring **P***U*∗*n* = |**P**|*Un*−1, **P** is

the Jacobien matrix of system (11) and |**P**| = **V**−1|**D**|**V**, **D** and **V** are respectively, eigenvalues diagonal matrix and eigenvectors matrix of **P**. |**P**| is then given by,

$$\begin{aligned} \left| \mathbf{P} \right| = \begin{pmatrix} \mathbf{L}\_{3} & \frac{\mathbf{n}\_{y}}{2\mathbf{c}\_{0}}(-L\_{1} + L\_{2}) & \frac{\mathbf{n}\_{y}}{2\mathbf{c}\_{0}}(-L\_{1} + L\_{2}) & \frac{-1}{\mathbf{c}\_{0}^{2}}L\_{3} + \frac{1}{2\mathbf{c}\_{0}}(L\_{1} + L\_{2}) \\ 0 & \frac{\mathbf{n}\_{x}^{2}}{2}(L\_{1} + L\_{2}) + \mathbf{n}\_{y}^{2}L\_{4} & \frac{\mathbf{n}\_{x}\mathbf{n}\_{y}}{2}(L\_{1} + L\_{2} - 2L\_{4}) & \frac{\mathbf{n}\_{x}}{2\mathbf{c}\_{0}}(-L\_{1} + L\_{2}) \\ 0 & \frac{\mathbf{n}\_{x}\mathbf{n}\_{y}}{2}(L\_{1} + L\_{2} - 2L\_{4}) & \frac{\mathbf{n}\_{y}^{2}}{2}(L\_{1} + L\_{2}) + \mathbf{n}\_{x}^{2}L\_{4} & \frac{\mathbf{n}\_{y}}{2\mathbf{c}\_{0}}(-L\_{1} + L\_{2}) \\ 0 & \frac{\mathbf{n}\_{x}\mathbf{c}\_{0}}{2}(-L\_{1} + L\_{2}) & \frac{\mathbf{n}\_{y}\mathbf{c}\_{0}}{2}(-L\_{1} + L\_{2}) & \frac{1}{2}(L\_{1} + L\_{2}) \end{pmatrix} \end{aligned} \tag{18}$$

where,

18 Will-be-set-by-IN-TECH

10 00

0 −2*nxny* 1 − 2*n*2

Using this condition, the mass flux computed by the Riemann solver is zero and the

Constructing absorbing (non-reflecting) boundary conditions for CAA is pretty delicate because of the high sensitivity of the accuracy to the small spurious wave reflexions at far field boundaries. Approaches based on the characteristics theory are not suited for CAA problems, other approaches, such as Perfectly Matched Layers (PML), presented by Hu (1996), and radial boundary condition of Tam & Webb (1993) are more indicated and widely discussed in the

In this work we employ upwinding technique used by Bernacki et al. (2006) with DG to select only outgoing waves at the outer boundaries. Intuitively, it means that whole waves energy is dissipated at boundaries but unfortunately nothing prove that energy is actually dissipated and no spurious wave reflexions persist. To overcome this inconvenient, we join to the above procedure a grid stretching zone, see Nogueira et al. (2010). Grid stretching transfers the energy of the wave into increasingly higher wavenumber modes and the numerical scheme removes this high-frequency content. With this process most of the energy of the wave is dissipated before reaching the boundaries. For a high wavenumber the numerical method

At the grid stretching zone, it is possible to use the MLS method as a filter in unstructured grids. The filtering process is developed by the application of a MLS reconstruction of the

where, *U* is the reconstructed variable, *U*¯ is the filtered variable and *N* is the MLS shape function. This reconstruction is performed by using a kernel with shape parameters favoring dissipative behavior that the ones used to the approximation of the variables. The value of

2

*U*∗*n* is the fictitious state corresponding to the absorbing side ensuring **P***U*∗*n* = |**P**|*Un*−1, **P** is

*nxI* ∑ *j*=1

*U*¯ (*x*) =

these parameters determines the range of frequencies to be filtered.

At the outer boundaries, we propose the following explicit numerical flux,

*H*(*Un*,*U*∗*n*,*n*) = 1

00 01

*x* −2*nxny* 0

*y* 0

⎞

⎟⎟⎠

0 1 − 2*n*2

*U*∗− = *RU*+ (14)

*U*(*x*)*Nj*(*x*) (16)

(**F**(*Un*) · *n* + |**P**|*Un*−1) (17)

(15)

The external state is then expressed as

non-permeability condition is satisfied.

**3.4.2 Absorbing boundary conditions**

literature for finite differences schemes.

introduces more dissipation.

variables, i.e:

with,

where *R* is a transition matrix function of *n*'s components, it reads

*R* =

⎛

⎜⎜⎝

$$\begin{aligned} L\_1 &= |\mathbf{V\_0} \cdot \mathfrak{n} - c\_0| \\ L\_2 &= |\mathbf{V\_0} \cdot \mathfrak{n} + c\_0| \\ L\_3 &= L\_4 = |\mathbf{V\_0} \cdot \mathfrak{n}| \end{aligned} \tag{19}$$

with, *V*0 = (*u*0, *v*0) and *c*0 the speed of sound.

#### **3.5 Application**

#### **3.5.1 Validation - convected monopole**

This case reproduces the example of Bailly & Juvé (2000). The radiation of a monopole source is computed in a subsonic mean flow, with Mach number *Mx* = 0.5. The source is located at *xs* = *ys* = 0, and is defined as:

$$S\_p = \frac{1}{2} \exp\left(-\ln(2)\frac{\left(\mathbf{x} - \mathbf{x}\_s\right)^2 + \left(y - y\_s\right)^2}{2}\right) \sin\left(\omega t\right) \times \left[1, 0, 0, 1\right]^T \tag{20}$$

where the angular frequency is *ω* = 2*π*/30 and *t* is the time coordinate. The wave length is *λ* = 30 units, and the computational domain is a square with 200 units for each side. The source term is made dimensionless with � *ρ*0*c*0/Δ*x*, 0, 0, *ρ*0*c*3 0/Δ*x* �*T* . With the aim of testing the stability and the behavior of the proposed method for the boundary conditions, an unstructured grid absorbing layer has been added. The absorbing layer is placed from the boundary of the computational domain to *x* = ±300 and *y* = ±300. In figure 17 (top) it is shown the grid used for the resolution of this problem. To build this grid, 800 equally spaced nodes at the circumference of the computational domain are used and 120 nodes at the outer boundaries circumference.

In addition to the absorbing boundary condition given by equation (17), the shape filter parameters of the absorbing layer are *κx* = *κy* = 8, see Nogueira et al. (2010) paper.

A 5th order mass matrix based FV-MLS solver is used for this example.

Two acoustic waves propagate upstream and downstream of the source, and due to the effect of the mean flow, the apparent wavelength is modified and it is different upstream (*λ*1 = (1 − *Mx*) *λ*) and downstream (*λ*2 = (1 + *Mx*) *λ*) of the source.

In figure 18 pressure isocontours for different non-dimensional times *t* are shown. The pressure profile along axis *y* = 0 at time *t* = 270 is reproduced in figure 19, and also match the results in Bailly & Juvé (2000).

Fig. 19. Convected monopole - M=0.5, 200 × 200

Fig. 20. 3D and 2D centrifugal fan geometric

used in this problem is shown in figure (21).

are linked downstream by a vaned diffuser , see figure (20).

the diffuser blades at their interface.

The centrifugal fan noise is usually dominated by tones produced by the impeller blade passage. The resulted tonal noise corresponds to the blade passage frequency (BPF) and its higher harmonics. This is a consequence of the strong interaction between the impeller and

Aerodynamic and Aeroacoustic Study of a High Rotational Speed Centrifugal Fan 147

Shrouded impellers are usually used in high-rotational speed centrifugal fans. The impellers

A methodology based on a hybrid modeling of the aeroacoustic behavior of a high-rotational speed centrifugal fan is presented in this section. The main objective is to visualize the wave propagation into this machine and demonstrate, then, the power of the proposed methodology. Linearized Euler's equations are used to propagate noise radiated by the rotor/stator interaction. The fluctuating forces at the interaction zone are obtained by an aerodynamic study of the centrifugal fan presented by Khelladi et al. (2005a; 2008). In this section we calculate the acoustic wave propagation of a centrifugal fan with a 9-bladed rotor and a diffuser with 17 blades, as shown in Fig. 13. For the computations we use an unstructured grid, with at least 10 points per wavelength. A detail of the unstructured grid

**3.5.2 Acoustic waves propagation into a centrifugal fan**

Fig. 17. Convected monopole - M=0.5, 200 × 200 grid, *t* = 270 at *y* = 0

Fig. 18. Convected monopole - M=0.5, 200 × 200 grid. a) t=60, b) t=90, c) t=150, d) t=210, e) t=270

In order to check the stability of the boundary conditions, we let the computations to continue until 180 periods of the source. This time corresponds to the travel of the wave until the four outer boundaries. Comparing the pressure field with the one corresponding to *t* = 270 (9 source periods), it is observed that there is no change in the solution. The acoustic wave is completely dissipated by teh buffer zone when it lefts the computational domain, see figure (17) (top).

Fig. 19. Convected monopole - M=0.5, 200 × 200

20 Will-be-set-by-IN-TECH

a) b) c)

d) e)

Fig. 18. Convected monopole - M=0.5, 200 × 200 grid. a) t=60, b) t=90, c) t=150, d) t=210, e)

In order to check the stability of the boundary conditions, we let the computations to continue until 180 periods of the source. This time corresponds to the travel of the wave until the four outer boundaries. Comparing the pressure field with the one corresponding to *t* = 270 (9 source periods), it is observed that there is no change in the solution. The acoustic wave is completely dissipated by teh buffer zone when it lefts the computational domain, see figure

t=270

(17) (top).

Fig. 17. Convected monopole - M=0.5, 200 × 200 grid, *t* = 270 at *y* = 0

#### **3.5.2 Acoustic waves propagation into a centrifugal fan**

The centrifugal fan noise is usually dominated by tones produced by the impeller blade passage. The resulted tonal noise corresponds to the blade passage frequency (BPF) and its higher harmonics. This is a consequence of the strong interaction between the impeller and the diffuser blades at their interface.

Shrouded impellers are usually used in high-rotational speed centrifugal fans. The impellers are linked downstream by a vaned diffuser , see figure (20).

A methodology based on a hybrid modeling of the aeroacoustic behavior of a high-rotational speed centrifugal fan is presented in this section. The main objective is to visualize the wave propagation into this machine and demonstrate, then, the power of the proposed methodology. Linearized Euler's equations are used to propagate noise radiated by the rotor/stator interaction. The fluctuating forces at the interaction zone are obtained by an aerodynamic study of the centrifugal fan presented by Khelladi et al. (2005a; 2008). In this section we calculate the acoustic wave propagation of a centrifugal fan with a 9-bladed rotor and a diffuser with 17 blades, as shown in Fig. 13. For the computations we use an unstructured grid, with at least 10 points per wavelength. A detail of the unstructured grid used in this problem is shown in figure (21).

and interference. Thus, all these effects will be explicitly represented in the far field. Note that

Aerodynamic and Aeroacoustic Study of a High Rotational Speed Centrifugal Fan 149

they are not represented when other approach is used (FW&H, for example).

Fig. 22. Acoustic pressure history

Fig. 21. 2D centrifugal fan mesh

#### **Sources modeling**

If we refer to FW&H analogy, one can identify three acoustic sources of three different natures:

When the quadrupole source is included, substantially more computational resources are needed for volume integration. However, in many subsonic applications the contribution of the quadrupole source is small. Thus, we have neglected it in this calculation. Moreover, the monopolar source is also neglected at low Mach numbers and small surface thickness.

In our case, the interaction between the impeller and the diffuser blades is the main source of noise radiated by the centrifugal fan, see Khelladi et al. (2008). It is expressed by a pres- sure fluctuation on impeller and diffuser blades. It is, then, of a dipolar nature. This study takes into account only sources located at trailing edge of impeller blades and at the leading edge of blades of diffuser. The rotation of the impeller blades is modeled by rotating sources. Impeller blades are not taken into account in the propagation zone. Thus, we place 17 stationary bipolar acoustic sources located at the leading edge of the blades of diffuser and 9 additional rotating impeller sources located at the trailing edge of each impeller blade.

As for FW&H analogy, the source terms introduced in the LEE are constructed from the momentum equations, and defined by:

$$S\_{p\_i}(\mathbf{x}, y, t) = e^{-\frac{\ln(2)}{2} \left[ \left( \mathbf{x} - \mathbf{x}\_{\boldsymbol{\theta}\_i}(\mathbf{x}, y, t) \right)^2 + \left( y - y\_{\boldsymbol{\theta}\_i}(\mathbf{x}, y, t) \right)^2 \right]} \times p\_i(t) \times \left[ \mathbf{0}, \boldsymbol{n}\_{\mathbf{x}\_{\boldsymbol{\theta}\_i}} \boldsymbol{n}\_{\mathbf{y}\_{\boldsymbol{\theta}\_i}} \mathbf{0} \right]^T \tag{21}$$

the subscribe *i* corresponds to blades id, the position of the of the sources are defined by there coordinates (*xs*(*x*, *y*, *t*), *ys*(*x*, *y*, *t*)). For the impeller the sources moves following a circle path, the diffuser sources are static. *pi*(*t*) is the aerodynamic static pressure and (*nx*, *ny*) are the components of the unit radial vector at sources (*xs*, *ys*). The exponential term of equation 21 models the punctual nature of the considered sources.

Acoustic pressure history is presented in figure 22. At the beginning of the simulation we can observe clearly the position of sources. But soon we lose track of them because of reflections 22 Will-be-set-by-IN-TECH

If we refer to FW&H analogy, one can identify three acoustic sources of three different natures: • Monopole or thickness source: it is a surface distribution due to the volume displacement

• Dipole source or loading source: it is a surface distribution due to the interaction of the

When the quadrupole source is included, substantially more computational resources are needed for volume integration. However, in many subsonic applications the contribution of the quadrupole source is small. Thus, we have neglected it in this calculation. Moreover, the monopolar source is also neglected at low Mach numbers and small surface thickness. In our case, the interaction between the impeller and the diffuser blades is the main source of noise radiated by the centrifugal fan, see Khelladi et al. (2008). It is expressed by a pres- sure fluctuation on impeller and diffuser blades. It is, then, of a dipolar nature. This study takes into account only sources located at trailing edge of impeller blades and at the leading edge of blades of diffuser. The rotation of the impeller blades is modeled by rotating sources. Impeller blades are not taken into account in the propagation zone. Thus, we place 17 stationary bipolar acoustic sources located at the leading edge of the blades of diffuser and 9 additional rotating

As for FW&H analogy, the source terms introduced in the LEE are constructed from the

the subscribe *i* corresponds to blades id, the position of the of the sources are defined by there coordinates (*xs*(*x*, *y*, *t*), *ys*(*x*, *y*, *t*)). For the impeller the sources moves following a circle path, the diffuser sources are static. *pi*(*t*) is the aerodynamic static pressure and (*nx*, *ny*) are the components of the unit radial vector at sources (*xs*, *ys*). The exponential term of equation 21

Acoustic pressure history is presented in figure 22. At the beginning of the simulation we can observe clearly the position of sources. But soon we lose track of them because of reflections

(*x*,*y*,*t*)) 2

] × *pi*(*t*) ×

0, *nxi* , *nyi* , 0*T* (21)

• Quadrupole source: it is a volume distribution due to the flow outside the surfaces.

impeller sources located at the trailing edge of each impeller blade.

(*x*,*y*,*t*)) 2 +(*y*−*ysi*

− *ln*(2) 2 [(*x*−*xsi*

models the punctual nature of the considered sources.

Fig. 21. 2D centrifugal fan mesh

of fluid during the motion of surfaces.

flow with the moving bodies

momentum equations, and defined by:

*Spi*(*x*, *y*, *t*) = *e*

**Sources modeling**

and interference. Thus, all these effects will be explicitly represented in the far field. Note that they are not represented when other approach is used (FW&H, for example).

Fig. 22. Acoustic pressure history

Ishida, M. & Surana, T. (2005). Suppression of unstable flow at small flow rates in a centrifugal

Aerodynamic and Aeroacoustic Study of a High Rotational Speed Centrifugal Fan 151

Khelladi, S., Kouidri, S., Bakir, F. & Rey, R. (2005a). Flow study in the impeller-diffuser interface of a vaned centrifugal fan, *Journal of Fluid Engineering* 127: 495–502. Khelladi, S., Kouidri, S., Bakir, F. & Rey, R. (2005b). A numerical study on the aeroacoustics

Khelladi, S., Kouidri, S., Bakir, F. & Rey, R. (2008). Predicting tonal noise from a high rotational speed centrifugal fan, *Journal of Sound and Vibration* 313, Issues 1-2: 113–133. Khelladi, S., Nogueira, X., Bakir, F. & Colominas, I. (2011). Toward a higher order unsteady

Krain, H. (2005). Review of centrifugal compressor's application and developement, *Journal*

Lazzaretto, A. (2003). A criterion to define cross-flow fan design parameters, *Journal of Fluids*

Lazzaretto, A., Toffolo, A. & Martegani, A. D. (2003). A systematic experimental approach to

Lighthill, M. J. (1952). On sound generated aerodynamically, i. general theory, *Proceedings of*

Lighthill, M. J. (1954). On sound generated aerodynamically, ii. turbulence as a source of

Meakhail, T. (2005). A study of impeller-diffuser-volute interaction in a centrifugal fan, *Transaction of the ASME Journal of Turbomachinery, Vol.127, Jan., 84-90* . Menter, F. (1993). Zonal two equation *k* − *ω* turbulence models for aerodynamic flows, *AIAA*

Moon, Y.-J., Cho, Y. & Nam, H.-S. (2003). Computation of unsteady viscous flow and aeroacoustic noise of cross flow fans, *Computers & Fluids, Vol. 32, 995-1015* . Moon, Y. & Seo, H. J. (2006). Linearized perturbed compressible equations for low mach

Nogueira, X., Cueto-Felgueroso, L., Colominas, I., Khelladi, S., Navarrina, F. & Casteleiro,

Schleer, M. & Hong, S. (2004). Investigation of an inversely designed centrifugal compressor stage - part ii: Experimental investigation, *Journal of Turbomachinery* 126: 82–90. Seo, S., Kim, K. & Kang, S. (2003). Calculations of three-dimentional viscous flow in a

M. (2010). Resolution of computational aeroacoustics problem on unstructured grids with high-order finite volume scheme, *Journal of Computational and Applied*

multiblade centrifugal fan by modeling blade forces, *Journal of Power and Energy*

Lynam, E. & Webb., H. (1919). The emission of sound by airscrews., *R. & M., No. 624* . Maaloum, A., Kouidri, S. & Rey, R. (2004). Aeroacoustic performances evaluation of axial

cross-flow fan design, *Transaction of the ASME Journal of Fluids Engineering, Vol. 125,*

fans based on the unsteady pressure field on the blades surface, *Applied Acoustics*

URL: *http://www.sciencedirect.com/science/article/pii/S0045782511001344* Kim, K. & Seo, S. (2004). Shape optimisation of forward-curved-blade centrifugal fan with

navier-stokes analysis, *Journal of Fluid Engineering* 126: 735–742.

sound, *Proceedings of the Royal Society of London* A222: 1–32.

number aeroacoustics., *Comput. Phys.* 218(2): 702–719.

*Applied Mechanics and Engineering* 200(29-32): 2348 – 2362.

*of Turbomachinery, Vol. 127, Jan., 76-83* .

*Houston, Texas* FEDSM2005-77134.

*of Turbomachinery* 127: 25–34.

*Engineering, Vol.125, 680-683* .

*the Royal Society of London* A211: 564–587.

*July, 684, 693* .

65: 367–384.

*Paper* 93-2906.

217: 287–297.

*Mathematics* 234 (7): 2089–2097.

fan by controlling tip leakage flow and reverse flow, *Transaction of the ASME Journal*

of a vaned centrifugal fan, *ASME Heat Transfer/Fluids Engineering Summer Conference,*

finite volume solver based on reproducing kernel methods, *Computer Methods in*

#### **3.6 Conclusion**

In this section, we presented the aeroacoustic analysis of the centrifugal fan studied aerodynamically above. For this kind of machines, the effect of the near-field (casing and confinement) is very important which makes it very difficult to use acoustic analogy techniques. Linearized Euler equations seem to be a good alternative. In this section we presented the first developments and results concerning this approach with complex geometries such as the fan we analyzed. Nevertheless, the effect of rotating blades of the impeller was modeled by rotating sources, only the stator is taken into account. Taking into account the moving parts into the propagation domain is a our challenge for the near future. A technique based on sliding mesh method is under-development, the first results were very promising but some difficulties subsist to maintain the overall space order of accuracy due to the rotor/stator interface stencil reconstructions.

#### **4. References**

FLUENT (1998). Fluent, inc.

24 Will-be-set-by-IN-TECH

In this section, we presented the aeroacoustic analysis of the centrifugal fan studied aerodynamically above. For this kind of machines, the effect of the near-field (casing and confinement) is very important which makes it very difficult to use acoustic analogy techniques. Linearized Euler equations seem to be a good alternative. In this section we presented the first developments and results concerning this approach with complex geometries such as the fan we analyzed. Nevertheless, the effect of rotating blades of the impeller was modeled by rotating sources, only the stator is taken into account. Taking into account the moving parts into the propagation domain is a our challenge for the near future. A technique based on sliding mesh method is under-development, the first results were very promising but some difficulties subsist to maintain the overall space order of accuracy due to

Bailly, C. & Juvé, D. (2000). Numerical solution of acoustic propagation problems using

Bardina, J., Huang, P. & Coakley, T. (1997). Turbulence modeling, validation, testing and

Bernacki, M., Lanteri, S. & Piperno, S. (2006). Time-domain parallel simulation of

Bryan, G. H. (1920). The acoustics of moving sources with application to airscrews., *R. & M.*

Colonius, T. & Lele, S. (2004). Computational aeroacoustics: progress on nonlinear problems

Colonius, T., Lele, S. & Moin, P. (1993). Boundary conditions for direct computation of

Curle, N. (1955). The influence of solid boundaries upon aerodynamic sound, *Proceedings of*

Eum, H.-J. & Kang, S.-H. (2002). Numerical study on tip clearance effect on performance of

Farassat, F. & Myers, M. K. (1988). Extension of kirchhoff's formula to radiation from moving

Ffowcs Williams, J. & Hawkings, D. (1969). Sound generation by turbulence and

Gutin, L. (1948). On the sound field of a rotating propeller, *NACA TM1195, (Traduction de "Uber* ¨

Hu, F. (1996). On absorbing boundary conditions for linearized Euler equations by a perfectly

matched layer, *Journal of Computational Physics* 129: 201–219.

a centrifugal compressor, *Proceedings of ASME FEDSM'02, Montreal, Quebec, Canada,*

surfaces in arbitrary motion, *Philosophical Transactions for the Royal Society of London*

*das Schallfeld einer rotierenden Luftschraube", Physikalische Zeitschrift der Sowjetunion,*

discontinuous Galerkin methods, *J. Computational Acoustics* 14(1): 57–82. Bogey, C., Bailly, C. & JuvÃl', D. (2002). Computation of flow noise using source terms in

heterogeneous wave propagation on unstructured grids using explicit, non-diffusive,

linearized Euler equations, *AIAA Journal* 38(1): 22–29.

development, *NASA Technical Memorandum 110446* .

of sound generation., *Prog Aerosp Sci* 40: 345–416.

*the Royal Society of London* A231: 505–514.

surfaces, *J. Sound and Vib.* 123: 451–560.

*Band 9, Heft 1, 1936).* pp. 57–71.

linearized eulerŠs equations., *AIAA Journal* 40(2): 235–243.

aerodynamic sound generation, *AIAA Journal* 31(9): 1574–1582.

**3.6 Conclusion**

**4. References**

the rotor/stator interface stencil reconstructions.

*No. 684, British A.R.C.* .

*July 14-18* .

A264: 321–342. FLUENT (1998). Fluent, inc.

navier-stokes analysis, *Journal of Fluid Engineering* 126: 735–742.

**Section 4**

**Nonstationary Aerodynamics**

**and CFD in Compressible Flow**

## **Section 4**

**Nonstationary Aerodynamics and CFD in Compressible Flow**

26 Will-be-set-by-IN-TECH

152 Applied Aerodynamics

Tam, C. & Webb, J. (1993). Dispersion-relation-preserving finite difference schemes for

Toffolo, A. (2004). On cross-flow fan theoretical performance and efficiency curves: An

Zangeneh, M. & Schleer, M. (2004). Investigation of an inversely designed centrifugal

energy loss analysis on experimental data, *Transaction of the ASME Journal of Fluids*

compressor stage - part i: Design and numerical verification, *Journal of Turbomachinery*

computational aeroacoustics, *J. Comput. Phys.* 107: 262–281.

*Engineering, Vol. 126, Sept., 743-751* .

126: 73–81.

**8**

*2CONICET*

*Argentina*

**An Approximate Riemann Solver**

A fundamental subject leading to numerical simulations of Euler equations by the Finite Volume (FV) method, is the calculation of numerical fluxes at cells interfaces. The degree of accuracy of the FV numerical scheme, its ability to capture discontinuities and the correct prediction of the velocity of propagating waves, are all flow properties strongly dependent

In many numerical schemes, the fluxes between cells are computed using truncated series expansions which based only on numerical considerations. These considerations to be strictly valid, must account for some degree of continuity in the functions and in their derivatives, but clearly, these continuity conditions are not satisfied when discontinuous solutions as shock waves or contact surfaces are present in the flow. This type of flow problems nevertheless, were solved with relative success until 1959. In that year Godunov published his work "A finite difference method for the computation of discontinuous solutions of the equations of Fluid Dynamics" (Godunov, 1959) in which an alternative approach for solving the system of Euler equations is presented. This new approach, in striking difference to previous ones, is basically supported by physical considerations and

The excellent results obtained with the Godunov technique, prompted several researches to develop new FV numerical schemes for two and three dimensional applications, achieving second order accuracy and total variation diminishing (TVD) properties (Toro, 2009; LeVeque, 2004; Yee, 1989). These new schemes were built around the use of Riemann solvers, making them generally very accurate but computationally expensive. Such high computational cost is attributed to the iterative technique required to solve the system of five nonlinear algebraic equations needed to find, in all cells, an exact solution of the Riemann problem. Alternative schemes computationally less demanding, although less accurate and less robust, were then

built based on approximate solutions of the Riemann problem (Toro, 2009).

**1. Introduction**

on the evaluation of numerical fluxes.

the essential part of it is the so called Riemann solver.

**for Euler Equations**

Oscar Falcinelli1, Sergio Elaskar1,2,

*1Aeronautical Department, FCEFyN, Nacional University of Cordoba*

*Nacional University of La Plata*

José Tamagno1 and Jorge Colman Lerner3

*3Aeronautical Department, Engineering Faculty,*

## **An Approximate Riemann Solver for Euler Equations**

Oscar Falcinelli1, Sergio Elaskar1,2, José Tamagno1 and Jorge Colman Lerner3 *1Aeronautical Department, FCEFyN, Nacional University of Cordoba 2CONICET 3Aeronautical Department, Engineering Faculty, Nacional University of La Plata Argentina*

#### **1. Introduction**

A fundamental subject leading to numerical simulations of Euler equations by the Finite Volume (FV) method, is the calculation of numerical fluxes at cells interfaces. The degree of accuracy of the FV numerical scheme, its ability to capture discontinuities and the correct prediction of the velocity of propagating waves, are all flow properties strongly dependent on the evaluation of numerical fluxes.

In many numerical schemes, the fluxes between cells are computed using truncated series expansions which based only on numerical considerations. These considerations to be strictly valid, must account for some degree of continuity in the functions and in their derivatives, but clearly, these continuity conditions are not satisfied when discontinuous solutions as shock waves or contact surfaces are present in the flow. This type of flow problems nevertheless, were solved with relative success until 1959. In that year Godunov published his work "A finite difference method for the computation of discontinuous solutions of the equations of Fluid Dynamics" (Godunov, 1959) in which an alternative approach for solving the system of Euler equations is presented. This new approach, in striking difference to previous ones, is basically supported by physical considerations and the essential part of it is the so called Riemann solver.

The excellent results obtained with the Godunov technique, prompted several researches to develop new FV numerical schemes for two and three dimensional applications, achieving second order accuracy and total variation diminishing (TVD) properties (Toro, 2009; LeVeque, 2004; Yee, 1989). These new schemes were built around the use of Riemann solvers, making them generally very accurate but computationally expensive. Such high computational cost is attributed to the iterative technique required to solve the system of five nonlinear algebraic equations needed to find, in all cells, an exact solution of the Riemann problem. Alternative schemes computationally less demanding, although less accurate and less robust, were then built based on approximate solutions of the Riemann problem (Toro, 2009).

An Approximate Riemann Solver for Euler Equations 157

where *Δu = uR - uL* . Once Eq.(3a) is solved for *p\** the solution for *u\** is obtained from Eq.(3b) and the remaining unknowns are found by means of standard gas dynamic relations. The functions *fL* and *fR* represent relations across the left non-linear wave and across the right

> \* \* \* \*

*A f p U p p if pp p B*

*KK K K*

*K K K K a p f p U if pp p*

1 2 \* \* \*

 1 1 *BK K p*

*K*

, *K*

2 , 1 1 *K*

>

 2

1 *K*

*K A*

where *γ* is the ratio of specific heats, *K* may be *L* or *R* depending on if the flow changes are evaluated across the nonlinear left or right waves, and *a*K is the speed of sound in the left or right states. The Eq.(3a) is an implicit algebraic nonlinear equation on the only unknown *p\**, and it must be solved using an iterative scheme. Once the pressure in the star region has been obtained, it is possible to calculate by means of explicit expressions the velocity *u\** and

In Eq.(3a), neither *uL* nor *uR* are explicitly written, but only its difference. Therefore the

In this article, the dimensional analysis is used to reduce the number of independent variables necessaries for describing the gas-dynamics Riemann problem. From this point of view, it is possible to consider as reference variables the density and pressure from one side

The Eq.(6) represents a relationship between pressures, velocities and densities. By means of the dimensional analysis, the densities can be written in non-dimensional form using a

*R pR*)1/2. Then, the solution of a particular Riemann problem is determined by the

1 23 *L L R R R*

1 1 \* , , 2 2 *u u u f pU f pU LR R R L L* (3b)

(4a)

(4b)

(5)

,,, , *LL RR* (6)

*ref R r ef R p p* (7)

*R*, the pressures using as reference *pR*, and the velocities respect of

(8)

\* \*

non-linear wave respectively, and are given by

the density at each side of the contact discontinuity.

pressure in the star zone becomes only function of five variables:

*R*

*u p p p*

*p*\* 1 *fu p p*

of the Riemann problem, for example, the right side.

(

reference value such as

following three parameters:

In this article, a new non-iterative way of solving the full Riemann problem, applicable to three-dimensional and time dependent Euler equations, is presented. This non-iterative solution requires at the beginning of the computation and only once, the generation in situ of tabulated exact Riemann solutions from which the information needed can be retrieved using multiple-linear interpolation. To reduce the time to access during successive calculations, the original five independent variables of the exact Riemann problem are, by way of dimensional analysis, reduced to only three thus allowing the build-up of an easiest to handle data-based matrix with three degrees of freedom.

#### **2. Description of the proposed Riemann Solver**

Compressible and inviscid flow problems are governed by the Euler equations. In one dimension (1D), these equations can be written as:

$$\left. \mathcal{U}\_t + F\left(\mathcal{U}\right)\_x = 0 \right. \tag{1}$$

$$\mathcal{U} = \begin{bmatrix} \rho \\ \rho u \\ E \end{bmatrix}; \quad \mathcal{F}(\mathcal{U}) = \begin{bmatrix} \rho u \\ \rho u^2 + p \\ (E+p)u \end{bmatrix} \tag{2}$$

*U* is the vector of conservative variables, *F* the vector of convective flows, *t* and *x* are the temporal and spatial coordinates respectively, is the density, *u* the velocity in the x direction, *p* the pressure and *E* the total energy per unit volume. The subscript means differentiation with respect to time and space.

The Riemann problem solves de Euler equations in a 1D domain where initial conditions are given by two different constant states separated by a discontinuity. The solution of this 1D Riemann problem depends directly on the ratio *x/t*, and it will consist of three types of waves: two nonlinear, shock or expansion fan, and one linearly degenerate, the contact discontinuity. These waves are separating four constant states where the conservative vector *U* acquires from the left to the right the following values, *UL, UL\*, UR\** and *UR*. The subscripts ''*L*'' and ''*R*'' indicate left and right, respectively, and the symbol "*\**" identify points located in the state between the nonlinear waves (star region).

To obtain the flow variations produced by the waves, the Riemann invariants technique for expansion and contact waves and of the Rankine-Hugoniot relationship for shock waves, must be implemented. In the Riemann problem, this generally leads to an algebraic system of nine equations with nine unknowns. The unknowns are the velocities of the three waves plus six variables necessary to characterize the states *UL\*, UR\**. However, an analysis of the eigenstructure of the Euler equations allows to establish that both pressure *p\** and particle velocity *u\** between the left and right waves are constant, while the density take the two constant values *ρ\*L* and *ρ\*R* (Toro, 2009). Based on these considerations and after some algebraic calculations, the two equations listed below are obtained and used for solving the Riemann problem.

$$f\_{\perp} \left( p\_\*, \mathcal{U}\_{\perp} \right) + f\_{\mathbb{R}} \left( p\_\*, \mathcal{U}\_{\mathbb{R}} \right) + \Delta \mathfrak{u} = \mathbf{0} \tag{3a}$$

156 Applied Aerodynamics

In this article, a new non-iterative way of solving the full Riemann problem, applicable to three-dimensional and time dependent Euler equations, is presented. This non-iterative solution requires at the beginning of the computation and only once, the generation in situ of tabulated exact Riemann solutions from which the information needed can be retrieved using multiple-linear interpolation. To reduce the time to access during successive calculations, the original five independent variables of the exact Riemann problem are, by way of dimensional analysis, reduced to only three thus allowing the build-up of an easiest

Compressible and inviscid flow problems are governed by the Euler equations. In one

2 ;

*E Ep u*

*U u FU u p*

*U* is the vector of conservative variables, *F* the vector of convective flows, *t* and *x* are the

direction, *p* the pressure and *E* the total energy per unit volume. The subscript means

The Riemann problem solves de Euler equations in a 1D domain where initial conditions are given by two different constant states separated by a discontinuity. The solution of this 1D Riemann problem depends directly on the ratio *x/t*, and it will consist of three types of waves: two nonlinear, shock or expansion fan, and one linearly degenerate, the contact discontinuity. These waves are separating four constant states where the conservative vector *U* acquires from the left to the right the following values, *UL, UL\*, UR\** and *UR*. The subscripts ''*L*'' and ''*R*'' indicate left and right, respectively, and the symbol "*\**" identify points located in

To obtain the flow variations produced by the waves, the Riemann invariants technique for expansion and contact waves and of the Rankine-Hugoniot relationship for shock waves, must be implemented. In the Riemann problem, this generally leads to an algebraic system of nine equations with nine unknowns. The unknowns are the velocities of the three waves plus six variables necessary to characterize the states *UL\*, UR\**. However, an analysis of the eigenstructure of the Euler equations allows to establish that both pressure *p\** and particle velocity *u\** between the left and right waves are constant, while the density take the two constant values *ρ\*L* and *ρ\*R* (Toro, 2009). Based on these considerations and after some algebraic calculations, the two equations listed below are obtained and used for solving the

  *u*

0 *t x U FU* (1)

(2)

(3a)

is the density, *u* the velocity in the x

to handle data-based matrix with three degrees of freedom.

**2. Description of the proposed Riemann Solver**

dimension (1D), these equations can be written as:

temporal and spatial coordinates respectively,

the state between the nonlinear waves (star region).

*f pU f pU u L LR R* \* \* , ,0

Riemann problem.

differentiation with respect to time and space.

$$u^\* = \frac{1}{2}(u\_\perp + u\_k) + \frac{1}{2} \left[ f\_\hbar \left( p\_\cdot, lL\_k \right) - f\_\hbar \left( p\_\cdot, lL\_\perp \right) \right] \tag{3b}$$

where *Δu = uR - uL* . Once Eq.(3a) is solved for *p\** the solution for *u\** is obtained from Eq.(3b) and the remaining unknowns are found by means of standard gas dynamic relations. The functions *fL* and *fR* represent relations across the left non-linear wave and across the right non-linear wave respectively, and are given by

$$f\_{\chi} \left( p\_\*, \mathcal{U}\_{\chi} \right) = \left( p\_\* - p\_{\chi} \right) \cdot \sqrt{\frac{A\_{\chi}}{p\_\* + B\_{\chi}}} \qquad \text{if} \qquad p\_\* \ge p\_{\chi} \tag{4a}$$

$$f\_{\mathbb{K}}\left(p\_{\ast}, \mathcal{U}\_{\mathbb{K}}\right) = \frac{2 \cdot a\_{\mathbb{K}}}{\gamma - 1} \cdot \left[ \left(\frac{p\_{\ast}}{p\_{\mathbb{K}}}\right)^{\frac{\gamma - 1}{2\gamma}} - 1 \right] \qquad \text{if} \qquad p\_{\ast} < p\_{\mathbb{K}} \tag{4b}$$

$$A\_{\mathbb{K}} = \frac{2}{\left(\gamma + 1\right) \cdot \rho\_{\mathbb{K}}} \quad B\_{\mathbb{K}} = \frac{\left(\gamma - 1\right)}{\left(\gamma + 1\right)} \cdot p\_{\mathbb{K}} \tag{5}$$

where *γ* is the ratio of specific heats, *K* may be *L* or *R* depending on if the flow changes are evaluated across the nonlinear left or right waves, and *a*K is the speed of sound in the left or right states. The Eq.(3a) is an implicit algebraic nonlinear equation on the only unknown *p\**, and it must be solved using an iterative scheme. Once the pressure in the star region has been obtained, it is possible to calculate by means of explicit expressions the velocity *u\** and the density at each side of the contact discontinuity.

In Eq.(3a), neither *uL* nor *uR* are explicitly written, but only its difference. Therefore the pressure in the star zone becomes only function of five variables:

$$p\_\* = f\_1\left(\Delta t, \rho\_{\perp'} p\_{\perp'}, \rho\_{\mathbb{k}'} p\_{\mathbb{k}'}\right) \tag{6}$$

In this article, the dimensional analysis is used to reduce the number of independent variables necessaries for describing the gas-dynamics Riemann problem. From this point of view, it is possible to consider as reference variables the density and pressure from one side of the Riemann problem, for example, the right side.

$$
\rho\_{ref} = \rho\_{\mathbb{R}} \qquad \qquad \qquad p\_{ref} = p\_{\mathbb{R}} \tag{7}
$$

The Eq.(6) represents a relationship between pressures, velocities and densities. By means of the dimensional analysis, the densities can be written in non-dimensional form using a reference value such as *R*, the pressures using as reference *pR*, and the velocities respect of (*R pR*)1/2. Then, the solution of a particular Riemann problem is determined by the following three parameters:

$$
\pi\_1 = \frac{\Delta u}{\sqrt{\frac{p\_k}{\rho\_k}}} \qquad \qquad \qquad \pi\_2 = \frac{p\_L}{p\_R} \qquad \qquad \qquad \pi\_3 = \frac{\rho\_L}{\rho\_R} \tag{8}
$$

An Approximate Riemann Solver for Euler Equations 159

 *R L A B R L A B A B p p u*

Unnecessary storage in symmetrical cases can be avoided if it is adopted as the selection

density, but those of the higher pressure side. Thus, the reference state is the higher pressure

*ref ref ref*

*u u p a*

*2* value will always be less than or equal to one.

a change of the Mach number because *a*ref is not the sound velocity neither that of the left

Finally, it is clear that no interpolation is necessary when the solution of the Riemann problem possesses left and right expansion waves, because in this case the Eq.(3) can be

To analyze the accuracy and computational efficiency of the proposed Riemann solver, comparisons with others solvers available in the literature (Toro, 2009), are made. Three were selected: one that, iteratively searches for the exact solution of Eq.(3), and others two, that they try to solve the Riemann problem with approximate schemes. To build-up the comparison, Riemann problems were generated randomly with values of the parameters

The following sub-sections explain how each one of the solvers selected for comparison

This solver searches by means of an iterative process the solution of Eq.(3a). This equation

*f pU f pU u R p L LR R* , ,

Usually to iteratively solve Eq.(13), the Newton-Raphson method is implemented. This method requires the calculations of the function as well as of its derivative, which should

where *R(p)* is the residual to cancel; *fL* and *fR*, are calculated according to Eqs.(4 and 5).

2

3 0.05 5.05 (12)

(13)

0.05 1

change between the left and right states in the Riemann problem. However

*ref*) in Eq.(7), not the left or right pressure and

1 variable it is convenient to re-define it as:

(11)

1 would be like a Mach number

1 is not strictly

1,

*u* (10)

Two Riemann problems will be symmetric if the following conditions are satisfied:

criterion of the reference variables (*pref ,*

1

To give a physical sense to the non-dimensional

where *a*ref is the sound velocity in the reference state and

**3. Comparison with other Riemann Solvers**

3 ranging between the limits set below:

1 10.05 4.95

initial state and the

state or that of the right.

solved analytically.

2 and

works.

**3.1 Iterative Riemann Solver**

for any value of *p* shall be written as:

Different cases of Riemann problems with identical values of *π*1, *π*2, *π*3 must have similar behavior, and the calculated relation *p\*/pR* must have the same value for all cases. To ascertain this behavior for the non-dimensional approach proposed, several numerical examples are analyzed. Seven Riemann problems whose left and right conditions are indicated in the following table, and all of them satisfying the parameters 1 = -0.78262, 2 = 50, 3 = 10, were tested:

Table 1. Test cases.

\_\_\_\_\_\_\_

The relation *p\*/pR* is obtained solving Eq.(3), and for the seven test cases considered it is found *p\*/pR* = 13.312, which proves that the proposed non-dimensional analysis works properly. Then, it is possible to write the Eq.(3) only as function of 1, 2, 3 and *p\*/pR*.

Since in solving the Eq.(3) there are involved only three independent variables, a data-base matrix with three degrees of freedom containing the values of *p\*/pR* for N values of *1*, M of *2*, and Q of *3* is, in situ generated. Then, to find the solution of a particular Riemann problem, it is only necessary to calculate the corresponding values of 1, 2, 3 and to interpolate for *p\*/pR* in the N*x*M*x*Q matrix (from now, simply called A-matrix). Finally the pressure in the star zone can be calculated as:

$$p\_\* = \overline{\left(\frac{p\_\*}{p\_\mathbb{R}}\right)} \cdot p\_\mathbb{R} \tag{9}$$

where \* *R p p* is the interpolated value from the A-matrix. After calculating the pressure in

the star region, the rest of the variables can be explicitly calculated using the same equations of the exact solver (Toro, 2009).

The previously described procedure involves the use of an interpolated value of pressure to calculate the density and velocity changes across each wave; however it is not the only possible procedure. Others alternatives are for instance, to develop arrangements with dimensionless density or dimensionless velocity variations across each wave.

To increase the approximate solution accuracy, it is desirable that the variation range of the parameters 1, 2, 3 be as small as possible. One way is to reduce the range of 2 or 3 avoiding unneeded storage of data in symmetrical Riemann problems.

Different cases of Riemann problems with identical values of *π*1, *π*2, *π*3 must have similar behavior, and the calculated relation *p\*/pR* must have the same value for all cases. To ascertain this behavior for the non-dimensional approach proposed, several numerical examples are analyzed. Seven Riemann problems whose left and right conditions are

1.2250E+00 1.0000E+02 1.0000E+05 1.2250E-01 0.0000E+00 2.0000E+03 4.9071E+01 8.4770E+02 8.7460E+06 4.9071E+00 6.9994E+02 1.7492E+05 6.7304E+00 6.5231E+02 1.0554E+07 6.7304E-01 2.1402E+02 2.1108E+05 4.1503E+00 7.8027E+02 1.1631E+07 4.1503E-01 1.9437E+02 2.3261E+05 9.4504E+00 6.4262E+02 1.5976E+07 9.4504E-01 1.8756E+02 3.1952E+05 3.0289E+01 2.9038E+02 2.9757E+06 3.0289E+00 1.8067E+02 5.9514E+04 3.6284E+01 3.0129E+02 6.5687E+05 3.6284E+00 2.5420E+02 1.3137E+04

The relation *p\*/pR* is obtained solving Eq.(3), and for the seven test cases considered it is found *p\*/pR* = 13.312, which proves that the proposed non-dimensional analysis works

Since in solving the Eq.(3) there are involved only three independent variables, a data-base

interpolate for *p\*/pR* in the N*x*M*x*Q matrix (from now, simply called A-matrix). Finally the

\_\_\_\_\_\_\_ \* \* *R R p p p p*

the star region, the rest of the variables can be explicitly calculated using the same equations

The previously described procedure involves the use of an interpolated value of pressure to calculate the density and velocity changes across each wave; however it is not the only possible procedure. Others alternatives are for instance, to develop arrangements with

To increase the approximate solution accuracy, it is desirable that the variation range of the

3 be as small as possible. One way is to reduce the range of

dimensionless density or dimensionless velocity variations across each wave.

avoiding unneeded storage of data in symmetrical Riemann problems.

*3* is, in situ generated. Then, to find the solution of a particular Riemann

is the interpolated value from the A-matrix. After calculating the pressure in

matrix with three degrees of freedom containing the values of *p\*/pR* for N values of

problem, it is only necessary to calculate the corresponding values of

*R uR R p*

1, 2,

3 and *p\*/pR*.

1, 2,

(9)

*1*, M of

2 or 3

3 and to

1 = -0.78262,

2 =

indicated in the following table, and all of them satisfying the parameters

properly. Then, it is possible to write the Eq.(3) only as function of

50,

where

parameters

\_\_\_\_\_\_\_ \* *R p p*

Table 1. Test cases.

*2*, and Q of

of the exact solver (Toro, 2009).

1, 2,

pressure in the star zone can be calculated as:

3 = 10, were tested:

*L uL L p*

Two Riemann problems will be symmetric if the following conditions are satisfied:

$$\left(\left(p\_{\:k}\right)\_{A} = \left(p\_{\!\!\perp}\right)\_{\!\!B} \qquad\qquad \qquad \left(\rho\_{\!\!k}\right)\_{A} = \left(\rho\_{\!\!\perp}\right)\_{\!\!B} \qquad\qquad \qquad \Delta u\_{A} = \Delta u\_{\!\!\perp} \tag{10}$$

Unnecessary storage in symmetrical cases can be avoided if it is adopted as the selection criterion of the reference variables (*pref , ref*) in Eq.(7), not the left or right pressure and density, but those of the higher pressure side. Thus, the reference state is the higher pressure initial state and the *2* value will always be less than or equal to one.

To give a physical sense to the non-dimensional 1 variable it is convenient to re-define it as:

$$
\pi\_1 = \frac{\Delta u}{\sqrt{\nu \cdot \frac{p\_{nf}}{\rho\_{nf}}}} = \frac{\Delta u}{a\_{nf}}\tag{11}
$$

where *a*ref is the sound velocity in the reference state and 1 would be like a Mach number change between the left and right states in the Riemann problem. However 1 is not strictly a change of the Mach number because *a*ref is not the sound velocity neither that of the left state or that of the right.

Finally, it is clear that no interpolation is necessary when the solution of the Riemann problem possesses left and right expansion waves, because in this case the Eq.(3) can be solved analytically.

#### **3. Comparison with other Riemann Solvers**

To analyze the accuracy and computational efficiency of the proposed Riemann solver, comparisons with others solvers available in the literature (Toro, 2009), are made. Three were selected: one that, iteratively searches for the exact solution of Eq.(3), and others two, that they try to solve the Riemann problem with approximate schemes. To build-up the comparison, Riemann problems were generated randomly with values of the parameters 1, 2 and 3 ranging between the limits set below:

$$-10.05 \le \pi\_1 \le 4.95 \quad 0.05 \le \pi\_2 \le 1 \quad 0.05 \le \pi\_3 \le 5.05 \tag{12}$$

The following sub-sections explain how each one of the solvers selected for comparison works.

#### **3.1 Iterative Riemann Solver**

This solver searches by means of an iterative process the solution of Eq.(3a). This equation for any value of *p* shall be written as:

$$f\_{\perp} \left( p\_{\prime} \mathcal{U}\_{\perp} \right) + f\_{\mathbb{R}} \left( p\_{\prime} \mathcal{U}\_{\mathbb{R}} \right) + \Delta u = R \left( p \right) \tag{13}$$

where *R(p)* is the residual to cancel; *fL* and *fR*, are calculated according to Eqs.(4 and 5).

Usually to iteratively solve Eq.(13), the Newton-Raphson method is implemented. This method requires the calculations of the function as well as of its derivative, which should

An Approximate Riemann Solver for Euler Equations 161

the function given by Eq.(4), and thus slightly shortening the time needed to reach the

A particular solution of the gas-dynamics Riemann problem is given when both non-linear waves are rarefaction waves. In this case the pressure in the star region can be obtained

> 1 1 2 2 \* \*

2 2 1 1 0

*a a p p u*

*L R*

*p p*

*L R*

*aa u p a a*

error incurred is at most of second order if instead of Eq.(4a), the Eq.(4b) is utilized.

1 2

>

1 1 2 2

>

*L R*

When the pressure exceeds the value of *p*min, and the Eq.(4a) instead of Eq.(4b) to calculate the *fK* functions is used, there are no discontinuities in the residual function and in its derivatives. This particular behavior of the Eq.(4a) and of its first derivative implies that the

However, if the Eq.(4b) is used to calculate the function *fK* the ability to obtain analytically the value *p*\* , is lost. The TRRS method, suggests to calculate the pressure in the star region always using the Eq.(17), no matter what kind of Riemann problem is studied. The error of this method will be null for the Riemann problems with two rarefaction waves ( *p*\* min *p* )

The TSRS is the opposite case of the TRRS. In the TSRS the solution of the Riemann problem is obtained considering that both non-linear waves are shocks. Then the functions *fK* are

> 2 2 1 1

1 1 1 1

*L R*

*L R*

, ,0

0

(18)

*L R*

*p p*

2 1

(17)

correct solution.

analytically.

where

**3.2 Two-Rarefaction Riemann Solver (TRRS)**

\*

and will increase as *p*\* becoming higher and moves away from *p*min.

given by the Eq.(4b) and the Eq.(3a) can be written as shown bellow:

*L R*

\* \*

*LL RR*

*p p gpU p p gpU u*

*p p p p u pp pp*

\* \*\* \*

\* \*

**3.3 Two-Shock Riemann solver (TSRS)**

1 1

*L R*

increase the computational cost of the method, but this increase is not significant because the evaluation of the derivative demands simple computations once the function has been evaluated. Since the derivative of the residual function, Eq.(13), is calculated by means of the functions *fL* and *fR* without considering *UR* and *UL* that are constant, it can then be written as:

$$\frac{d\left[\int f\_L\left(p, \mathcal{U}\_L\right)\right]}{dp} + \frac{d\left[\int f\_R\left(p, \mathcal{U}\_R\right)\right]}{dp} = \frac{d\left[\int R\left(p\right)\right]}{dp} \tag{14}$$

The derivatives ' , *L L L d f pU f dp* y ' , *R R R d f pU f dp* are calculated as:

$$f\_{\chi} = \begin{cases} \sqrt{\frac{A\_{\chi}}{B\_{\chi} + p}} \cdot \left(1 - \frac{p - p\_{\chi}}{2 \cdot (B\_{\chi} + p)}\right) & \text{if} \quad p \ge p\_{\chi} \\\\ \frac{1}{\rho\_{\chi} \cdot a\_{\chi}} \cdot \left(\frac{p}{p\_{\chi}}\right)^{\frac{-(\gamma + 1)}{2\gamma}} & \text{if} \quad p < p\_{\chi} \end{cases} \tag{15}$$

From the computational point of view it is noted that the most expensive steps involved in the numerical evaluation of Eq.(15) are the powers with fractional exponents. However, Eqs.(4a and 4b) and Eq.(15) show that these computational steps have already been made when the residual *R(p)* is computed, and as a result of this its derivative calculation is relatively fast.

The Newton-Raphson algorithm applied to Eq.(13) can be written as:

$$p\_{i+1} = p\_i - \frac{R\_{(p\_i)}}{R\_{(p\_i)}} \tag{16}$$

where *pi* and *pi*+1 are the pressure for the iteration *i* and *i*+1 respectively, *R* and *R'* are the residual function and its derivative.

In solving iteratively Eq.(13) difficulties may appears because, as Eq.(16) suggests, the pressure can becomes negative. To avoid this problem, the residual function is evaluated (considering *p = pmin* and *p = pmax*, where *pmin = min[pL,pR]* and *pmax = max[pL,pR]*. If both residuals are positive, the pressure which cancels the residual is less than *pL* and *pR*, and the Riemann problem has an explicit solution consisting of two rarefaction waves. If the residual corresponding to the maximum pressure is greater than zero and the corresponding to minimum pressure is less than zero, the Riemann problem have as solutions a shock wave and a rarefaction fan, and it is adopted as a first iteration the value of the minimum pressure (which will undoubtedly be lower than the pressure to cancel the residual). When both residuals are negative, the Riemann problem has two shock waves and the sought pressure will be greater than both, so as first iteration the maximum pressure value is adopted.

Using the above described procedure, it is possible to know beforehand the kind of solution expected for the Riemann problem at each iterations, which allows an effective selection of

increase the computational cost of the method, but this increase is not significant because the evaluation of the derivative demands simple computations once the function has been evaluated. Since the derivative of the residual function, Eq.(13), is calculated by means of the functions *fL* and *fR* without considering *UR* and *UL* that are constant, it can then be written as:

, , *LL RR d f pU d f pU d R p*

*d f pU f dp*

2

*K K*

*K K*

1 if

From the computational point of view it is noted that the most expensive steps involved in the numerical evaluation of Eq.(15) are the powers with fractional exponents. However, Eqs.(4a and 4b) and Eq.(15) show that these computational steps have already been made when the residual *R(p)* is computed, and as a result of this its derivative calculation is

1 ´

In solving iteratively Eq.(13) difficulties may appears because, as Eq.(16) suggests, the pressure can becomes negative. To avoid this problem, the residual function is evaluated (considering *p = pmin* and *p = pmax*, where *pmin = min[pL,pR]* and *pmax = max[pL,pR]*. If both residuals are positive, the pressure which cancels the residual is less than *pL* and *pR*, and the Riemann problem has an explicit solution consisting of two rarefaction waves. If the residual corresponding to the maximum pressure is greater than zero and the corresponding to minimum pressure is less than zero, the Riemann problem have as solutions a shock wave and a rarefaction fan, and it is adopted as a first iteration the value of the minimum pressure (which will undoubtedly be lower than the pressure to cancel the residual). When both residuals are negative, the Riemann problem has two shock waves and the sought pressure will be greater than both, so as first

Using the above described procedure, it is possible to know beforehand the kind of solution expected for the Riemann problem at each iterations, which allows an effective selection of

*i i*

where *pi* and *pi*+1 are the pressure for the iteration *i* and *i*+1 respectively, *R* and *R'*

KK K

The Newton-Raphson algorithm applied to Eq.(13) can be written as:

*a p*

*L d f pU f dp* y ' , *R R*

'

*K*

*f*

The derivatives ' , *L L*

residual function and its derivative.

iteration the maximum pressure value is adopted.

relatively fast.

*dp dp dp* (14)

are calculated as:

*K*

*p p R* (16)

(15)

are the

*K*

*R*

1 if

*A p p p p Bp Bp*

1 2

>

*R*

*p*

*i i p*

*p p p*

the function given by Eq.(4), and thus slightly shortening the time needed to reach the correct solution.