Dong-Ho Rhee
email: rhee@kari.re.kr

Jeong-Seek Kang

Sanga Lee

Young-Seok Kang

Bong Jun Cha

Experimental Study of Film Cooling Performance on Turbine Guide Vane with Optimized Hole Arrangements in a Transonic Annular Sector Cascade

INTRODUCTION

Film cooling is one of the most effective ways to cool down the hot components in gas turbine engines. Since early 1970's, the extensive works on the film cooling have been conducted and various concepts of the film cooling have been developed and used in the gas turbine engines, such as crater hole configuration and the holes with a trenched slot [START_REF] Fric | Method for improving the cooling effectiveness of a gaseous coolant stream which flows through a substrate, and related articles of manufacture[END_REF][START_REF] Dorrington | Film Effectiveness Performance for Coolant Holes Embedded in Various Shallow Trench and Crater Depressions[END_REF][START_REF] Kalghatgi | Improved Film Cooling Effectiveness with a Round Film Cooling Hole Embedded in a Contoured Crater[END_REF][START_REF] Bunker | Film Cooling Effectiveness Due To Discrete Holes within A Transverse Surface Slot[END_REF]. Recently, to overcome the limitation of the conventional design schemes, the optimization method has been introduced and developed as a possible way to improve the effectiveness and reduce the aerodynamic losses. While most of the researchers have dealt with the hole shape optimization [START_REF] Lee | Optimization of a Cylindrical Film Cooling Hole using Surrogate Modeling[END_REF][START_REF] Lee | Shape optimization of a fanshaped hole to enhance film-cooling effectiveness[END_REF][START_REF] Nita | Film Cooling Hole Shape Optimization Using Proper Orthogonal Decomposition[END_REF], only a few researchers have worked on the hole arrangement optimization. Johnson et al. [START_REF] Johnson | Design Optimization Methods for Improving HPT Vane Pressure Side Cooling Properties Using Genetic Algorithms and Efficient CFD[END_REF][START_REF] Johnson | Genetic Algorithm Optimization of a High-Pressure Turbine Vane Pressure Side Film Cooling Array[END_REF] studied the optimization of the hole patterns as well as the hole shapes on the high pressure turbine vane by CFD, but their study has a limitation in that the optimized configuration was determined among the previously selected designs and the film hole geometry is not exactly modeled due to the computation resource matter.

In this study, the film cooling performance on the pressure side surface of the vane is evaluated experimentally for the optimized hole arrangement proposed by Lee et al. [START_REF] Lee | Optimization of the Array of Film Cooling Holes on a High Pressure Turbine Nozzle[END_REF][START_REF] Lee | Optimization of Film Cooling Hole Array on Pressure Side of HPT Nozzle for Enhanced Cooling Performance[END_REF] and compared with the baseline configuration.

EXPERIMENTAL APPARATUS & PROCEDURES 1.1 EXPERIMENTAL APPARATUS

Figure 1 shows the schematic view of the transonic turbine cascade test facility in Korea Aerospace Research Institute. The mainstream is supplied by 500 HP compressor into the annular turbine cascade test section through the 300kW heater and a venturi flow meter. To simulate the coolant flow, 50 HP compressor supplies the air through the manifold and

Abstract

The present study investigated the film cooling performance on the pressure side surface of the nozzle guide vane with the optimized film cooling hole arrangement. An annular sector transonic turbine cascade test facility in Korea Aerospace Research Institute was used. The mainstream and the secondary air for cooling are supplied by 500hp and 50hp compressors, respectively. To measure the film cooling effectiveness on the surface, a transient measurement method was used with the FLIR infrared camera system and the mainstream temperature was maintained 15 to 20C above the secondary flow temperature using 300kW heater. Five nozzle guide vanes were installed in the test section. The configuration with the cooling hole arrangement proposed by the optimization procedure was tested and compared with the baseline configuration. For the cooling hole arrangement optimization, nine design variables are taken into account for a single objective function of film cooling effectiveness on the pressure side surface. Design of Experiment, Kriging surrogate model and Genetic Algorithm are used for the optimization. The tests were conducted at the mainstream exit Reynolds number based on the chord of 2.2x10 6 and the total coolant mass flow rate was set to be 10% of the inlet main flow. The flow periodicity in the cascade passage was ensured by the surface static pressure measurements. The results showed that considerable improvement of the film cooling effectiveness is observed in the whole region with the optimized configuration, which gives the possibility of coolant mass flow reduction. The present study used the annular sector turbine cascade with five vanes for the test model. With the annular turbine cascade, one can simulate the radial pressure gradient inside the passage as well as can study the effect of three dimensional shape of the guide vane. For that reason, some of the recent studies adopted the annular sectional cascade for turbine heat transfer investigation [START_REF] Luque | A Novel Technique for Assessing Turbine Cooling System Performance[END_REF][START_REF] Saha | Aerodynamic Investigations of a High Pressure Turbine Vane with Leading Edge Contouring at Endwall in a Transonic Annular Sector Cascade[END_REF].

BASELINE NOZZLE GUIDE VANE GEOMETRY

The external profile of the test vane is derived from the previous study on the aerodynamic design [START_REF] Kang | Aerodynamic Optimization of Axial Turbine Tip Cavity with Approximation Model[END_REF]. To meet the similarity condition of Reynolds number and to obtain sufficient resolution of the image for the detailed measurements, the test vanes are three-times scaled up. The vane has a chord length of 123.2 mm and an exit angle of 74. The vane model has a contoured endwall on the shroud side so that the span at leading edge is 30% longer than that at trailing edge. The pitch to chord ratio is 0.76 based on the midspan geometry. The details are presented in Table 1.

The cooling design features are presented in Fig. 2. The vane has two inner cooling plenums, plenum #1 and #2, having the internal cooling features such as pedestals, rib turbulators and an insert for jet impingement cooling. For the external surface cooling, seven rows of film cooling holes are located on the surface: three rows of circular film cooling holes at the leading edge region, a single row of fan-shaped holes on the suction side surface and three rows of fanshaped holes on the pressure side surface. The inlet diameter of the film cooling hole is 1.5 mm and the fan-shaped holes expand laterally with an angle of 10 at one-third of hole length from the inner surface of the vane wall. The trailing edge has a pressure side cutback with four cooling slots. In this study, the focus was on the film cooling of the pressure side surface as mentioned above. The geometric features of the pressure side film cooling holes are listed in Table 2.

Three rows of the cooling holes are named as PS1, PS2 and PS3, respectively as shown in Fig. 2. The first row, PS1, has an injection angle of 34 with respect to the surface while PS2 and PS3 have 30. The chordwise location of each row for the baseline configuration is described in Table 2.

OPTIMIZATION OF FILM COOLING HOLE ARRANGEMENT AND NUMERICAL SIMULATION

As mentioned above, two configurations of film cooling hole arrangements are considered in the present study to examine the cooling performance: the one is the baseline and the other one is the optimized case, of which the objective function is to maximize the film cooling effectiveness. The optimization of the hole arrangements was conducted numerically under the actual operating engine conditions. For example, the average total temperature and pressure at mainstream inlet is 1,673K and 30.7bar, respectively while the temperature of the coolant is 837K [START_REF] Lee | Optimization of the Array of Film Cooling Holes on a High Pressure Turbine Nozzle[END_REF][START_REF] Lee | Optimization of Film Cooling Hole Array on Pressure Side of HPT Nozzle for Enhanced Cooling Performance[END_REF]. The numerical model for the optimization was simplified to have the film cooling holes only on the pressure side surface as shown in Fig. 3 because this study only focused on the hole arrangement on the pressure side surface. Also, this can reduce computing resources and time significantly. To simulate the coolant mass flow rate distributions of the fully cooled nozzle vane, the inlet boundary conditions of two coolant plenums were carefully determined to have the same mass flow rate through the cooling holes on the pressure side surface.

The optimization procedure is shown in Fig. 4. For the optimization, Design of Experiment (DOE), Kriging surrogate model and Genetic Algorithm (GA) with Estimated Improvement (EI) are used. First, the design variables, the design space and the objective function(s) are determined for the cooling geometry. In this study, to define and optimize the hole arrangements on the pressure side surface, nine design variables such as chordwise and spanwise locations, arrangement curvature, gradient and spacing are controlled within a confined design space (Fig. 3(a)) while the number of holes, the hole inlet diameter and the injection angle are fixed. Also the position of PS1 is fixed to give a reference position to PS2 and PS3. The objective function is to minimize the adiabatic wall temperature on the film-cooled pressure side surface while the coolant mass flow rate remains to be equal to that for the baseline. A set of sampling points, that is, various combinations of the design variables for the cooling design, is established by Design of Experiment (DOE) and the numerical simulations are conducted for the sampling points to evaluate the film cooling effectiveness on the pressure side surface. The number of the sampling points was approximately 150 in this study.

For the numerical analysis, compressible steady Reynolds averaged Navier-Stokes (RANS) equations in ANSYS CFX 14.5 were used, where the Element-based Finite Volume Method (FVM) and Algebraic Multi-grid (AMG) coupled solver were used. For turbulence modeling, the SST model was chosen. The SST model is known to yield better prediction for near-wall turbulent flow and Ayoubi et al. [START_REF] Ayoubi | Aerothermal shape optimization for a double row of discrete film cooling holes on the suction surface of a turbine vane[END_REF] concluded that SST model shows good performance in the turbine heat transfer and film cooling hole problem. The number of grid for the numerical simulation is approximately 8 million. The average y + at the wall ranges approximately from 5 to 30, and the maximum y+ at the wall is kept under 100 as depicted in Fig. 5.

Once the Kriging model is constructed from the effectiveness results of the initial sampling points, Genetic Algorithm (GA) is performed to find the optimum values of design variables from the Kriging model. Then, additional numerical simulations are conducted for the geometry of the optimum design variables to ensure reliability of the Kriging model and to refine the Kriging model. Details of the optimization procedure are well described in the references [START_REF] Lee | Optimization of the Array of Film Cooling Holes on a High Pressure Turbine Nozzle[END_REF][START_REF] Lee | Optimization of Film Cooling Hole Array on Pressure Side of HPT Nozzle for Enhanced Cooling Performance[END_REF].

Figure 6 presents the baseline and the optimized hole arrangements on the pressure side surface. The chordwise locations of the rows of holes are listed in Table 3 as well as in Fig. 6. Comparing the baseline and the optimized configurations reveals that PS2 moves upstream and radially outward to make a staggered array with PS1. PS3 also moves far upstream. The interesting feature is that the hole spacing of PS2 and PS3 near the hub is slightly smaller for the optimized case than that for the baseline while the hole spacing near the shroud is larger. In addition, PS3 does not have a straight array of holes along the spanwise direction like the baseline configuration but a C-shaped curved pattern of the hole position.

EXPERIMENTAL CONDITIONS

To measure the film cooling effectiveness on the surface, a transient method was used. Table 4 describes the detailed test conditions in this study.

At first the air supplied from the compressor is heated by a heater and discharged through the bypass line as depicted in Fig. 1. Then the main valves are switched to open within 0.4 seconds after the temperature reaches the steady state, which can assure the step change in the mainstream temperature (Fig. 7). The mainstream temperature is controlled to be 15 to 20C higher than the coolant temperature, which results in the density ratio(c/m) of 1.06. The exit Reynolds number of the mainstream based on the chord is 2.2x10 6 for the experimental test. The corresponding average exit Mach number is 0.8.

The total coolant mass flow rate supplied into two plenums is 10% of the inlet mainstream in this study so that each cooling plenum has 5% of the coolant. Based on the previous CFD results [START_REF] Seo | Conjugate Heat Transfer Analysis of a Single Stage High-Pressure Turbine[END_REF], the coolant mass flow rate through the PS holes is estimated to be 26.9% of the total coolant mass flow: the coolant mass flow rate through PS1, PS2 and PS3 is 0.81% 0.71% and 1.17% of the mainstream, respectively. The corresponding effective blowing ratio at the hole exit ranges from 1.33 to 2.51 based on the inlet mass flux condition of the mainstream. When the pitchwise averaged mass flux at each row position is used for the mainstream, the effective blowing ratio ranges from 0.52 to 1.17 as presented in Table 5. The test models shown in Fig. 8 are manufactured by 3D printer, CMET EQ-1, which uses a stereo-lithography method. The solid properties of the test model, such as density, specific heat and thermal conductivity were carefully measured for a transient test data reduction. Two thermocouples were embedded on the surface from the inside of the cooling channel cavity and are used as the reference temperature for the in-situ calibration. In addition, the models with the pressure taps were manufactured for the static pressure measurements.

TEST PROCEDURE & DATA REDUCTION

For the measurements of the adiabatic film temperature (Tf) on the vane surface, 1-D semi-infinite conduction assumption is applied to the vane surface. The governing equation is defined as Eqn. [START_REF] Fric | Method for improving the cooling effectiveness of a gaseous coolant stream which flows through a substrate, and related articles of manufacture[END_REF].

𝑇 𝑤 -𝑇 𝑖 𝑇 -𝑇 𝑖 = 1 -exp ( ℎ 2  𝑘 2 ) 𝑟 𝑐 ( ℎ √ 𝑘 ) (1)
The solid properties (, k) in Eqn. (1) were obtained from the separate property measurement test. The initial temperature, Ti and the surface temperature, Tw on the pressure side surface are measured during the test using an IR camera at a sampling rate of 12Hz.

To calibrate the measured value from the IR camera, an insitu calibration was conducted in the test section for various mainstream temperature conditions. As mentioned above, the data from two thermocouples were used as the reference temperature.

The mainstream temperature was controlled by 300kW heater and the calibration data were taken when the steady state was assured. Figure 9 presents the calibration curve of IR camera data showing a quite linear pattern. Then, two unknown variables, hf and Tf can be obtained from the data at two different time frames in a single experiment. The details of the data reduction procedure follow the reference given by Ekkad et al. [START_REF] Ekkad | A Transient Infrared Thermography Method for Simultaneous Film Cooling Effectiveness and Heat Transfer Coefficient Measurements From a Single Test[END_REF]. In this study, the nominal test duration is within 60 seconds and the data of two time frames, for example 5 and 30 seconds, are taken from the surface temperature measurement results. The film cooling effectiveness is then calculated from Eqn. [START_REF] Dorrington | Film Effectiveness Performance for Coolant Holes Embedded in Various Shallow Trench and Crater Depressions[END_REF]. The recovery temperature is used as a reference temperature of the mainstream in Eqn. [START_REF] Dorrington | Film Effectiveness Performance for Coolant Holes Embedded in Various Shallow Trench and Crater Depressions[END_REF]. The definition of the recovery temperature is presented in Eqn. [START_REF] Kalghatgi | Improved Film Cooling Effectiveness with a Round Film Cooling Hole Embedded in a Contoured Crater[END_REF]. In Eqn. (3), recovery factor, r, is set to be Pr 1/3 and the incoming Mach number (Min) of the mainstream is used for this study. The mainstream static temperature (Ts) is measured upstream of the center nozzle and the coolant temperature (Tc), the temperature inside the coolant plenums, is measured by four thermocouples inside the plenums during the test.

 = 𝑇 -𝑇 𝑟 𝑇 𝑐 -𝑇 𝑟

(2) Estimated uncertainties in surface temperature measurement using IR camera and the solid properties are within 1.0C and 1%, respectively. The test duration has 1.6% uncertainty based on the time frame used in this study. Using the method by Kline and McClintock [START_REF] Kline | Describing Uncertainties in Single Sample Experiments[END_REF], the uncertainties for film cooling effectiveness are estimated 7.8% at =0.6. However, the uncertainty becomes much higher with lower film cooling effectiveness. For example, the uncertainty is 48% at =0.1 while 6% at =0.8.

𝑇 = 𝑇 1 + 𝑟  -1 𝑀 𝑖𝑛 (3)

RESULTS AND DISCUSSION 2.1 STATIC PRESSURE MEASUREMENTS

To check the periodicity of the flow through the passages, the static pressure on the surface was measured and compared with the numerical simulation results [START_REF] Lee | Optimization Framework using Surrogate Model for Aerodynamically Improved 3D Turbine Blade Design[END_REF] for the nozzle guide vane without cooling holes. Figure 10 shows the distributions of the pressure coefficient on the midspan for the center and adjacent vanes. It is noted that the axial locations of the rows of the cooling holes for the baseline and the optimized cases are indicated in Fig. 10 for the reference.

The pressure coefficient is defined in Eqn. [START_REF] Bunker | Film Cooling Effectiveness Due To Discrete Holes within A Transverse Surface Slot[END_REF].

Cp = (Pt -Ps,surface) / (Pt -Ps,exit) (4)
The results show that the static pressure distributions are quite well coincident with the CFD results all over the surface including the throat area (x/Cx0.65 on the suction side) and the adjacent nozzles also have quite similar distributions to the center nozzle vane.

FILM COOLING EFFECTIVENESS MEASUREMENTS

Figure 11 presents the contour plots of film cooling effectiveness on the pressure side surface for the baseline and the optimized cases. It is noted that some of the effectiveness data on the trailing edge slot surface are scattered because (Tw-Ti) in Eqn. ( 1) is maintained almost constant during the test, which makes Tf solution scattered.

The baseline has an in-line array for PS1 and PS2, which results in very low effectiveness values between the holes as shown in Fig. 11(b). Although the interaction of the coolant ejected from PS1 and PS2 increases the effectiveness in the downstream region of PS2, still low effectiveness values are found between the holes. High film cooling effectiveness values are observed in the downstream region of PS3 when compared with the values in the upstream region. There are several reasons for that. A favorable pressure gradient formed in the downstream region suppresses the coolant to the surface and consequently improve the film coverage on the surface. Low local blowing ratio at PS3 with the upstream coolant from PS1 and PS2 also improves the film cooling effectiveness on that region. From this result, it can be concluded that improvement of the local effectiveness in the region between PS2 and PS3 is essential to increase the overall film cooling performance while maintaining high effectiveness values in the downstream region of PS3. The optimized design in this study seems to follow this guideline when the characteristics of the optimized design are examined. As shown in Fig. 11(c), the optimized configuration has more uniform and higher film cooling effectiveness all over the surface as expected. In the upstream region, the coverage of film is significantly better because PS2 moves upstream and shifted upward in the spanwise direction so that PS1 and PS2 have a staggered array. Consequently, the coolant effectively protects the inter-hole region and high effectiveness values are maintained until the flow reaches PS3. When the regions downstream of PS3 are compared, the levels of the effectiveness are similar. However the region just downstream of PS3, especially the near-hub region, has much larger area of high effectiveness with the optimized configuration. This is due to the coalescence of the coolant resulting from smaller hole spacing with the upstream coolant protection effect in that region.

Figure 12 shows the spanwise distributions of film cooling effectiveness downstream of rows of holes. Downstream of PS2 (Fig. 12(a)), the highest values of effectiveness are around 0.7 for both cases. However, the regions between holes have quite different values. This is due to the relative position of the holes at PS1 and PS2 as mentioned earlier.

For the optimized configuration, PS2 moves upward in spanwise direction and the hole spacing is rearranged, that is, has a relatively small spacing near the hub. The superposition of the coolant from PS1 and PS2 makes higher and more uniform effectiveness values in the downstream region of PS2. Similar patterns are observed in the downstream region of PS3 as shown in Fig. 12(b). In the region of z/H0.5, the effectiveness values between the holes are much higher for the optimized case while the peak values are almost the same. In the region upstream of the trailing edge slot (Fig. 12(c)), the mean values become very close to each other configuration showing effectiveness of around 0.5 in the overall region, but much better uniformity is obtained with the optimized case.

Figure 13 shows the streamwise distributions of the spanwise averaged film cooling effectiveness for the baseline and the optimized configurations. The effectiveness around PS1 shows the similar values for both the cases. However the downstream regions of PS2 and PS3 are found to have much higher effectiveness values with the optimized configuration as expected from the contour and local plots. This reveals that the optimization process modifies the cooling design to strengthen favorable interaction of coolant ejecting from the upstream and downstream rows. As a result, for the optimized case, relatively high effectiveness values are found even far downstream region of PS3.

The film cooling effectiveness contours downstream of the first row obtained by the CFD are shown in Fig. 14. It is noted that the contour range is different from that of the experimental results to show the coolant behaviors on the surface more clearly. The CFD cases shows the elevated level of effectiveness values when compared with the experimental results. There are two things to be considered on these results. First, the boundary conditions for the CFD were based on the real engine condition as mentioned earlier, which means the difference in density ratio of the coolant to the mainstream. Second, due to the limitation of computing resources, the grid quality, especially near the surface, is not as good as that of fully dedicated three dimensional numerical simulations having y + value of around 1 for a single calculation. These may cause over-predicted film cooling effectiveness on the pressure side surface when comparing with the experimental results. However the overall patterns are quite similar to those of experiments qualitatively. This means that the optimization procedure and its methodology are quite reasonable to use for improvement of the cooling performance. 15 presents the average film cooling effectiveness for two cases. In the experimental cases, the optimized case shows 16.5% higher effectiveness value than the baseline case. The CFD results also consistent trend with the experimental ones showing 22.3% increase in the average effectiveness. However, as mentioned above, the CFD overpredicts 14.7%p to 20.5%p above the experimental values.

SUMMARY

The present study investigated the film cooling performance on the pressure side surface of the nozzle guide vane with the optimized film cooling hole arrangements experimentally and compared the results with the numerical simulation results from the optimization.

The results showed that significant improvement of film cooling effectiveness is obtained with the optimized hole arrangement. The optimization process changes the hole arrangement to strengthen the superposition effect of the coolant ejecting from the rows of the holes. From those results, it is clearly shown that the optimization procedure and methodology is well established and can be used for not only improvement of film cooling effectiveness but also the coolant mass flow reduction achievement.

The numerical simulations were conducted in the optimization process and their results were compared with the experimental results. The comparison reveals that the experimental results are qualitatively coincident with the CFD results, which ensures the reliability of the optimization procedure and its output.

Future works will focus on the performance of various optimization configurations so that the optimization procedure and methodology can be validated and revised.

Figure 1 .Figure 2 .

 12 Figure 1. Schematic diagram of test facility

3 .

 3 (a) film cooling holes with design space (b) numerical domain with boundary conditions Figure Simplified model for the optimization process and the numerical analysis

Figure 4 .Figure 5 .

 45 Figure 4. Optimization procedure[START_REF] Lee | Optimization of the Array of Film Cooling Holes on a High Pressure Turbine Nozzle[END_REF][START_REF] Lee | Optimization of Film Cooling Hole Array on Pressure Side of HPT Nozzle for Enhanced Cooling Performance[END_REF]

Figure 7 . 8 .Figure 9 .

 789 Figure 7. Typical example of mainstream temperature variation during the transient test

Figure 10 . 11 .

 1011 Figure 10. Pressure coefficients at mid-span on the nozzle surface

12 .Figure 13 .

 1213 Figure 13. Streamwise distributions of spanwise averaged film cooling effectiveness for different hole arrangements

Figure 14 .Figure 15 .

 1415 Figure 15. Average film cooling effectiveness comparison

Figure

  Figure15presents the average film cooling effectiveness for two cases. In the experimental cases, the optimized case shows 16.5% higher effectiveness value than the baseline case. The CFD results also consistent trend with the experimental ones showing 22.3% increase in the average effectiveness. However, as mentioned above, the CFD overpredicts 14.7%p to 20.5%p above the experimental values.

Table 1 .

 1 Nozzle guide vane geometry

	Parameters	Values
	Number of vanes	5
	Scale-up factor	3
	Chord, C	123.2 mm
	Inlet / Exit angle []	0 / 74
	Aspect ratio (l/C)	Inlet: 0.53, Exit: 0.40
	Spacing (P/C) at midspan	0.76 at inlet
		• Two coolant plenums with
	Cooling design	Internal cooling features • Seven rows of film cooling
	(depicted in Fig. 2)	holes • Four slots at pressure side
		cutback

Table 2 .

 2 Pressure side surface film cooling geometry based on the baseline configuration

	Parameters	Values
	Number of rows of holes	3 (PS1 and PS2 at plenum
	on the pressure side surface	#1, PS3 at plenum #2)
	Hole inlet diameter	1.5 mm
	Injection angle	34 for PS1 30 for PS2 and PS3
	Laterally Expansion angle	10
	Hole spacing (p/d)	3.6
	Axial positions of PS1, PS2 and PS3	PS1: x/Cx=15.4% PS2: x/Cx=27.6% PS3: x/Cx=66.5%

Table 3 .

 3 Chordwise locations of film cooling rows

	Parameters	Baseline	Optimized
	PS1	x/Cx=15.4%	x/Cx=15.4%
	PS2	x/Cx=27.6%	x/Cx=18.3%
	PS3	x/Cx=66.5%	x/Cx=39.5%*
	Injection	34 for PS1
	angle	30 for PS2 and PS3

*nominal value for the curved array of PS3

Table 4 .

 4 Test conditions

	Parameters	Values
	Exit Reynolds number, ReC,ex	2.2 x 10 6
	Mainstream temperature	316K
	Inlet turbulence intensity at midspan	2.6%
	Exit Mach number	0.8
	Total coolant mass flow rate, mc	10% of mainstream
	mc at PS1	0.81%
	mc at PS2	0.71%
	mc at PS3	1.17%
	Coolant temperature	298K
	Density ratio (2/)	1.06

Table 5 .

 5 Effective blowing ratio at hole exit for the baseline configuration

	Position	Meffective based on the inlet main flow	Meffective based on the pitchwise averaged main flow
	PS1	1.38	1.17
	PS2	1.33	0.52
	PS3	2.51	0.54

ACKNOWLEDGMENTS

This project is funded by Ministry of Industry, Trade & Energy through the Aerospace Component Technology Development Program (grant No. KA000157)


Tezhuan Du
email: dutezhuan@imech.ac.cn

Chenguang Huang

Yiwei Wang

A Numerical Model for Evolution of Internal Structure of Cloud Cavitation

Keywords:

Interaction between bubbles in cloud cavitation play an important role in cavitating flows. A numerical model was developed to study the internal structure and collapse of cloud cavitation. The model consists of : (1) an evolution model of bubble number density taking bubble breakup effects into account; (2) a modified cavitation model based on dynamics of bubble cluster; and (3) the multi-phase Reynolds Averaged Navier-Stokes equations (RANS) for background flow. The evolution model of bubble number density is governed by a transport type equation and the expression of source term is derived from a bubble breakup model. The condensation rate is constructed through dimension analysis and direct simulation of collapse of bubble cluster. The proposed model was tested by flows over a projectile. The results show that the evolution of internal structure of the cloud cavitation is closely related to the development of re-entrant jet. Moreover, the performance of proposed model and the effects of bubble number density on cavitating flows are discussed.

INTRODUCTION

Cavitation is the key issue in the high-speed underwater propulsion, which can result in structure failure and acoustic. Cloud cavitation, as one of the most common type of cavitation for the marine vessels and projectile, is consists of a large amount of small bubbles. The pressure pulse generated by the bubble collapse in the cloud cavitation is usually considered to be the major cause of the structure failure and the noise radiation.

Numerical simulations are commonly used to study cavitating flows in many applications, such as hydrofoils, projectiles, and turbo machines, etc. Most of the researches were focused on the homogeneous flow modelling which is based on the Navier-Stokes equations of mixture phase. Kubota [START_REF] Kubota | Unsteady structure measurement of cloud cavitation on a foil section[END_REF],Merkle [START_REF] Merkle | Computational modeling of the dynamics of sheet cavitation[END_REF], Kunz [START_REF] Kunz | A preconditioned Navier-Stokes method for two-phase flows with application to cavitation prediction[END_REF] and Singhal [START_REF] Singhal | Mathematical Basis and Validation of the Full Cavitation Model[END_REF] proposed a number of cavitation models based on the phase-change method respectively. However, the limitation of the homogeneous cavitation modeling is apparent, since it is only capable of providing macro-solutions of vapor volume fraction while the micro mechanism of cavitation evolution remains unknown. The evolution of the cloud cavitation includes various micro-physical processes such as the bubble collapse and the development of the re-entrant jet, and is highly related to the distribution of bubbles in terms of the size and number. For example, bubble clusters with different distributions may lead to different collapse pressures due to the interaction of bubbles while the vapor volume fractions can be identical. Therefore, the numerical model taking account of the internal structure of cloud cavitation will be essential for the understanding of the bubble evolution, and is expected to improve the prediction of the bubble collapse pressure.

Similar considerations of modeling small particle distributions can be found in many chemical engineering simulations. A statistical model of the dispersed phase is usually adopted in dispersed flows to describe the distribution of particle size [START_REF] Martínez-Bazán | On the breakup of an air bubble injected into a fully developed turbulent flow. Part 1: Breakup frequency[END_REF]- [START_REF] Sosinovich | Modle of gas bubble breakup in turbulent liquid flow[END_REF], in which the evolution of the distribution function is governed by a Boltzmann-type equation. Unlike the chemical engineering applications, in which most of the models were based on the steady state condition and fully developed homogeneous turbulence, the cavitating flows with the evolution of cloud cavitation are highly unsteady and far more complex in terms of the flow physics.

In the present paper, the evolution model of bubble number density is proposed to simulate the process of the bubble breakup and transportation. A condensation rate based on bubble cluster is deduced by dimensional analysis and direct simulation. A numerical strategy is established to solve the flow with cloud cavitation combining our model with the homogeneous cavitation model, and was applied to the simulation of the cavitating flow around a projectile. The solutions can predict the evolution of both the outline and the internal structure of the cloud cavitation, which are in good agreement with the experimental observations.

Evolution model of the bubble number density

The original homogeneous cavitation models can provide the vapor volume fraction only, but the information of the internal structure of the cloud cavitation such as the bubble number density and the averaged bubble size remains unknown. In the present work, we introduce a new transport equation of the bubble number density to the homogeneous cavitation model, so that the internal structure of cavitation can be solved. The transport equation of the bubble number density can be constructed as follows

() ( ) n un S m m n t      n un S n
un S and coalescence. In the present work, only the effect of bubble breakup is taken into account.

Cavitation cloud is normally a bubble cluster which consists of large amounts of small bubbles with various sizes. Bubble-bubble interaction plays an essential role in the collapse of the bubble cluster, and will become stronger while the bubble cluster gets more compact. Given a homogeneous bubble cluster, the vapor volume fraction  can be related to the bubble number density n and the averaged bubble size r as:

4 3 3 nr   (2)
Assuming the vapor volume fraction being constant in the process of bubble breakup, the source term of Eq. ( 1) can be written as follows

The averaged bubble size is given by: (, ) 0

r f t r rdr    (4)
where (, ) f tr is the probability density function of bubble size at time t .

Bubble breakup model

Bubbles break up in turbulent flows when the velocity of turbulent fluctuations exceeds a threshold value. In other word, the bubble breaks up in turbulent flows when the bubble size exceeds a critical value, which is related to the turbulent fluctuations and the surface tension. Refer to Sosinovich [START_REF] Levich | Physicochemical Hydrodynamics (English translation by Scripta Tchnica)[END_REF], the evolution equation of probability density function of bubble size can be written as:

(, ) (, ) (, ) ( ,) (, ) 11 () ( ) 0 r ftr ftr ftR r r dr r R dR tr R r           (5)
where () r  represents the breakup time of a bubble of radius r , and the breakup function (, ) rR  represents the probability of a bubble of size R breaking into bubbles of size r . Neglecting the convection effect, the effects of bubble breakup on the time rate of probability density function can be divided into two parts: The first term on right hand side of Eq. ( 5) represents the decrease in probability function due to the bubbles of size r breaking into smaller ones, the second term indicates the increase in probability function due to bubbles of radii 1 rr  break into ones of size r .

It is known that a bubble will break into two bubbles near the critical Weber number [START_REF] Sosinovich | Model of fractionation and coalescence of gas bubbles in a turbulent liquid flow[END_REF], based on which the breakup function can be given as follows:

(, ) ( ) ( ) ( ) 3 2 R rR Rr Ra r cr       (6)
where () x  is the step function, and a cr is the critical size, which is determined by the following equation [START_REF] Sosinovich | Model of fractionation and coalescence of gas bubbles in a turbulent liquid flow[END_REF]:

2 13 13 ( , ) (1 )
( )( ) 0 a cr a P t r dr C cr f l l a      (7)
where

0.5 C f  is drag coefficient of a spherical bubble, 0.074  
is surface tension, l  and a  denote density of liquid and gas, respectively. (, ) Ptr is the turbulent kinetic energy distribution over length scales, which can be expressed as follows The bubble breakup time  is the key parameter in the bubble breakup model, which determines the speed of the bubble breakup process. Assuming that a bubble breakup along the midsection into two bubbles, the breakup time can be written as :

1 2 ( , ) exp ( ) 2 2 2 k Ptr r L p p          (8)
() r ra t k   (9)
Considering that surface tension will impede the process of bubble breakup, an additional term 1 ab tt kr is introduced to model the effect of surface tension and avoid divergence near r =0. So the breakup duration can be modeled as follows:

11 () r a rb tt r k       (10)
where r and k represent the characteristic length and characteristic velocity, respectively. a t and b t are case dependent user defined constants . substituting Eq.( 6) into Eq.( 5) we obtain :

    3 (, ) (, ) (, 2 ) 3 () 2 () 3 () ( 2) ftr ftr ft r ra t ra t cr cr t r r           (11)
Using Eqs. ( 4), [START_REF] Levich | Physicochemical Hydrodynamics (English translation by Scripta Tchnica)[END_REF], we can derive the evolution equation of the averaged bubble size:

    3 (, ) (, 2 ) 3 () 2 () 3 () ( 2) 0 r ftr ft r ra t ra t r d r cr cr t r r                     (12)
where

  (, ) (, ) () () () 0 f tr f tr r a t rdr rdr cr r r a cr                 (13)   3 (, 2) 1 (, ) 3 2 () 3 3 () ( 2) 4 0 ft r ftr r a t rdr rdr cr r r a cr                  (14)
Thus Eq. ( 11) can be rearranged as follows

1 (, ) ( 1) 3 () 4 r f tr rdr t r a cr        (15)
substituting Eqs. ( 14) and ( 4) into Eq.( 3) we obtain the source term in the transport equation of the bubble number density:

4 9 1 (, ) ( , ) ( 1) 3 4
() 4 0

f tr m S f t r rdr rdr n r a cr                              (16)
During the simulation, the probability density function of bubble size is given by:

2 1 1 4 3 ( , ) exp 2 3 2 2 f tr r n f f                         ( 17
)
where f  is the standard error and set to be

0.3r f   .
The evolution of the bubble number density can be simulated by solving Eqs. ( 1), [START_REF] Zhang | Investigation of bubble-bubble interaction effect during the collapse of multi-bubble system[END_REF].

Condensation rate based on bubble cluster

In the present work, we are focusing on the bubble-bubble interaction in the collapse stage of cavitation cloud. A condensation rate based on bubble cluster will be constructed in this section.

Dimension analysis of collapse of bubble cluster

Consider a bubble cluster collapsing under pressure p  (Figure 1). Assuming that all bubbles are spherical and in the state of equilibrium at the beginning.

Figure 1 Schematic of bubble cluster

The related parameters are listed as followed: Time: [] tT ; Characteristic length of bubble cluster: [] aL ; Initial bubble radii:

[] 0 RL ;
Population of bubbles per unit volume:

1 [] 3 n L Material parameters: liquid density [] 3 M L L  , liquid viscosity [] M L LT  , surface tension [] 2 M S T
, reference density of gas

[] 0 3 M A L  , gas viscosity [] M A LT



, ratio of specific heats of gas [START_REF] Kubota | Unsteady structure measurement of cloud cavitation on a foil section[END_REF]  ;

Surrounding pressure

[] 2
M p LT  , pressure inside bubbles [] 2 M p B LT ;
The variation rate of volume of all bubbles can be expressed in a relationship: (; , ; ; , , , , , ; , )

0 0 V f tR an S p p LL A A B      0 0 V f t R
Va S L f nR R Rt Rt Rt LL p pp A A B Rt Rt Rt L L LL             V a S f nR ; ( 19
)
where 0 a R can be rewritten as population of bubbles

V S L fN RT RT RT L CL C C p t A A pp RT T L B LC C              V S f N ; (20)
The collapse of bubbles is controlled by pressure difference and inertial force, so viscosity and surface tension can be neglected. For given material

0 A L   and
 are constant. Thus, Equation can be reduced to: Non-dimensional pressure :0.6, 0.7, 0.8, 0.9 p The volume variation rate of bubble cluster is recorded during the calculation and averaged by time. The following discussions are carried out around non-dimensional time-averaged volume variation rate.

Influence of bubble population

Influence of volume fraction

Similarly, we can provide the power-law exponent of  using Curve fitting shows that the power-law exponent of  is 0.160 1 k  . Higher  corresponds to a compacter bubble cluster. Therefore, the inner bubbles will reduce the collapse velocity of outer bubbles. On the contrary, the outer layer of bubbles will shield the inner bubbles from collapsing. As a result, the collapsing velocity will slow down with the rise of volume fraction.

Influence of non-dimensional pressure

Non-dimensional pressure represents the driving force for collapse of bubble cluster. With the increase of p , the collapse will become more violent (Figure 4).

Numerical strategy

The basic approach for cavitating flows consists of unsteady Reynolds Averaged Navier-Stokes equations of the mixture phase and the continuity equation of the vapor phase.

The numerical procedure of the cloud cavitation computation is demonstrated in Figure 5. Unsteady RANS equations combined with modified RNG k- model [START_REF] Reboud | Two-phase flow structure of cavitation: experiment and modeling of unsteady effects[END_REF] are solved by FLUENT using SIMPLEC algorithm to provide background flow information to the continuity equation of vapor phase and transport equation of bubble number density, which are computed in sequence. The calculation of additional transport equation of bubble number density and its source term consumes more computing power, which depends on the mesh size and the resolution of the bubble size distribution function.

Validation

The simulation of the cavitating flow over a cylindrical body with 90° blunt conical head is presented in this section. The diameter of cylinder is d=37.5mm, and the free-stream velocity is v=18m/s. The corresponding cavitation number and Reynolds number are =0.612, Re=6.75×105, respectively. The 2-D axisymmetric geometry is utilized in the simulation, with a 350×150 structural mesh, depicted in Figure 6. The heigh of first grid layer is 2×10 -4 mm, which makes sure that the wall y + is of the order O(1). The numerical results are validated by the experimental data based on the Split Hopkins Pressure Bar (SHPB) launch system [START_REF] Wei | A scaled underwater launch system accomplished by stress wave propagation technique[END_REF]. Figure 9 The full cloud cavitation. Figure 10 The detachment of cloud cavitation.

Figure 11

The collapse of the detached cavitation.

The evolution of cavitation can be divided into four stages: the growth of the sheet cavitation, the development of re-entrant jet, the detachment of cavitation cloud, and the collapse of detached cavitation cloud. When the cavitation occurs, a transparent laminar void like a flat bubble starts to grow from the shoulder of the projectile (Figure 7).

The bubbles then start to break up as the re-entrant jet develops, with a simultaneous increase of the bubble number density. The breakup process turns the sheet cavitation to the cloud cavitation (Figure 8).

The whole cavitation area turns into a non-transparent cavitation cloud when re-entrant jet reaches the shoulder of the projectile 9), where the cavitation cloud detaches and rolls downstream with the mean until collapse (Figure 10, Figure 11). Meanwhile, new attached cavitation appears at the shoulder of projectile. The simulation results show that the formation of cloud cavitation is closely related to the re-entrant jet. Figure 13 presents the distribution of the bubble number density in the cavitation region along the wall of projectile, and the origin of X axis represents the leading edge of cavitation. The cavitation can be divided by re-entrant jet into two parts: the transparent laminar cavitation and the nontransparent cloud cavitation. It is shown that the bubble number density increases rapidly at the front of re-entrant jet and decreases close to the cavitation closure. Bubble size is not uniform in the region of the cloud cavitation (Figure 14), the minimum value of r is 0.202 mm, which corresponds to the maximum n of 3.06×10 10 .

Good agreement with the experimental data suggests that our numerical strategy is capable of predicting the evolution of cavitation with high precision. In addition, the distribution of the bubble number density solved by the proposed model can help us better understand the internal structure of cloud cavitation.

Disscussions

The bubble number density plays an important role in cavitating flows, especially in collapse process of cavitation [START_REF] Reisman | Shock wave measurements in cloud cavitation[END_REF]- [START_REF] Zhang | Investigation of bubble-bubble interaction effect during the collapse of multi-bubble system[END_REF]. There is an important characteristic named geometry focusing of collapse pressure in the collapse process of bubble cluster. The pressure waves emitted by collapse bubbles will propagate inward which leads to a higher magnitude of collapse pressure [START_REF] Wang | Distribution on the dynamics of a spherical cloud of cavitation bubbles[END_REF]. In this paper, we compare the present model with Singhal model to test the ability of predicting pressure focusing effect.

Figure 15 Comparation of condensation rates

Figure 15 shows that present model presents a higher condensation rate than Singhal model. That means present model will predict a higher collapse pressure thant Singhla model (Figure 16). Further more, the condensation rate of Singhal model is quit smooth while present model shows local effects of collapse. The region of cavitation cloud with higher bubble number density collapse more fiercely than other else.

Conlusions

Homogeneous cavitation models can only give volume fractions of different phases. It's not enough for cavitating flow, especially for collapse of cloud cavitation. Different bubble sizes and bubble number density may be corresponding to the same vapor fraction. Cavitation cloud of different internal structure can result in different collapse pressure due to interaction between bubbles. Internal information of cloud cavitation is needed to predict collapse pressure more accurately.

In the present work, a numerical model consists of modified homogenous cavitation model and evolution equation of bubble number density was proposed. The evolution model of bubble number density is governed by a transport type equation and takes account of bubble breakup effect. The condensation rate is constructed through Collapse pressure dimension analysis and direct simulation of collapse of bubble cluster. This model was implemented on cavitating flow over a projectile with blunt conical head and validated by the experiment base on SHPB launch system. The simulation results showed that present model can predict the evolution of cavitation and the distribution of bubble number density (or bubble size) well.

The comparation of present model and Singhal modle was discussed. It is shown that the collapse pressure of detached cavitation cloud is strongly affected by bubble number density. Non-uniform distribution of bubble number density may lead to different collapse velocities. The region of cavitation cloud with higher bubble number density collapse more fiercely than other else.

The proposed model can solve the distribution of the bubble number density and help us better understand the internal structure of cloud cavitation. But, continuous efforts should be made to improve this model. First, only bubble breakup was considered in present model, the effects of coalescence needs to be studied. Second, although the breakup model is simplified, the present model is still time consuming. Further optimization for the numerical method is needed to enhance the computational efficiency.

  length scale which describes the size of large energy-containing eddies in turbulent flows. k is the turbulent kinetic energy, and  is the turbulent dissipation rate.

Figure 2 

 2 Figure 2 Flow field and numerical model Volume of fraction (VOF) method and large eddy simulation (LES) are adopted in the numerical simulations. The parameters are listed as follows: Volume faction : 0.150, 0.268, 0.500  Population of bubbles 3333 : 3 , 4 , 5 , 6 N

Figure 3

 3 Figure 3 Collapse of bubble cluster (  =0.150, N=5 3 , p =0.8 )

Figure 4

 4 Figure 4 Curve of V  corresponding to p Curve fitting gives the power-law exponent 0.609 3 k  . p represents the driving pressure of collapse of bubble cluster. It affects not only the collapse velocity but also the procedure of collapse.

Figure 5

 5 Figure 5 Numerical procedure.

Figure 6

 6 Figure 6 Computational domain. Figure 7 ~ Figure 11 illustrate the evolution of the cloud cavitation. The yellow curves represent the iso-surface of vapor fraction with the value of 0.1, and the red zones in the cavitation indicate high bubble number density area ( 9 10 n ) which can be seen as cloud cavitation.

Figure 7 Figure 8

 78 Figure 7The growth of the sheet cavitation Figure8The development of the re-entrant jet.

Figure 12

 12 Figure 12 Formation of cloud cavitation Figure 12 illustrates the outline of cavitation and the

Figure 13

 13 Figure 13 Distribution of the bubble number density in the cavitation.

Figure 14

 14 Figure 14 Distribution of the bubble size in the cavitation cloud.

Figure 16

 16 Figure 16 Comparation of collapse pressure (at x/d=2)

Table 1~ Table 3

 1~3 give the V  under different bubble populations, where R B represents initial bubble radii, D B is the distance between the centers of two neighbor bubbles.

Table 1 .

 1 V  corresponding to different bubble populations ( =0.150)

	Parameters	1	2	3	4
	N	3×3×3 4×4×4 5×5×5 6×6×6
	R B mm	6.87	5.15	4.12	3.43
	D B mm	20.8	15.6	12.5	10.4
	V 	-10.7	-19.6	-31.4	-45.6

Table 2 .

 2

	Parameters	1	2	3	4
	n	3×3×3	4×4×4	5×5×5	6×6×6
	R B mm	8.33	6.25	5 .00	4.16
	D B mm	20.8	15.6	12.5	10.4
	V 	-9.97	-18.2	-29.0	-42.1

V  corresponding to different bubble populations (=0.268)

Table 3 .

 3

	Parameters	1	2	3	4
	n	3×3×3	4×4×4	5×5×5	6×6×6
	R B mm	10.3	7.69	6.15	5.13
	D B mm	20.8	15.6	12.5	10.4

V  corresponding to different bubble populations =0.500

Table 4 .

 4 V  corresponding to different volume fractions

	N 3×3×3	4×4×4	5×5×5	6×6×6
	0.15	-10.7	-19.6	-31.4	-45.6
	0.268	-9.97	-18.2	-29.0	-42.1
	0.5	-9.05	-16.4	-25.9	-37.3

Consider average variation rate of volume, the non-dimensional time will not appear in the formulation:

where

represents volume variation rate of single bubble, therefore

volume variation rate of bubble cluster. We can see from the Eq. ( 22) that three key parameters:  , N and p ,To provide the specific relation, we can assume that:

The empirical parameters ,, , 123 ck k k can be determined by direct simulation of collapse of bubble clusters.

Direct simulation of collapse of bubble cluster

Consider a cubic bubble cluster in which bubbles are arranged in nnn  order as shown in Figure 2. The initial pressure inside bubble cluster is set to be p B , and the pressure at boundary is p  .

Condensation rate

The constant coefficient c is closed to 1, combined with the influence of  , N and p , Eq. ( 23) becomes :

Equation ( 24) gives the no-dimensional collapse velocity of bubble cluster. Then the change rate of void fraction can be given as:


Renfang Huang

An Yu

Ji Bin

Jiajian Zhou

Zhihong Zhai

Xianwu Luo

Cavitating flow features in a water-jet pump under different upstream conditions

Keywords:

Turbulent cavitating flows in a water-jet mixed-flow pump were investigated by using the k- SST turbulence model and a cavitation model based on mass transport equation to give an insight into cavitating features under different upstream conditions. Both the hydraulic and cavitation characteristics predicted by the present numerical method agree well with the available experimental results. The numerical results depict that the cavity in the pump firstly occurs at the blade leading edge close to blade tip, and then grows both downstream and toward the impeller hub. With cavitation development, sheet cavities expand along the impeller shroud, and are likely to block the flow passage under the critical cavitation condition. Besides, the thrust of the water-jet pump is enhanced gradually during cavitation evolution, but the hydraulic performance of the diffuser remains unchanged. The hydraulic loss at the impeller suction surface increases with secondary flows and Boundary Vorticity Flux (BVF) peak during cavitation, which is located at the vapor-liquid interface, especially near the blade tip.

INTRODUCTION

Nowadays, more and more interest has been attracted by the water-jet propulsion which has a variety of advantages such as less vibration, low noise, high efficiency, etc. in the specified speed range [START_REF] Kandasamy | Optimization of water-jet propelled high speed ships[END_REF][START_REF] Carlton | Marine propellers and propulsion[END_REF]. In recent years, the water-jet propulsion has been rapidly developed and widely used in the marine field, which is expected to be the potential technology against the propeller due to the improvement in the progress of pump design, manufacture and fabrication [START_REF] Takai | Verification and validation study of URANS simulations for an axial water-jet propelled large high-speed ship[END_REF]. Many scholars [START_REF] Kandasamy | Simulation based design optimization of water-jet propelled Delft catamaran[END_REF][START_REF] Kandasamy | Numerical and experimental evaluation of water-jet propelled delft catamarans[END_REF][START_REF] Tahara | CFD-based multiobjective optimization of water-jet propelled high speed ships[END_REF][START_REF] Bulten | Numerical analysis of a water-jet propulsion system[END_REF] have investigated the water-jet propulsion in different ships using the computational fluid dynamics (CFD) method or the model test according to the procedure proposed by the ITTC Water-jet Performance Prediction Specialist Committee [START_REF] Van | Report of the specialist committee on validation of water-jet test procedures[END_REF][START_REF] Kim | The propulsion committee: final report and recommendations to the 25th ITTC[END_REF]. Kandasamy et al. [START_REF] Kandasamy | Simulation based design optimization of water-jet propelled Delft catamaran[END_REF][START_REF] Kandasamy | Numerical and experimental evaluation of water-jet propelled delft catamarans[END_REF], and Tahara et al. [START_REF] Tahara | CFD-based multiobjective optimization of water-jet propelled high speed ships[END_REF] optimized the water-jet propulsion high speed ships using both the simulation-based design method and model test. Bulten et al. [START_REF] Bulten | Numerical analysis of a water-jet propulsion system[END_REF] studied a water-jet propulsion system, and provided a reference for other researchers.

Up to now, many studies focused on the propulsion system optimization and characteristic prediction without considering cavitation in the water-jet pump. In fact, cavitation is inevitable when the water-jet pump is operated at a high speed which may result in pressure fluctuation, noise and erosion [START_REF] Huang | Numerical investigation of cavitation-vortex interaction in a mixed-flow waterjet pump[END_REF][START_REF] Luo | Design and optimization for fluid machinery[END_REF]. Thus, good cavitation performance as well as high hydrodynamic efficiency is considered as the necessary criterion to design and evaluate a water-jet pump.

In general, experiments and numerical simulation are very useful tools to evaluate the pump performance. Until now, compared with the long-time period and huge investment of the experiment, numerical simulation has been wide-spread used due to the progress of computer technology, which can provide a basic understanding for the evaluation and improvement of hydraulic machines [START_REF] Luo | Impeller inlet geometry effect on performance improvement for centrifugal pumps[END_REF][START_REF] Ji | Unsteady numerical simulation of cavitating turbulent flow around a highly skewed model marine propeller[END_REF][START_REF] Luo | Effect of impeller inlet geometry on centrifugal pump cavitation performance[END_REF][START_REF] Ji | Numerical analysis of cavitation evolution and excited pressure fluctuation around a propeller in nonuniform wake[END_REF].

In this paper, the turbulent cavitating flows in a water-jet mixed-flow pump was investigated under different upstream conditions by solving Reynolds Averaged Navier-Stokes (RANS) equations and a cavitation model [START_REF] Luo | Comparison of cavitation prediction for a centrifugal pump with or without volute casing[END_REF]. Based on the calculated results, cavitating flow features were further discussed.

Geometrical parameters of the water-jet pump

The water-jet mixed-flow pump is composed of a 6-blade impeller, a 7-blade diffuser, an inlet flow channel and an outlet flow channel. The pump was designed with flow coefficient of =0.183 at the rotational speed of n=1450 r/min.

The meridional schematic diagram for the water-jet pump is shown in Figure 1 with geometrical parameters listed in Table 1. The Reynolds number Re, is 2.47×10 6 , and the blade tip clearance relative to the impeller diameter is 0.86×10 -3 .

  m j m j 0 u tx       (1)       j i mm m ii j i jj j m i j t x x p u uu u x xx uu                  (2)
where p, u are the pressure and velocity; m, m are density and dynamic viscosity of the mixture, and defined as follows:

  l vv mv      1 (3)   l vv mv      1 ( 4
)
where subscript v, l represent the vapor and liquid, respectively.

The k- SST turbulence model is adopted to obtain highly accurate predictions of flow separation under adverse pressure gradients thanks to its consideration of the turbulent shear stress.

The cavitation is modeled by mass transfer equation: derived from the cavitation model [START_REF] Zwart | A two-phase flow model for predicting cavitation dynamics[END_REF].

  v v vi e c v umm t             (
Then, the Rayleigh-Plesset equation is applied to solve the source terms and control vapor generation and condensation.

    nuc v v v vap e nuc l max , 0 31 2 3 p p r m F R         (6)   v v v cond c nuc l max , 0 3 2 3 pp m F R        (7)
To obtain the interphase mass transfer rate, four model parameters are given as follows: bubble radius Rnuc = 1 m, the nucleation volume fraction rnuc = 510 -4 , the evaporation coefficient Fvap = 50 and the condensation coefficient Fcond = 0.01.

Mesh generation and independence test

In this study, one blade-to-blade flow channel are considered for simulation. The computational domain is divided into the suction pipe, the outlet pipe, the impeller and the diffuser. Two interfaces are set in the suction and outlet pipe to monitor the total pressure. The structured hexahedral mesh is applied to obtain accurate simulation results. The topological structure and mesh point distribution are the same in different models. What is more, mesh around the impeller blade are refined to satisfy the requirements for the turbulence model with the Yplus ranging from 4 to 90 in Figure 2.

Figure 3 shows the mesh indepandence test for mesh generation. Note that the pump head coefficient and efficiency are more sensitive to the mesh density of diffuser than that of impeller. When the mesh quantity in impeller and diffuser reach 400,000 respectively, the pump head coefficient and efficiency remain unchanged. Therefore, 460,000 in impeller and 400,000 in diffuser is decided as the final mesh scheme with considering the computing resources. The impeller is a rotating part and other parts are static, so the multiple reference frame (MRF) is used in this simulation and the interfaces between the rotating part and its neigboring static part are set to be frozen-rotor under steady condition.

For boundary conditions, the mass flow rate is given at the outlet plane based on the mass equilibrium. The total pressure is applied at the inlet plane and gradually reduced when cavitation develops. All solid walls are set as the nonslip wall condition.

RESULTS AND DISCUSSION

Figure 4 shows the characteristic curves of the pump. The calculation results are marked as "cal.", and the experimental results are marked as "exp.". The results demonstrates that the calculated head coefficients agree well with the experimental data although the calculated efficiency is always a little higher than that by experiments. When the water-jet pump is operated at the design point with flow coefficient =0.183 and rotational speed n=1450 r/min, the head coefficient is 0.33 in calculation and 0.31 in experiment, and the efficiency is 0.89 in calculation and 0.83 in experiment. This discrepancy in efficiency is reasonable due to the fact that the efficiency in numerical simulation is the hydraulic efficiency, the effects of the mechanical and volumetric losses of the pump are not considered.

Net Positive Suction Head (NPSH) is usually used to describe the cavitation performance for pumps. The definition of NPSH is shown subsequently:

tv ll pp NPSH g g   ( 8
)
where pt is the total pressure at the inlet plane in Figure 2. Thoma's cavitation number i.e.  is a dimensionless cavitation performance parameter, and can be defined as the following:

NPSH H   (9)
where H is the pump head at the design point. The critical NPSH, i.e., NPSHc, is defined as the NPSH value when the pump head drops 3% from its original value without cavitation.

Figure 5 shows cavitation characteristics, where both the numerical and experimental results are included for comparison. The critical cavitation number is 0.36 predicted by calculation and 0.49 by experiment.There is a difference in predicting cavitation characteristics by present numerical method since many factors in experiments, such as gas in local water and water temperature, can't be taken into account. However, the simulation shows the same tendency of cavitation effects on pump performance as the experiment does.

Usually, the marine water-jet propulsion vessel is required to move at the desirable speed, so the thrust is an important performance for the water-jet pump. Figure 6 shows the thrust T and vapor volume under different cavitation conditions with flow coefficient =0.183 and rotational speed n=1450 r/min. The thrust firstly slightly changes and then increases in a jump near the critical cavitation condition. Table 2 presents the mass-averaged velocity at the inlet plane, impellet outlet plane and diffuser outlet plane to analyze the pump thrust T=Q(V2z.diffuser-Vin). As the inlet velocity Vin is basicaly constant, the axial velocity at the diffuser outlet planeV2z.diffuser directly affects the thrust. It is indicated that the hydraulic performance of the diffuser hardly changes with cavitation development, and the thrust is essentially connected with cavitating flows in impeller.

Cavitation evolution in impeller is shown in Figure 7. Cavity illustrated by iso-surface of vapor volume fraction=0.1 firstly occurs at the blade leading edge and close to the blade tip, then grows both downstream and toward the impeller hub. When cavitation number is 0.43, the maximum cavity extends to 2/3 blade length with scattered cavities attaching to the impeller shroud before the blade inlet. Gradually, the scattered cavities develop into sheet cavities surrounding the shroud and extend toward the blade-to-blade passage under the critical cavitation condition. The flow passage blockage due to cavitation development causes an increase in the velocity at the impeller outlet plane as shown in Table 2, and results in a thrust jump.

Figure 8 shows the secondary flow contours and streamline near impeller suction surface. Boundary Vorticity Flux (BVF) [START_REF] Wu | Vorticity and vortex dynamics[END_REF] is introduced to diagnose the flow separation, and the results are shown in Figure 9. Note that BVF is defined as

zl / BVF n p   
, where z n is local axial normal direction at the wall. Figure 10 shows loss distribution at the impeller suction surface. Note that the loss is defined as

  2 T li n =/ / 2 pp u  
, where l p is the local total pressure, in p is the total pressure at the impeller inlet plane and T u is tangential velocity.

Based on those results, the following features can be seen:

(1) In non-cavitating flows, i.e. the case of  =2.72, the secondary flow appears at the leading edge (LE) near the impeller hub, and BVF has peak values at both leading edge (LE) and trailing edge (TE). Large loss occurs at TE-hub corner. It is noted that this loss is not affected by cavitation.

(2) When cavitation occurs, i.e.  =0.89, the secondary flow seems to expand from LE-corner toward blade shroud, and BVF peak line at the shroud side is likely to bend downstream compared with that in non-cavitating condition. A new loss area appears at the vapor-liquid interface near the blade hub.

(3) As cavity reaches the hub and the maximum cavity extends 1/2 blade length (=0.43), the secondary flow seems enhanced at the LE-hub corner, and another secondary flow occurs near the cavity frontier. The streamlines mixes along the cavity-liquid interface.

(4) When cavitation fully develops near the critical condition (=0.36), the secondary flow as well as BVF peak disappears at the LE-hub corner due to the enclosure by cavity, but the secondary flow and BVF are strong near the cavity-liquid interface. The mixed streanline approaches the blade midspan, and causes a jump of hydraulic loss increase.

Conclusions

The cavitating flows in a water-jet mixed-flow pump are simulated by using the k- SST turbulence model and a mass transport equation based cavitation model. The numerical results were compared with the available experimental data, and the cavitating flow features together with loss mechanism were dicussed in different cavitating conditions. According to those results, the following can be concluded:

(1) Both the hydraulic and cavitation performance were reasonably predicted by the present numerical method. Thus, the numerical method is suitbale for the cavitating flow analysis in a water-jet pump.

(2) The thrust of water-jet pump increases gradually with cavitation development since the flow velocities in impeller as well as jet velocities increased. However, hydraulic performance of the diffuser hardly change due to cavitation development.

(3) Cavity in the pump firstly occurred at the blade leading edge near the blade tip, and then grew both downstream and toward the impeller hub. For cavitation at the critical condition, sheet cavities expand along the impeller shroud, extend toward the blade-to-blade flow passage, and cause the flow passage blockage.

(4) Secondary flows as well as BVF peak cause no loss at the leading edge under non-cavitation condition. Loss at the TE-hub corner is not affected by cavitation and partially results from the BVF peak at TE.

(5) As cavitation develops, hydraulic loss at the impeller suction surface increases along with secondary flows and BVF peak, which is located at the vapor-liquid interface, especially near the blade tip.

!Figure 1 .

 1 Figure 1. Meridional schematic diagram for the pump

  5) wherev, v respectively represent the vapor phase density and volume fraction. The source terms e m   and c m   are

Figure 2 .Figure 3 .

 23 Figure 2. (a) Computation domainand (b) Yplus distribution at the impeller blade (left: suction surface (SS); right: presure surface (PS)) SS PS

Figure 4 .Figure 5 .Figure 6 .Figure 7 .Figure 8 .Figure 9 .Figure 10 .

 45678910 Figure 4. Characteristic curves for the pump (n=1450 r/min)

m

  Densityof the homogenous phase v Density of the vapor phase  Thoma's cavitation number, =NPSH/H  Flow coefficient  Head coefficient

Table 1 .

 1 Geometrical parameters for the impeller

	Parameters	Symbols Values
	Blade width ratio at the inlet	b1/D2	0.318
	Blade width ratio at the exit	b2/D2	0.247
	Blade inlet diameter at the tip	D1/D2	0.856
	Blade inlet diameter at the hub	D1h/D2 0.245
	Blade exit diameter at the hub	D2h/D2 0.604
	Hub ratio	Dh/D2	0.209
	Blade tip clearance	c/D2	0.86e-3

2. Numerical method 2.1 Governing equations and cavitation model

A homogenous assumption is used to solve the cavitation flow, where the liquid and vapor are considered as a mixture, and share the same velocity. The mixture flow can be described by the Reynolds Averaged Navier-Stokes (RANS) equations as follows:

ACKNOWLEDGMENTS

This work was financially supported by the National Natural Science Foundation of China (Project Nos. 51376100, 51306018 and 51536008), and State Key Laboratory for Hydroscience and Engineering (Project Nos. 2015-E-03 and 2014-KY-05).

NOMENCLATUR b1

Blade width at the impeller inlet b2

Blade


Lorenzo Codecasa
email: codecasa@elet.polimi.it

Dario D'amore
email: damore@elet.polimi.it

Paolo Maffezzoni

Nadia Spennagallo
email: spennaga@elet.polimi.it

A Novel Approach for Generating Dynamic Compact Models of Thermal Networks Having Large Numbers of Power Sources

A novel approach is presented for generating compact dynamic thermal networks having large numbers of power sources. The achievable accuracy of the compact model is controllable. The number of elements of the compact model linearly increases with the number of power sources.

Introduction

Various techniques are reported in literature for generating compact models whose responses approximate the responses of large thermal networks having multiple power sources [START_REF] Lasance | Creation and Evaluation of Compact Models for Thermal Characterization Using Dedicated Optimisation Software[END_REF][START_REF] Batty | Electrothermal CAD of power devices and circuits with fully physical time-dependent compact thermal modeling of complex non linear 3-D systems[END_REF][START_REF] Rencz | Dynamic Thermal Multiport Modeling of IC Packages[END_REF][START_REF] Sabry | Dynamic Compact Thermal Models: An Overview of Current and Potential Advances[END_REF][START_REF] Codecasa | Compact modeling of electrical devices for electrothermal analysis[END_REF]. With such techniques, a compact model can be determined by modeling the thermal problem as a whole. Usually, very accurate compact models are generated in this manner. However they have a number of components quadratically increasing with the number n of power sources in the thermal network. That is because distinct components are introduced in the compact models in order to represent each interaction among the n power sources. Thus the resulting compact models are actually compact only for small n. These techniques cannot be applied to generate compact models of thermal networks having large numbers of power sources. For instance, these techniques cannot be applied for generating compact models, at the transistor level, for electronic circuits with many electrical devices and, at the package levels, for modules with many dies.

However, a compact model of a thermal network having a large number of power sources can be obtained as proposed in [START_REF] Batty | Electrothermal CAD of power devices and circuits with fully physical time-dependent compact thermal modeling of complex non linear 3-D systems[END_REF][START_REF] Sabry | Dynamic Compact Thermal Models: An Overview of Current and Potential Advances[END_REF]. Firstly the spatial domain of the thermal problem is partitioned into subdomains in such a way that only few power sources are in each sub-domain. Secondly the thermal problem in each of these sub-domains is represented by a boundary conditions independent compact model. Lastly, by interconnecting the obtained boundary conditions independent compact models, a compact model of the whole thermal network is obtained. Since each boundary condition independent compact model has a number of components independent from n, and since the number of the boundary condition independent compact models is proportional to n, the composite model has a number of components proportional to n. However, this com-pact model is not assured to be accurate, since the boundary conditions accurately modeled by each boundary condition independent compact model, in general cannot model the actual boundary conditions of the thermal problem.

In this paper a novel technique is proposed for generating dynamic compact models of discretized thermal networks having large numbers of power sources. A compact model generated in this manner assures a desired accuracy. Moreover it has a number of components proportional to the number n of the power sources. Such a compact model is the interconnection of boundary condition independent compact thermal networks. Each boundary condition independent compact thermal network is generated in such a way that a desired accuracy is assured. This is achieved by exploiting the projection vectors of the Multipoint Moment Matching Method. In this sense the novel method is an hybrid of the approach based on boundary condition independent models and of the Multipoint Moment Matching approach.

The rest of this paper is organized as follows. In section 2, dynamic thermal networks are introduced. The Multipoint Moment Matching method and the use of boundary condition independent models are revised in sections 3 and 4. The novel hybrid method is presented in section 5. The application of this technique to a power DMOS is presented in section 6.

Dynamic Thermal Networks

We refer to a generic heat diffusion problem in the bounded spatial region Ω. The relation between the power density F (r, t) and the temperature rise x(r, t) with respect to the ambient temperature is modeled by the heat diffusion equation

c(r) ∂x ∂t (r, t) + ∇ • (-k(r)∇x(r, t)) = F (r, t) (1)
in which k(r) is the thermal conductivity and c(r) is the volumetric heat capacity. Boundary conditions are

-k(r) ∂x ∂ν (r, t) = h(r)x(r, t) (2)
in which h(r) is the heat exchange coefficient and ν(r) is the unit vector outward normal to the boundary ∂Ω.

A dynamic thermal network N is introduced as follows. The power density is limited to the form

F (r, t) = n i 1 f i (r)P i (t) (3)
in which P i (t), i-th element of the n×1 vector P(t), is the power of the i-th power source measured at the i-th port of the thermal network and f i (r) is a volumetric density function of support Ξ i . As in [START_REF] Codecasa | Compact modeling of electrical devices for electrothermal analysis[END_REF] it is assumed

T i (t) = Ω f i (r)x(r, t) dr, (4)
in which T i (t), i-th element of the n×1 vector T(t), is the temperature rise of the i-th power source measured at the i-th port of the thermal network.

The only assumption on the f i (r) volumetric density functions is that they are linearly independent.

In this manner the dynamic thermal network has the generalized form introduced in [START_REF] Sabry | Higher order compact thermal models[END_REF][START_REF] Codecasa | Multivariate moment matching for generating boundary condition independent compact dynamic thermal networks of packages[END_REF]. This dynamic thermal network preserves the main physical properties of the heat diffusion problem, in particular passivity. The discretization of Eqs. (1), (2) by any proper discretization technique has the form

C dx dt (t) + K x(t) = f (t) (5)
in which x(t) is the M ×1 vector of freedom degrees, f (t) is the M ×1 source vector, C and K are M ×M matrices. The discretized forms of Eqs. ( 3) and ( 4) are

f (t) = FP(t), (6)
T(t) = F T x(t) (7)
in which F is an M × n matrix. Eqs. ( 5)-( 7) define a discretized dynamic thermal network. Also this discretized dynamic thermal network preserves the main physical properties of the heat diffusion problem, in particular passivity.

A compact dynamic thermal network is a lumped network that is simpler than a discretized dynamic thermal network but is able to accurately approximate the relation between the powers and the temperature rises at the ports of a dynamic thermal network. Hereafter the problem of generating a compact dynamic thermal network from a discretized dynamic thermal network is tackled.

Compact Dynamic Thermal Networks by Multipoint Moment Matching

A very accurate and efficient technique for generating a compact dynamic thermal network is the multipoint moment matching method [START_REF] Codecasa | Compact modeling of electrical devices for electrothermal analysis[END_REF], hereafter revised. In this technique, x(t) is approximated by

x(t) ≈ m j 1 v j xj (t) (8)
in which v j are basis vectors and xj (t) are freedom degrees, with j = 1, . . . m. Equivalently

x(t) ≈ Vx(t) (9)
in which

V = [v 1 , . . . , v m], x(t) =    x1 . . . x m    .
Firstly the basis vectors v j , with j = 1, . . . , m, are determined. To this aim the q problems

(α r C + K)V r = F
are solved for the M × n rectangular matrices V r , with r = 1, . . . q. Here the α r matching points are

α r = λ M dn 2r -1 2l K, k (10)
in which K is the complete elliptic integral of first kind of modulus k, dn is the elliptic dn function of modulus k. Besides

k = 1 -k 2 k = λ 1 /λ M ,
in which λ 1 and λ M are the minimum and maximum eigenvalues of matrix

C -1 K. It is set V = [V 1 , . . . , V q ]
having m = nq columns. Secondly the freedom degrees xj (t) with j = 1, . . . , m are determined by means of Galerkin's method. To this aim Eqs. ( 5)-( 7) are projected onto the space spanned by the vectors v j with j = 1, . . . , m. It results in

Ĉ dx dt (t) + Kx(t) = f (t), (11)
with

f (t) = FP(t) (12)
T(t) = FT x(t), (13)
in which

Ĉ = V T CV, (14)
K = V T KV, (15) F = V T F. (16)
Eqs. ( 11), (13) determine x(t).

As proved in [START_REF] Codecasa | Parameters for multi-point moment matching reduction of thermal networks[END_REF], with this choice of the basis vectors v j and of the freedom degrees xj , with j = 1, . . . , m, the expansion given by Eq. ( 9) approximates x(t), with a relative error in the energy norm ε exponentially decreasing with q. Thus q is of order O(-log ε).

As a result multipoint moment matching can be interpreted as an effective method for generating basis vectors v j , with j = 1, . . . , m, in such a way that a linear combination of such basis vectors can approximate x(t) with an error exponentially converging to zero with the number m of basis vectors.

Eqs. ( 11), (13) accurately approximate also the relation between the port powers P i (t) and the port temperature rises T i (t), with i = 1, . . . , n. Thus they define a compact dynamic thermal network N of state-space dimension m. This compact dynamic thermal network preserves the main physical properties of the heat diffusion problem, in particular passivity. Also such compact dynamic thermal network can be synthesized by means of passive elements proceeding as follows. The generalized eigenvalue problem

K Û = Ĉ Û Λ, ( 17
) ÛT Ĉ Û = 1, (18)
is solved for the eigenvectors matrix Û and for the diagonal eigenvalue matrix Λ. The transformation of variables is then introducing

x(t) = Ûâ(t),
by which Eqs. ( 11), ( 13) are transformed into the equivalent form

dâ dt (t) + Λâ(t) = ΓP(t), (19)
T(t) = ΓT â(t) (20)
where

Γ = ÛT F.
By interpreting Eqs. ( 19), (20) as the modified nodal analysis equations of a thermal network, the compact dynamic thermal network is synthesized by the canonical form [START_REF] Codecasa | Compact modeling of electrical devices for electrothermal analysis[END_REF] shown in Fig. m × n ideal transformers. Thus the number of elements in such a thermal network is O(n 2 ), preventing its use for large n.

Compact Dynamic Thermal Networks by Boundary Condition Independent Models

As shown in [START_REF] Codecasa | Compact Thermal Networks for Modeling Packages[END_REF] the dynamic thermal networks introduced in section 2 can be generated in such a way that they mimic boundary condition independent models. This is achieved by assuming that the F (r, t) function in Eq. ( 1) does not only model a volume power density but also models a surface power density over the boundary ∂Ω and by assuming that the heat exchange coefficient h(r) in Eq. ( 2) is zero on the boundary ∂Ω. In this way some of the f i (r) functions of Eqs. ( 3), ( 4) are assumed to be surface power density functions whose supports are on the boundary ∂Ω. With a proper choice of the surface density functions such a dynamic thermal network mimics a boundary condition independent model. It is thus referred to as a boundary condition indepedent dynamic thermal network. Such a boundary condition indepedent dynamic thermal network can still be discretized as shown in section 2. Besides the multipoint moment matching method can still be used with a slight modification for generating a boundary condition independent compact dynamic thermal network from the boundary condition independent discretized dynamic thermal network. Precisely the V matrix is assumed as

V = [e, V 1 , . . . , V q ]
in which e is an M × 1 vector of ones. Moreover k = λ 2 /λ M , being λ 1 = 0. This compact thermal network still admits the canonical form of section 3. One resistor has conductance λ1 = 0.

By generalizing the approach proposed in [START_REF] Batty | Electrothermal CAD of power devices and circuits with fully physical time-dependent compact thermal modeling of complex non linear 3-D systems[END_REF], such boundary condition independent compact dynamic thermal networks can be exploited for generating a compact model of the dynamic thermal network N of Eqs. ( 1)-( 4). Precisely 1. The Ω region is partitioned into the parts Ω i , with i = 1, . . . , v. The non-empty intersections of these parts define the surfaces Σ i , with i = 1, . . . , f , as shown in Fig. 2.

2. The volume power densities in Eqs. ( 3), ( 4) are introduced. Proper surface power densities having supports on the Σ i surfaces, with i = 1, . . . , f , are also introduced.

3.

The Ω i parts, with i = 1, . . . , v, are modeled by boundary condition independent compact dynamic thermal networks Ni , with i = 1, . . . , v.

Each Ni , with i = 1, . . . , v, has one port for each volume power density whose support is in Ω i and has one port for each surface power density whose support is on the boundary ∂Ω i .

4. The ports which correspond to one surface power density are interconnected, as shown in Fig. 2.

Since the connection of passive thermal networks is passive the determined composite thermal network is passive and thus also stable.

If the partition is such that each of the Ω i parts, with i = 1, . . . , v, contains only few of the supports Ξ i , with i = 1, . . . , n, and has on the boundary only few of the surfaces Σ i , with i = 1, . . . , f , then the generated compact dynamic thermal network has only O(n) components and thus can be used also for large n. As it is well known, the drawback of this approach is that the achievable accuracy is not well controllable.

Compact Dynamic Thermal Network by Hybrid Approach

A compact dynamic thermal network composed of boundary condition independent models is accurate only if the restrictions of the heat flux density of Eqs. ( 1)-( 4) to the Σ i surfaces, with i = 1, . . . , f , are accurately approximated by linear combinations of the surface power density functions. This PSfrag replacements is assured if the surface power density functions are the restrictions to the Σ i surfaces, with i = 1, . . . , f , of the heat flux densities due to the temperature rises v j , with j = 1, . . . nq, determined by means of the multipoint moment matching method for Eqs.

P1(t) T1(t) N1 P2(t) T2(t) P3(t) T3(t) N2 N3 P4(t) T4(t) P5(t) T5(t) P6(t) T6(t) P7(t) T7(t) Ξ 1 Ξ 2 Ξ 3 Σ 1 Σ 2 Σ 3 Ω 1 Ω 2 Ω 3 ∂Ω (a) (b)
(1)-( 4). Let Q k be the matrix whose columns are the restriction of such heat flux densities to the Σ k surface, with k = 1, . . . , f . The question of approximating the linear combinations of the nq columns of Q k by linear combinations of p k nq vectors is here considered. It is well known [START_REF] Golub | Matrix computations[END_REF] that by determining the Singular Value Decomposition of

Q k , Q k = U k Σ k W T
k , and by choosing the p k vectors as the first p k columns of U k , the relative error in the energy norm reaches its minimum given by the p k + 1-th singular value. The following procedure for choosing the surface power density functions is thus deduced. The first p k step of the Singular Value Decomposition of Q k are performed, stopping when the p k +1th singular value is smaller that the desired accuracy ε. The surface power density functions are then chosen equal to the first p k columns of U k . The Multipole Algorithm [START_REF] Greengard | The Rapid Evaluation of Potential Fields in Particle Systems[END_REF] suggests that for a heat diffusion problem the singular values exponentially decrease with p k . As a result only O(-log ε) surface power density functions are introduced in this manner for each Σ k surface, with k = 1, . . . , f . Similar considerations can be made for determining the projection vectors for generating the boundary condition independent compact dynamic thermal networks. Let V k be the matrix whose columns are the restriction to the Ω k region, with k = 1, . . . , v, of the temperature rises v j , with j = 1, . . . nq, determined by means of the multipoint moment matching method for Eqs. ( 1)-( 4). The question of approximating the linear combinations of the nq columns of V k by linear combinations of r k nq vectors is considered. Again by detemining the Singular Value Decomposition of V k ,

V k = U k Σ k W T
k , and by choosing the r k vectors as the first r k columns of U k , the relative error in the energy norm reaches its mimimum given by the r k + 1-th singular value. The following procedure for choosing the projection vectors is thus deduced. The first r k step of the Singular Value Decomposition of V k are performed, stopping when the r k + 1-th singular value is smaller that the desired accuracy ε. The projection vectors are chosen equal to the first r k columns of U k . Again the Multipole Algorithm [START_REF] Greengard | The Rapid Evaluation of Potential Fields in Particle Systems[END_REF] suggests that for a heat diffusion problem the singular values exponentially decrease with r k . As a result only O(-log ε) projection vectors are introduced in this manner for each Ω k region, with k = 1, . . . , v. Since only O(-log ε) surface power density functions are introduced for each Σ k surface, with k = 1, . . . , f , and only O(-log ε) projection vectors are introduced for each Ω k region, with k = 1, . . . , v, the generated compact dynamic thermal network has only O(n) components and thus can be used also for large n. Moreover the achievable accuracy ε is controllable.

Numerical Example

The heat diffusion problem for two fingers of a power DMOS is considered. All 96 channel re- gions are assumed as independent power sources. Firstly a compact dynamic thermal network is generated by means of the Multipoint Moment Matching method with ε = 10 -4 . The relative error, both in the time and in the frequency domains, is below 0.1%. The number of elements in the thermal network is 112 896.

Secondly a compact dynamic thermal network is generated by means of a boundary condition independent model. The Ω region is partitioned into the parts Ω k , with k = 1, . . . , 96 as shown in Fig. 3. A single uniform surface power density function is introduced for each Σ k surface, with k = 1, . . . , 164. The relative error, both in the time and in the frequency domains, is above 15%, as shown in Figs. 4,5. The number of elements in the thermal network is 33 024.

Lastly a compact dynamic thermal network is generated by means of the hybrid method with ε = 10 -4 . The partition of Fig. 3 is still used. The relative error, both in the time and in the frequency domains, is below 0.1%, as shown in Figs. 6,7. The number of elements in the thermal network is 12 832.

Conclusion

In this paper a novel approach has been presented for generating compact dynamic thermal networks modeling heat diffusion problems with large numbers of power sources. As in the multipoint moment matching approach, the achievable accuracy of the model is controllable. As in the approach based on boundary condition independent approach, the number of elements of the compact thermal network linearly increases with the number of power sources.

Figure 1 :

 1 Figure 1: Equivalent network of an n-port compact dynamic thermal network N . Parameter γij is the coefficient at the i-th row and j-th column of Γ.

Figure 2 :

 2 Figure 2: (a) Partition of the Ω region. (b) Compact Dynamic Thermal Network Composed of Boundary Condition Independent Models.

1 Ξ 1 Figure 3 :

 113 Figure 3: Partition of the Ω region.

Figure 4 :

 4 Figure 4: Power step thermal responses of the compact dynamic thermal network composed of boundary condition independent models.

Figure 5 :

 5 Figure 5: Nyquist plot of thermal impedances of the compact dynamic thermal network composed of boundary condition independent models.

Figure 6 :

 6 Figure 6: Power step thermal responses of the compact dynamic thermal network derived by the hybrid method.

Figure 7 :

 7 Figure 7: Nyquist plot of thermal impedances of the compact dynamic thermal network derived by the hybrid method.