Zbigniew Suszyński

Radosław Duer
email: 2rduer@man.koszalin.pl

Mateusz Kosikowski

THE METHOD OF NON-LINEAR DISTORTIONS ELIMINATION IN PHOTOACOUSTIC INVESTIGATION OF LAYERED SEMICONDUCTOR STRUCTURE

The method of non-linear distortions elimination in photoacoustic investigation of layered semiconductor structure

Z. Suszynski, R. Duer, M.

INTRODUCTION

Magnificent tools for the quality of adhesion control are thermal wave methods. Among them, the photo-acoustic (PA) technique in broad-band modulation mode can be distinguished. For this technique, The photo-acoustic phenomenon is being utilitized within. The modulated laser beam is used to excite a variable heat flux inside the structure so it changes the pressure of constant volume of gas closed inside photo-acoustic chamber. The sound generated as a result of pressure changes is being registered with a microphone. This technique allows measuring the temperature response fast and with relatively low costs [START_REF] Quimby | Real-time photo-acoustic microscopy[END_REF][START_REF] Tam | Applications of photoacoustic sensing techniques[END_REF][START_REF] Mandelis | Generalized methodology for thermal diffusivity depth pofile reconstruction in semi-finite and finitely thick in-homogeneous solids[END_REF][START_REF] Rosencwaig | Thermal wave imaging with thermo-acoustic detection[END_REF].

However, the PA technique, as any other indirect method of temperature measurement is characterized by high level of linear distortions [START_REF] Suszyński | Termofalowe metody badania materiałów i przyrządów elektronicznych[END_REF]. Thus the correct analysis of thermal properties basing upon the photoacoustic measurements is possible under two conditions: the reduction of linear distortion influence and elimination of non-linear distortions.

Fig.1. General scheme of photo-acosutic technique in broad-band excitations mode.

The first condition can be assured using a comparative method of analysis, for example the contrast method. Regarding the second condition, consideration of the presence and the influence of non-linear distortions as well as the elimination of that are presented in this contribution.

THE OBJECT OF INVESTIGATION AND THE REFERENCE OBJECT

High power thyristor structure was used as an object of investigation and a sample of passivated glass was assumed as a reference object. The investigated object consisted of 3 layers: silicon -eutectics -molybdenum. The thickness and thermophysical properties of each layer are presented in Tab. 1.

In the case of weak adhesion between layers, an additional layer of high thermal impedance and strongly diversified thermal parameters can be formed. It was assumed that it is a layer of air.

THE EXPERIMENT

The experiment was performed using the photo-acoustic measurement set-up presented in Fig. 3. Photo-acoustic (PA) images were acquired in pulse excitations mode. This mode allows significantly reduce the time of experiment because the spectra of the temperature response contain in this case complete information about thermal properties of object. Laser diode was used as a source of excitation. The intensity of a laser beam was modulated by a rectangular pulse (duration t 0 = 0.1s, period T = 2.5s). PA signal was registered synchronously with switching the laser on with time resolution of 1ms. The registration of PA signal was performed for both the investigated and the reference objects. Fig. 4 presents the PA signal registered during the experiment. In order to investigate the presence of non-linear distortions as well as their possible influence on PA signal, PA responses were registered for different values of laser diode current ranging from 0.86A to 1.26A. The signal, registered for the excitation with 1.26A laser diode current, is highly distorted. In this regime the amplifier was saturated, however it was necessary to investigate the signal in this range either because it was assumed that non-linear distortions should be seen in this regime for sure. In effect spectra of all signals are slightly different and it will be evidenced farther.

THE ANALYSIS OF EXPERIMENTAL RESULTS

As it was mentioned before, the successful elimination of the influence of PA signal linear distortion is the crucial condition of the analysis of thermal properties of investigated structure. One of the most effective methods for that is a use of contrasts characteristics in frequency domain [START_REF] Suszyński | Thermal Diagnostic of Power Thyristor[END_REF]. The contrast's definitions presented (1, 2) were evaluated basing upon the spectra of registered PA signal. In order of that it is necessary to determine the spectrum of photo-acoustic signal using DFT transformation [START_REF] Suszyński | The method of linear distortions elimination in photoacoustic investigation of high power thyristor structure[END_REF].

The comparison of the spectrum characteristic of the PA signal for different values of supply current is presented in Fig. 5. The changes of supply power results in the relocation of low frequency harmonic energy into the range of high frequencies. The changes of the spectra caused by non-linear distortions seem to be not so big but they may cause the substantial differences in the contrasts characteristics. At the same time, it can be observed that the spectra are unstable in time.

Zbigniew Suszyński, Radosław Duer, Mateusz Kosikowski

The method of non-linear distortions elimination in photoacoustic investigation... Fig. 4. The photoacoustic signal registered for the reference object.

©TIMA Editions

(1)

(2)

A) amplitude spectrum B) phase spectrum Fig. 5. Frequency spectra of registered photo-acoustic signal.

Zbigniew Suszyński, Radosław Duer, Mateusz Kosikowski

The method of non-linear distortions elimination in photoacoustic investigation...

©TIMA Editions/THERMINIC 2006

-page-ISBN: 2-916187-04-9

It seems to be caused by a relation between the level of non-linear distortions and low frequency acoustic disturbances registered by the microphone. Relatively high amplitude of low frequency disturbance influences on the range of all photo-acoustic signal changes and for this its spectrum is unstable with time. This value is in the denominator of eq.( 1) and for this a small fluctuation of its can result in a huge changes of amplitude contrast. This phenomenon is all the more significant that the signal-to-noise ratio nearby these frequencies is essentially smaller then for other harmonics. The level of disturbances is growing with frequency. It comes from the fact that the PA signal decreases with the growing frequency as 1/f. It means that signal-to-noise ratio decreases with frequency too. Simultaneously, the amplitude and the phase contrast characteristics change with power of excitation, which is caused by non-linear distortions of PA signal.

CONCLUSIONS

The results presented above lead to the conclusion that the analysis of thermal properties of the investigated structure is possible under two conditions:

1. Increase of the noise-to-signal ratio; 2. The elimination of non-linear distortions of the pa signal by performing the measurements in the range of exciation's power for which the non-linear distortions are not present.

  in ne ea ar r d di is st to or rt ti io on ns s Zbigniew Suszyński, Radosław Duer, Mateusz Kosikowski The method of non-linear distortions elimination in photoacoustic investigation... ©TIMA Editions/THERMINIC 2006 -page-ISBN: 2-916187-04-9

Fig. 2 . 1 .

 21 Fig. 2. A cross-section of the investigated structure and the reference object

Fig. 6 Fig. 6 .Fig. 7 .

 667 Fig.6presents the amplitude and phase contrasts characteristics obtained from (1) and (2). Both of contrast characteristics are strongly distorted especially for frequencies being multiple of 10Hz. It is caused by near to zero value of spectra for multiple of 10Hz.

Zbigniew Suszyński, Radosław Duer, Mateusz Kosikowski

The method of non-linear distortions elimination in photoacoustic investigation...

©TIMA Editions/THERMINIC 2006

-page-ISBN: 2-916187-04-9

In authors opinion, fulfilling these two opposite requirements is possible only with the use of special kind of excitation signal. Such a signal should be characterised by discrete and finite spectrum consisted of selected frequencies. Simultaneously, the amplitude of each harmonic should be distributed aversely to the frequency characteristics of components of the measurement setup. That way, the energy of excitation will be congregated in specific frequencies and the amplitude of signal will be distributed more efficiently. An example of the spectrum of desired excitation signal is presented in Fig. 7. Such modification of photoacoustic technique is cheap and seems to be effective. However it is necessary to point that the successful application of the broad-band excitation signal proposed is conditioned by the performance of the light source.


D A Benoy

THERMAL MODELING OF HIGH POWER LED MODULES

This paper presents a study of accuracy issues in thermal modeling of high power LED modules on system level. Both physical as well as numerical accuracy issues are addressed. Incorrect physical assumptions may result in seemingly correct, but erroneous results. It is therefore important to motivate the underlying key physical assumptions of a thermal model. In this paper thermal measurements are used to calibrate a computational fluid dynamics (CFD) model of a high power LED module model at a reference application condition, and to validate it at other application conditions.

INTRODUCTION

One of the key advantages of LED's are the high energy and optical system efficiencies and the product design freedom, due to their small form factor. Some illumination applications require white high power LED modules under a broad and versatile range of ambient boundary conditions. A prototype of a passively cooled high power density LED module is shown in picture A and B of Figure 1. For these applications thermal management is a major issue for both optical and reliability properties of the LED module shown in Figure 2. For thorough analysis of the thermal performance, and further optimization of these high LED modules a detailed thermal model has been developed for performing thermal simulations. For future design and product development work it is also important to know how to perform predictive thermal model calculations when prototype modules are not yet available. This implies that the accuracy of a thermal model much be such that it becomes predictive, or that the source of model inaccuracy can be identified. This paper presents a study of accuracy issues in thermal modeling of these high power LED modules. Thermal measurements are used to calibrate and validate a thermal computational fluid dynamics (CFD) model by comparing temperatures profiles of various components of the module. The IcePak ™ software package is used to implement a CFD thermal model on system level, which is shown in Figure 2. The CFD thermal comprises the whole LED module including the geometric details of the internal material interface layers, the heat sink geometry and environmental boundary conditions.

THERMAL MODELING

D.A. Benoy Thermal modeling of high power LED modules

The finite volume IcePak ™ software package is used to perform package level thermal model calculations of the LED module shown in Figure 2. It is well known that both physical and numerical issues have a large impact on the accuracy of the computed results [1]. Both issues are considered in the calibration and validation phases of the CFD thermal model of the high power LED module. Concerning physical issues CFD packages requires 1) the specification of the flow regime: i.e. laminar, or turbulent flow, 2) thermal radiation properties, and 3) specification of heat generation in heat dissipating components. Concerning the specification of the heat generation a part of the input electric power is converted into light, which is described by the wall plug efficiency (WPE) of the LED components. The WPE is defined as the ratio of the optical visible and electrical power and is an important characteristic performance number for LED devices. The heat dissipation in the LED's is then given by P heat = I×V×(1-WPE), where the WPE must be treated as an input parameter for the model calculations. The WPE is thus an additional source of uncertainty when it comes to a quantitative comparison between thermal CFD model calculations and experimental data.

Different types of model calculations are performed: Laminar and turbulent flow, In case of a turbulent flow regime turbulence is described by the zero-equation approach. Wall plug efficiencies of 10% and 15% are assumed and the electric input powers of the LED modules are 5.25 and 8.26W.

A combination of errors in the aforementioned physical issues may result in temperature distributions and temperature values at different module components, which are seemingly in good agreement with thermal measurements. For example the assumptions of a turbulent flow regime including thermal radiation with a "too low" WPE of 10% and an electric input power of 8.26W result in a temperature distribution and component temperatures at the die, the ceramic substrate, and the heat sink, which agree well with measured values. In case of IcePak ™ 's default grid settings the maximum temperature differences between measured and calculated temperatures is 4°C at the die (measured 169°C) and the heat sink (measured 124°C, see Table 1). Average surface temperatures of the high power LED module are measured with an IR-camera. The emissivities of the various different materials are calibrated at an elevated temperature of 40°C. But is the assumption of a turbulent flow justified? The nature of the natural convection flow is predicted by the dimensionless Rayleigh number Ra = Gr×Pr, where Gr and Pr are the Grashof and Prandtl numbers, respectively. The Grashof number is defined as According to [2] transition to turbulence in buoyant flows in vertical enclosures occurs at Ra ~ 4×10 6 (for enclosures with aspect ratios of the order 8). This critical value is much larger than the estimated Ra value (2.38×10 5 ) for the LED module so that the flow regime is expected to be laminar. This means that if we do include turbulence in the thermal calculations the heat transfer will be overestimated and the temperatures of the LED module will be too low.

Component

Switching from turbulent to laminar flow, and do the IcePak ™ calculation again then the maximum differences between measured and calculated temperatures increases up to 12°C at the die and the heat sink (see Table 2). These are systematic temperature differences of the order of 10%. The observation that the calculated temperatures at all module components are systematically too high brings us to the plausible conclusion that a WPE of 10% is too low. If the LED components are more efficient than anticipated, the WPE must then exceed 10%. Increasing the WPE from 10 to 15% will result in a decreased heat dissipation from 90%×8.26W = 7.46W to 85%×8.26W = 7.00W. Thermal model calculations with a laminar flow, thermal radiation, and WPE = 15% result in maximum temperature differences between measured and calculated temperatures of 4°C, and 7°C at the die and heat sink, respectively. With the satisfactory agreement between calculated and measured temperature values we have in fact calibrated the thermal model calculation by fitting the WPE to a value of 15% resulting in a heat dissipation of 7.00W. By comparing the temperature data measured at electric input power of 5.25W with calculated results using the same WPE (15%) resulting in a heat dissipation power of 85%×5.25W=4.46W the thermal CFD model is validated. This comparison yields maximum temperature differences of 6.7% at the hottest module parts (i.e. the die) and 6.4% at the coolest module part (the heat sink) so that the thermal CFD model for this high power LED module is indeed validated. The calculated LED chip and heat sink temperature are both overestimated (7.6% and 5.6%, respectively) at the validation condition of P=5.25W. For the submount there is a very good agreement between the calculated and measured temperatures.

Component

Component

Component

The thermal model calculations have been done with the default coarse mesh settings of IcePak TM . The resulting grid near the edges of the fins of the heat sink is shown in Figure 4A. It is well known that in terms of solution accuracy a denser mesh is more accurate than a coarse mesh. When the laminar flow case including thermal radiation with a heat dissipation of 7.00W, i.e. a WPE of 15%, is solved using a finer meshes near the heat sink boundaries as shown in Figure 4B, then there is almost a perfect match between the measured and the calculated temperatures at the die (hottest module part) and ceramic substrate. Calculated temperatures of the main components of the LED module using a fine mesh are given in the third column of Table 5. For these LED module parts the differences are of the order of 0.5%, while for the heat sink the difference is 4%.

Figure 1 .

 1 Figure 1. The pictures show the top and isometric view of a high power LED module. Some parts are painted black for IR-measurements.

Figure 2 .

 2 Figure 2. IcePak TM solid model of the high power LED module.

  β is the thermal expansion coefficient of air, g is the gravitational constant. ΔT is temperature difference between LED module and ambient temperature. For a safe estimation the largest possible ΔT is considered by taking a typical die temperature for the LED module temperature. L and ν are a characteristic module length and the kinematic viscosity coefficient, respectively. The Prandtl number is defined as Pr ν κ = where κ is the thermal diffusion coefficient. For the current application the approximate values for the Grashof and Prandtl numbers are

Figure 3 .

 3 Figure 3. Comparison between calculated and measured temperatures of the LED chips, substrate (top surface), and heat sink. The comparison is done for 2 powers (experimental loads of 5.25, and 8.26W).

Table 1 . Calculated versus measured temperatures at different module components. The flow is assumed to be turbulent, the WPE = 10%, and thermal radiation is included. For an electric power input of 8.26W this implies a power dissipation of 7.46W.

 1

		Calculated	Measured
		temperatur	temperature
		e
	LED chip	173.0°C	169.0°C
	Submount	131.6°C	134.1°C
	Ceramic substrate	127.6°C	129.0°C

Table 2 . Calculated versus measured temperatures for a laminar flow, including thermal radiation, and WPE = 10%.

 2

	Calculated	Measured

Table 3 . Calculated versus measured temperatures for a laminar flow, including thermal radiation, and WPE = 15%.

 3

		Calculated	Measured
		temperatur	temperature
		e
	LED chip	173.0°C	169.0°C
	Submount	134.4°C	134.1°C
	Ceramic substrate	131.6°C	129.0°C
	Heat sink	131.0°C	124.0°C

Table 4 . Validation of calculated temperatures for a laminar flow, including thermal radiation, and WPE = 15%. The power dissipation in LED module is 4.46W. Measured versus calculated temperatures in high power LED module

 4

		Calculated	Measured
		temperatur	temperature
		e
	LED chip	127.0°C	118.3°C
	Submount	101.9°C	101.3°C

©TIMA Editions/THERMINIC 2006-page-ISBN: 2-916187-04-9

D.A. Benoy Thermal modeling of high power LED modules

Wall plug efficiency

In the thermal model calculations it was assumed that the LED components have a WPE of 15%. The question arises whether this WPE is consistent with available data of WPE of the used LED components. It is important to note that we can only estimate the WPE of these components. From luminous efficacy measurements the WPE is determined to be ±10%, which is considerably lower than the WPE = 15% used in the thermal model calculations. The discrepancy between the optically estimated WPE, and the WPE used in the thermal model illustrates that the current thermal model for CFD calculations has a limited accuracy.

CONCLUSIONS

This study gives us confidence in how to perform thermal model calculations of future LED modules. Still an unknown, but important parameter in the thermal modeling of LED modules is the wall plug efficiency (WPE) of the LED components, which does not follow from the thermal CFD calculations. From the thermal point of view the WPE of the LED components can be fitted when measured temperature results are available. However, from optical measurements the WPE can be estimated and is lower than the thermally estimated WPE.

For the analysis and optimization of the thermal performance of LED modules the current thermal model is good enough.


Clément Cancès
email: clement.cances@cmi.univ-mrs.fr

Two-phase Flows Involving Discontinuities on the Capillary Pressure

Keywords: flows in porous media, finite volume method, discontinuous capillarity pressures, nonlinear PDE of parabolic type

We consider a immiscible two-phase flow through an one-dimensional heterogeneous porous medium made of an apposition of several homogeneous porous media. This leads to a nonlinear degenerate parabolic problem, with monotonous transmission conditions between the different homogeneous subdomains. We give an implicit finite volume scheme for such a twophase flow, and we prove the convergence of the inducted discrete solutions to a weak solution.

Under assumption on the initial condition, i.e. if it yields bounded flux, and if the total flowrate belongs to BV (0, T ), then the discrete solution obtained via the scheme converges toward a solution with bounded flux. We prove a L 1 -contraction principle for such a bounded flux solution. This uniqueness result is extended using a SOLA approach in the case where there are no particular regularity assumptions on the initial data, and we check that the discrete solution obtained via the implicit finite volume scheme converges to this SOLA.

Introduction: the model

We consider a heterogeneous porous medium Ω, which is an apposition of homogeneous porous media Ω i , representing the different geological layers, so that the physical properties of the medium only depend of the rock type and are piece-wise constant. We restrict our study to the one-dimensional case, even if all the results stated in section 2 can be quite easily adapted to larger dimensions (see [ENC 06]). For the sake of simplicity, we only deal with two geological layers with same size, because a generalisation to an arbitrary finite number of geological layers would only lead to notation difficulties. In the sequel, we denote by Ω =]-1, 1[ the heterogeneous porous medium, and by Ω 1 =] -1, 0[, Ω 2 =]0, 1[ the two homogeneous layers. The interface between the layers is thus {x = 0}. T is a positive integer.

The study of the flow in each Ω i , i = 1, 2 leads to the following degenerated parabolic equation:

φ i ∂ t u + ∂ x (qf i (u) -λ i (u)∂ x π i (u)) = 0, [1]
where, for i = 1, 2:

φ i is the porosity of Ω i , 0 < φ i < 1, -u(x, t) is the water-saturation, 0 ≤ u ≤ 1, -q(t) is the total flow rate, q ∈ L 1 (0, T ), q ≥ 0,

-f i is a continuous increasing function, f i (0) = 0, f i (1) = 1, -λ i is a continuous function, with λ i (0) = λ i (1) = 0, λ i (s) > 0 if 0 < s < 1, -the capillary pressure u → π i (u) ∈ C 1 ([0, 1], R) is increasing.
Focusing on the interface {x = 0}, we get the following transmission conditions: denoting u i the restriction of u to Ω i and

ϕ i (s) = s 0 λ i (a)π i (a)da, q(t)f 1 (u 1 )(0, t) -∂ x ϕ 1 (u 1 )(0, t) = q(t)f 2 (u 2 )(0, t) -∂ x ϕ 2 (u 2 )(0, t), [2] π1 (u 1 (0, t)) ∩ π2 (u 2 (0, t)) = ∅, [3] where πi (s) =    π i (s) if s ∈]0, 1[, ] -∞, π i (0)] if s = 0, [π i (1), +∞[ if s = 1.
We choose the boundary conditions:

q(t)f 1 (u 1 )(-1, t) -∂ x ϕ 1 (u 1 )(-1, t) = g(t) ∈ [0, q(t)],
[4]

∂ x ϕ 2 (u 2 )(1, t) = 0. [5]
and an initial data u 0 ∈ L ∞ (Ω).

Definition 1.1 (weak solution) A function u is said to be a weak solution if it fulfills:

1) u ∈ L ∞ (Ω × (0, T )), with 0 ≤ u ≤ 1, 2) for i = 1, 2, ϕ i (u) ∈ L 2 ((0, T ); H 1 (Ω i )), 3) for a.e. t ∈ (0, T ), π1 (u 1 (0, t)) ∩ π2 (u 2 (0, t)) = ∅, 4) for all ψ ∈ D(Ω × [0, T [), T 0 i=1,2 Ωi φ i u(x, t)∂ t ψ(x, t)dxdt + i=1,2 Ωi φ i u 0 (x)ψ(x, 0)dx + T 0 i=1,2 Ωi [q(t)f i (u)(x, t) -∂ x ϕ i (u)(x, t)] ∂ x ψ(x, t)dxdt + T 0 g(t)ψ(-1, t)dt - T 0 q(t)f 2 (u 2 )(1, t)ψ(1, t)dt = 0. [6]
2. Convergence of an implicit finite volume scheme

For the sake of simplicity, we choose to deal only with uniform discretization of Ω × (0, T ), which are defined by two positive integers M, N .

Discretization of Ω: we denote by δx = 1/N and

x j = j/N, ∀j ∈ [[-N, N]], x j+1/2 = j + 1/2 N , ∀j ∈ [[-N, N -1]].
Discretization of [0, T ]: we denote by δt = T /M , and

t n = nδt, for all n ∈ [[0, M]].
Discretization of the data:

-∀j ∈ [[-N, N -1]], u 0,D (x j+1/2 ) = u 0 j+1/2 = 1 δx xj+1 xj u 0 (x)dx, -q 0 = q 1 = 1 δt δt 0 q(t)dt, q n+1 = 1 δt t n+1 t n q(t)dt, -g 0 = g 1 = 1 δt δt 0 g(t)dt, g n+1 = 1 δt t n+1 t n g(t)dt, Equation [1] can be rewritten φ i ∂ t u + ∂ x F (x, t) = 0, with F (x, t) = q(t)f i (u) - ∂ x ϕ i (u).
We consider the following implicit scheme:

∀j ∈ [[-N, N -1]], ∀n ∈ [[0, M -1]], φ i u n+1 j+1/2 -u n j+1/2 δt δx + F n+1 j+1 -F n+1 j = 0, [7]
where the discrete fluxes are given by the following two-points scheme:

-

∀j ∈ J int = [[-N + 1, -2]] ∪ [[1, N -2]], ∀n ∈ [[0, M -1]], F n+1 j = q n+1 f i (u n+1 j-1/2 ) - ϕ i (u n+1 j+1/2 ) -ϕ i (u n+1 j-1/2 ) δx , [8] -∀n ∈ [[0, M -1]], F n+1 -N = g n+1 , [9] F n+1 N = q n+1 f 2 (u n+1 N -1/2 ). [10] -∀n ∈ [[0, M -1]], π1 (u n+1 0,1 ) ∩ π2 (u n+1 0,2 ) = ∅ [11] F n+1 0 = q n+1 f 1 (u n+1 -1/2 ) - 2(ϕ 1 (u n+1 0,1 ) -ϕ 1 (u n+1 -1/2 )) δx [12] = q n+1 f 2 (u n+1 0,2 ) - 2(ϕ 2 (u n+1 1/2 ) -ϕ 2 (u n+1 0,2 )) δx . [13]
The definition of F n+1 0 is justified by the following lemma, which is a consequence of the monotony of the transmission conditions.

Lemma 2.1 (discrete transmission conditions) Let (a, b) ∈ R 2 , let n ∈ [[0, M -1]],
then there exists a unique couple (c, d) ∈ [0, 1] 2 such that:

q n+1 f 1 (a) - 2(ϕ 1 (c) -ϕ 1 (a)) δx = q n+1 f 2 (d) - 2(ϕ 2 (b) -ϕ 2 (d)) δx , π1 (c) ∩ π2 (d) = ∅,
where for all s ≤ 0, f i (s) = 0, ϕ i (s) = 0, and for all s ≥ 1, We deduce from lemma 2.1 that [7] is a monotonous scheme, and thus we can apply a classical study (see e.g. [EYM 00, ENC 06]) to claim the following.

f i (s) = 1, ϕ i (s) = ϕ i (1).

Proposition 2.2 (discrete solution)

There exists a unique solution (u n+1 j+1/2 ) j,n to the scheme, fulfilling furthermore:

∀j ∈ [[-N, N -1]], ∀n ∈ [[0, M]], 0 ≤ u n j+1/2 ≤ 1.
We define the unique discrete solution u D almost everywhere in Ω × (0, T ) by:

u D (x, t) = u n+1 j+1/2 if (x, t) ∈]x j , x j+1 [×]t n , t n+1 ].
If u 0,D , v 0,D are two discrete initial data, and if u D , v D are the inducted discrete solution, we obtain:

∀t ∈]0, T ], i=1,2 Ωi φ i |u D (x, t) -v D (x, t)|dx ≤ i=1,2 Ωi φ i |u 0,D (x) -v 0,D (x)|dx. [14]
The discrete solution fulfils a discrete L 2 (0, T ; H 1 (Ω i )) estimate, i.e. there exists C(q, T, φ i , π i ) depending on the data, but not on the discretization of Ω × (0, T ) such that:

∀i = 1, 2, ∀(N, M ) ∈ (N ) 2 , T 0   j∈Jint δx ϕ i (u D )(x j+1/2 , t) -ϕ i (u D )(x j-1/2 , t) x j+1/2 -x j-1/2 2   dt ≤ C. [15]
We deduce from [15] some estimates on the space and time translates (see [EYM 00]), and so we obtain a compactness result on the set of the discrete solutions. Letting the discretization step tend to 0, one gets the following convergence result.

Proposition 2.3 (convergence to a weak solution) Let (M p ) p∈N , (N p ) p∈N be two sequences of positive integers tending to +∞, and let (u Dp ) p be the sequence of the inducted discrete solutions. Then there exists a weak solution u in the sense of definition 1.1 such that, up to a subsequence: ∀r ∈ [1, +∞[,

u Dp → u in L r (Ω × (0, T )).
All the results of this section still hold in larger dimensions. This can be proven by adapting the proof of [ENC 06] to transmission conditions [2],[3].

Uniform bound on the flux

As will be stated in section 4, we are not able to prove the uniqueness of the weak solution. Some additional regularity on the weak solution is necessary, so we define the notion of the bounded-flux solution.

Definition 3.1 (bounded-flux solution) A function u is said to be a bounded-flux solution if:

1) u is a weak solution in the sense of definition 1.1,

2) ∂ x ϕ i (u) belongs to L ∞ (Ω i × (0, T )).
In order to get an existence result for bounded-flux solutions, we need some more regularity on the data, as stated below.

Assumption 1 (additional regularity on the data) We assume that:

1) ∂ x ϕ i (u 0 ) ∈ L ∞ (Ω i ), 0 ≤ u 0 ≤ 1, 2) π1 (u 0,1 ) ∩ π2 (u 0,2 ) = ∅, where u 0,i is the trace of u 0|Ω i on {x = 0}, 3) q ∈ BV (0, T ), q ≥ 0, 4) g ∈ L ∞ (0, T ), 0 ≤ g ≤ q.
Let us define (F 0 j ) |j|≤N as in [8], [9], [10], [11], [12] and [13] where (n + 1) must be replaced by 0. Assumption 1 allows us to claim that:

max j=-N,...,N |F 0 j | ≤ q 0 + 2 max i=1,2 ∂ x ϕ i (u 0 ) L ∞ (Ωi) ≤ C. [16]
Lemma 3.1 (uniform bound on the discrete fluxes) Under assumption 1, there exists C > 0 only depending on u 0 , ϕ i , q, g and T such that:

max j∈[[-N,N]] max n∈[[0,M]] |F n j | ≤ C.

Sketch of the proof of lemma 3.1

We here perform an adaptation to the discrete case of the proof stated in [CAN ], which uses a maximum principle on the fluxes. For i = 1, 2, equation [1] can be rewritten

φ i ∂ t u + ∂ x F (x, t) = 0, [17] with F (x, t) = q(t)f i (u)(x, t) -∂ x ϕ i (u)(x, t). [18]
Formally, deriving [18] w.r.t. t and using [17] yields:

∂ t F + qf i (u) φ i ∂ x F -∂ x ϕ i (u) φ i ∂ x F = q f i (u).
[19]

This computations on the discrete case gives a monotonous scheme with (F n+1 j ) j,n as unknowns: we deduce from [7] and [8] that for all j ∈ J int , for all n ∈ [[0, M -1]], there exists (a n+1 j,j-1 , a n+1 j,j+1 ) ∈ (R + ) 2 such that:

F n+1 j-1 , F n+1 j , F n+1 j+1 •   -a n+1 j,j-1 1 + a n+1 j,j-1 + a n+1 j,j+1 -a n+1 j,j+1   ≤ F n j + q n+1 -q n + . [20]
Using strongly the monotony of the discrete transmission conditions [11], [12] and [13], we can also show that there exist (a n+1 0,-1 , a n+1 0,1 ) ∈ (R + ) 2 such that [20] also holds for j = 0. Since the discrete boundary conditions [9] and [10] insure that F n+1 N and F n+1 -N stay uniformly bounded by q ∞ , we can claim that

max j F n+1 j ≤ max q ∞ , max j F n j + q n+1 -q n + ,
and thus

max n max j F n+1 j ≤ max q ∞ , max j F 0 j + 2T V (q).
Thanks to assumptions 1 and estimate [16], the right-hand-side member is uniformly bounded. We conclude the proof of lemma 3.1 by showing the corresponding underestimate:

min n min j F n+1 j ≥ min 0, min j F 0 j -2T V (q).
Proposition 3.2 (convergence to a bounded-flux solution) Let (M p ) p , (N p ) p be two sequences of positive integers tending to +∞, and let (u Dp ) p be the inducted sequence of discrete solutions. Under assumption 1 there exists a bounded-flux solution u in the sense of defintion 1.1 such that, up to a subsequence: ∀r ∈ [1, +∞[,

u Dp → u in L r (Ω × (0, T )).

Uniqueness of the bounded-flux solution

Theorem 4.1 (L 1 -contraction for bounded-flux solution) Suppose that for i = 1, 2,

f i • ϕ -1 i ∈ C 0, 1 2 ([0, ϕ i (1)], R).
Let u,v be two bounded-flux solutions in the sense of definition 3.1 associated to initial data u 0 ,v 0 . Then u and v belong to C([0, T ]; L 1 (Ω)), and for all t ≥ 0,

i=1,2 Ωi φ i (u(x, t) -v(x, t)) ± dx ≤ i=1,2 Ωi φ i (u 0 (x) -v 0 (x)) ± dx.
Corollary 4.2 Keeping the notations of proposition 3.2, and under assumptions 1, the whole sequence (u Dp ) p converges in L r (Ω × (0, T )) towards the unique bounded flux solution.

Sketch of the proof of theorem 4.1

The proof is based on a doubling variable method w.r.t. the time variable, first in each Ω i , and on a transmission lemma 4.3 to deal with the interface {x = 0}.

Since u, v are weak solutions, they belong to C([0, T ]; L 1 (Ω)) (see lemma 7.41 in [M`96]). Some classical computations (see e.g. [OTT 96]) give: ∀θ ∈ BV (0, T ), θ ≥ 0, θ(T ) = 0, ∀ψ ∈ W 1,1 (Ω), ψ ≥ 0, ψ(0) = 0,

T 0 i Ωi φ i (u(x, t) -v(x, t)) + ∂ t θ(t)ψ(x)dxdt + i Ωi φ i (u 0 (x) -v 0 (x)) + ψ(x)θ(0)dx + T 0 θ(t) i Ωi q(t)(f i (u)(x, t) -f i (v)(x, t)) + ∂ x ψ(x)dxdt - T 0 θ(t) i Ωi ∂ x (ϕ i (u)(x, t) -ϕ i (v)(x, t)) + ∂ x ψ(x)dxdt ≥ 0.
[21]

We need the following key lemma, which is a consequence of the uniform bound on the flux and of the monotony of the transmission conditions [2] and [3].

Lemma 4.3 If ε > 0, we denote ξ ε (x) = max 0, 1 -|x| ε . ∀θ ∈ BV (0, T ), θ ≥ 0, θ(T ) = 0, lim inf ε→0 T 0 θ(t) i=1,2 Ωi q(t)(f i (u)(x, t) -f i (v)(x, t)) + -∂ x (ϕ i (u)(x, t) -ϕ i (v)(x, t)) + ∂ x ξ ε (x)dxdt ≥ 0. Let t ∈]0, T [. Taking ψ(x) = (1 -ξ ε (x)), θ t (s) = χ [0,t] (s) in [21]
, and letting ε tend to 0 gives:

i=1,2 Ωi φ i (u(x, t) -v(x, t)) + dx ≤ i=1,2 Ωi φ i (u 0 (x) -v 0 (x)) + dx.
Similar calculations leads to the same inequality with (•) -instead of (•) + .

Convergence to the unique SOLA

The goal of this part is to get the convergence of the whole sequence of discrete solutions despite u 0 "only" belongs to L ∞ (Ω).

Definition 5.1 (SOLA to the problem) A function u is said to be a SOLA if it fulfils:

u is a weak solution in the sense of definition 1.1, -there exists a sequence (u n ) n∈N of bounded flux solutions such that

u n → u in C([0, T ]; L 1 (Ω)), as n → +∞.
Theorem 5.1 (Convergence to the unique SOLA) Let u 0 ∈ L ∞ (Ω), 0 ≤ u 0 ≤ 1, and q ∈ BV (0, T ), g ∈ L ∞ (0, T ), 0 ≤ g ≤ q. There exists a unique SOLA u in the sense of definition 5.1. Let (M p ) p , (N p ) p be two sequences of positive integers tending to +∞, and let (u Dp ) p be the sequence of the inducted discrete solutions. Then the whole sequence (u Dp ) p converges to the unique SOLA u in L r (Ω × (0, T )), r ∈ [1, +∞[ as p tends to +∞.

Sketch of the proof of theorem 5.1

Let (u 0,ν ) ν∈N be a sequence of approximate initial data such that: for all ν ∈ N, (u 0,ν , q, g) fulfils assumptions 1, and lim ν→+∞ u 0,ν -u 0 L 1 (Ω) = 0. Let (u ν ) ν be the inducted sequence of bounded flux solutions. The L 1 -contraction principle stated in theorem 4.1 allows us to claim that (u ν ) ν is a Cauchy sequence in C([0, T ]; L 1 (Ω)). It is then easy to check that u is a weak solution, so it is a SOLA. The L 1 -contraction principle still holds for the SOLAs, which is thus unique.

Let (M p ) p , (N p ) p be two sequences of positive integers tending to +∞, (u Dp ) p (resp.(u ν,Dp ) p ) be the sequence of the inducted discrete solutions for u 0 (resp. u 0,ν ) as initial data. From [14], we obtain: Letting ν tend to +∞ shows that u is a SOLA in the sense of definition 5.1.

  Furthermore, (a, b) → c and (a, b) → d are continuous and non-decreasing w.r.t. each one of their arguments.

i=1, 2

 2 Ωi

φ

  i |u Dp (x, t) -u ν,Dp (x, t)|dx ≤ i=1,2 Ωi φ i |u 0,Dp (x) -u 0,ν,Dp (x)|dx.Letting p tend to +∞, u Dp (x, t) tends to a weak solution u thanks to proposition 2.3. u ν,Dp tends to a bounded flux solution u ν thanks to proposition 3.2, and so:i=1,2 Ωi φ i |u(x, t) -u ν (x, t)|dx ≤ i=1,2 Ωiφ i |u 0 (x) -u 0,ν (x)|dx.


Antoine Dozolme

Hossam Metwally

Thierry Marchal

ELECTRONICS COOLING FAN NOISE PREDICTION

Using the finite volume CFD software FLUENT, one fan was studied at a given flow rate (1.5m 3 /min) for three different operational rotating speeds: 2,000, 2,350 and 2,700 rpm. The turbulent air flow analysis predicts the acoustic behavior of the fan. The best fan operating window, i.e. the one giving the best ratio between noise emissions and cooling performance, can then be determined. The broadband noise acoustic model is used. As the computation is steady state, a simple Multiple Reference Frame model (MRF, also known as stationary rotor approach) is used to represent the fan. This approach is able to capture the effects of the flow non-uniformity at the fan inlet together with their impact on the fan performance. Furthermore, it is not requiring a fan curve as an input to the model. When compared to the available catalog data the simulation results show promising qualitative agreement that may be used for fan design and selection purposes.

INTRODUCTION

In order to evacuate the quickly increasing heat generated by the always more powerful computers, new fans are constantly designed. The performance of these new designs in term of thermal evacuation is obvious. However, the quality of the fan is not restrained to its capability to evacuate heat through intense air flow anymore. The competitiveness of the end product, such as a laptop, now includes the technical performances of the computer but also the comfort it is maintaining to the endusers. Hence, the best designed fan has little value if it is too noisy. More and more fan designers are now including this noise dimension in the proper design criteria Fans acoustics emissions are the main noise source in forced convection air cooled electronics systems. Improving the fan design for a given system will not only result in an increase in the amount of air flow it delivers, but may also lead to a reduction of noise emissions. A fan operating at a higher efficiency delivers more air for the same amount of torque; it is also less noisy within the same system. Fan noise is becoming an important issue especially with the increasing usage, and power, of portable electronics.

CFD analyses can help design the fan to improve both hydrodynamic and acoustical performance. In this work, steady state analysis of an axial fan is carried out at three different rpm's to estimate the broad band noise the fan emits. The surface acoustical power level is compared to available catalog data.

CFD MODELING

The GAMBIT package is used to generate the geometry and the mesh of the model. The flow is analyzed by solving the mass and momentum conservation equations. Momentum equation:

( ) ( ) F g p v v v t r r r r r + + -∇ = ∇ + ∂ ∂ ρ ρ ρ . ( 1
)
Continuity equation:

( ) m S v t = ∇ + ∂ ∂ r ρ ρ . ( 2
)
For a 3D model, the momentum equation is discretized as follow:

( )

V S A A v faces faces N f N f f n f f f f φ φ φ φ ρ + ∇ Γ = ∑ ∑ r r r . . (3)
In this analysis, the model is isothermal. Turbulence is taken into account through the standard k-epsilon model. It is mathematically described by the following equations:

( ) ( ) k M b k j k t j i i S Y G G x k x ku x k t + - - + +         ∂ ∂         + ∂ ∂ = ∂ ∂ + ∂ ∂ ρε σ µ µ ρ ρ (4) ( ) ( ) ( ) ε ε ε ε ε ε ε ρ ε ε σ µ µ ρε ρε S R k C G C G k C xj x u x t b k t j i i + - - + +         ∂ ∂         + ∂ ∂ = ∂ ∂ + ∂ ∂ 2 2 3 1

2

(5) Where G k represents the generation of turbulence kinetic energy due to the mean velocity gradients, G b is the generation of turbulence kinetic energy due to buoyancy. Y M represents the contribution of the fluctuating dilatation in compressible turbulence to the overall dissipation rate.

©TIMA Editions/THERMINIC 2006

-page-ISBN: 2-916187-04-9

The quantities σ k and σ ε are the inverse effective Prandtl numbers for k and ε, respectively. S k and S ε are userdefined source terms

The geometry consists of a tube. In the middle of the tube, a fan is explicitly modeled (figure 1). The mesh is made of 1,500,000 hexahedral elements. The inlet is represented in blue and the outlet is red. The walls are presented in light yellow with a mesh representation on them. The mesh is refined close to the fan in order to capture the flow field around it. The model is assumed isothermal, and the gravity is not taken into account. The inlet is defined as a constant velocity inlet. The air enters the domain at 2.5m/s, corresponding to the value of 1.5 m 3 /min. The outlet is specified as a free boundary condition. The entire domain walls except the fan walls have a zero velocity boundary condition. The fan is modeled using the MRF approach for each of the three velocities.

The standard k-epsilon model is used to describe the turbulence behavior. The value of k and epsilon will be used in the acoustics models to perform broadband acoustic calculation.

The air passing inside the fan has a density of 1.225kg/m 3 and a viscosity of 1.78940.10 -5 Pa.s. The air compressibility is neglected.

BROADBAND MODEL

In many practical applications involving turbulent flows, noise does not have any distinct tones, and the sound energy is continuously distributed over a broad range of frequencies. In those situations involving broadband noise, statistical turbulence quantities readily computable from RANS equations (mean flow field, turbulent kinetic energy (k) and the dissipation rate (ε)) can be utilized, in conjunction with semi-empirical correlations and acoustic analogies, to shed some light on the source of broadband noise. Unlike the computationally expensive direct method (which involves computationally expensive detailed unsteady simulation) and the integral method, the broadband-noise-source models do not require transient solutions.

The broadband model is however limited to the broadband noise characteristics prediction and does not provide any tonal performance data

Proudman [START_REF] Proudman | The Generation of Noise by Isotropic Turbulence[END_REF] ,using Lighthill's acoustic analogy, derived a formula for acoustic power generated by isotropic turbulence without mean flow. More recently, Lilley [START_REF] Lilley | The radiated noise from isotropic turbulence revisited[END_REF] re-derived the formula by accounting for the retarded time difference which was neglected in Proudman's original derivation. Both derivations yield acoustic power due to unit volume of isotropic turbulence.

The Proudman's formula gives an approximate measure of the local contribution to total acoustic power per unit volume in a given turbulence field. Proper caution, however, should be taken when interpreting the results in view of the assumptions made in the derivation, such as high Reynolds number, small Mach number, isotropy of turbulence, and zero mean motion.

MRF FORMULATION

The MRF model gives FLUENT the ability to model flows in an accelerating reference frame. In this situation, the acceleration of the coordinate system is included in the equations of motion describing the flow.

When problems such as axial fan are defined in a rotating reference frame, the rotating boundaries become stationary relative to the rotating frame, since they are moving at the same speed as the reference frame.

When the equations of motion are solved in a rotating frame of reference, the acceleration of the fluid is augmented by additional term that appears in the momentum equation. The left hand side of the momentum equation appears as follow for an inertial reference frame:

( ) ( ) v v v t r r r ρ ρ . ∇ + ∂ ∂ (6) ©TIMA Editions/THERMINIC 2006 -page- ISBN: 2-916187-04-9
For flow in rotating domains, the equation for conservation of mass, or continuity equation can be written as follows:

( ) m r S v t = ∇ + ∂ ∂ r ρ ρ . (7)

RESULTS

The computational results are compared with available data from the manufacturer. The comparison is done for the static pressure rise and the noise prediction. The experimental data are obtained from [4].

Each computational run takes approximatively 3 hours on a 2 processors Pentium 4 machine running Linux.

Pressure rise is taken from a point in front of the fan to a point right after the fan. The computed results are less than 3Pa different with measurements what represents an average difference of about 10%. The faster the fan gets, the bigger the pressure rise. The next step is to compare the noise level. The results are qualitatively comparable. The noise level is increasing when the rpm are increasing in both cases.

The faster the fan gets, the noisier it is, as expected.

The noise distribution within the fan is then studied. For all three fans, the noise repartition is analyzed on the fan itself. Detailed analyses are described in section 5.2 to 5.4

Flow field study.

The air enters the domain at a constant velocity (2.5m/s). The flow direction is on the Z axis, from the lower Z values to the higher ones. When entering the fan, the air gets a swirl velocity from the fan, accelerating it (green to red on figure 2). After the fan, it clearly appears that the highest velocities remain on the larger radius of the cylindrical domain whereas the lowest velocities are found with the lower radius. On the blade, the larger velocities are found on the outside of the blades. The pathlines representation illustrates the swirling flow induced by the fan. The air is going straight to the fan, after the device, the air is rotating until it gets out of the domain. Inside the fan, the air is rotating. Figure 4 shows the air entering the fan with a quasi normal velocity, when the air is leaving the fan, it has been given a rotational movement.

©TIMA Editions/THERMINIC 2006 -page-ISBN: 2-916187-04-9 The pressure distribution on the fan for all three cases has the same repartition (figure 5, figure 6, figure 7). The pressure is increasing inside the fan. The lower pressures are located at the inlet of the fan whereas the higher pressures are located after the fan. At 2,000 rpm (figure 5), the pressure drop is 15.3Pa.

The 2,000 rpm case

The noise is the most important at the fan inlet. The noise is nearly non existent after the fan, and decreases inside the fan. At the leading edge of the blades, the noise is the largest. The sound level repartition follows the same contours as the pressure, but it is the larger where the pressure is the lower. The same behaviors as for the 2,000 rpm case are observed. There is a slight increase of the noise level and pressure on the blades. The noise repartition remains the same, but at a higher level.

The 2,350 rpm case

Most of the noise is observed on the negative static pressure zones on the blades. Interestingly, the zones with no noise on the fan hub (blue part) are occurring where the pressure is zero.

The 2,700 rpm case

The noise repartition remains the same as for the 2,350 rpm case on the fan body. The slight difference is a decrease of the silent part in front of the fan. The global noise level has slightly increased as well. Increasing the fan rotating speed increases the noise level for a given mass flow rate.

CONCLUSIONS

FLUENT can predict the noise repartition inside a fan for several flow regime. The comparisons with experimental data for both global noise level and static pressure drop are comparable with the computed results. The broadband model gives accurate results for this study , for a limited CPU time, due to its physical formulation (steady state, standard k-epsilon turbulence model).

Fig 1 :

 1 Fig 1: Geometry and mesh.

Fig 2 :

 2 Fig 2: velocity magnitude (m/s) repartition in the computational domain for 2,700 rpm.

Fig 3 :

 3 Fig 3: Path line inside the domain for 2,350 rpm. Path lines are colored by ID.

Fig 4 :

 4 Fig 4: Velocity vector representation inside the fan for the 2,000 rpm case.

Fig 5 :

 5 Fig 5: Static pressure (Pa) repartition on the fan at 2,000 rpm At 2,350 rpm, negative static pressure locations are larger at the leading edges (figure 6). The maximum pressure has increased as well. The overall pressure drop as increased as well, up to 23 Pa.

Fig 6 :

 6 Fig 6: Static pressure (Pa) repartition on the fan at 2,350 rpm At 2,700 rpm (figure 7), the negative static pressure zones have increased, as well as the maximum pressure on the blades. The overall pressure drop is now 33.7 Pa.

Fig 8 :

 8 Fig 8: acoustic level (dB) repartition on the fan at 2,000 rpm

Fig 9 :

 9 Fig 9: acoustic level (dB) repartition on the fan at 2,350 rpm

Fig 10 :

 10 Fig 10: acoustic level (dB) repartition on the fan at 2,700 rpm

Table 1

 1

		shows the comparison of
	the results:
	rpm	Experimental	Computed
	2,000	12.5	15.3
	2,350	21.1	23
	2,700	35.6	33.7

Table 1 :

 1 Static pressure rise (Pa)

Table 2 shows

 2 the results.

	rpm	Experimental	Computed
	2,000	38	48
	2,350	42	55
	2,700	47	62
	Table 2: surface acoustic power level on the fan surface
		(dB)