Victor Birman

EFFECT OF TEMPERATURE ON VIBRATIONS OF PHYSICALLY NONLINEAR PIEZOELECTRIC RODS

The paper is concerned with a physically nonlinear behavior of piezoelectric rods. Such rods are often employed as transducers. We evaluate the effect of a uniformly distributed temperature on vibrations of a piezoelectric rod polarized in the axial direction and subjected to a harmonic (periodic-in-time) electric field. We conclude that temperature can significantly affect vibrations of a physically nonlinear transducer in a broad spectrum of frequencies. It is also shown that physical nonlinearity has a strong effect on the dynamic response of piezoelectric rods to alternating electric fields, especially in the vicinity of the resonance frequency and at high power fields.

INTRODUCTION

Physically linear relationships between the tensors of stress and strain, on one hand, and the vectors of electric field and electric displacement, on the other, can be employed only if the electric field is weak and, as a result of that, the nonlinear effects in piezoelectric materials may not be considered. While geometric nonlinearity can be nevertheless partially incorporated into the linear equations by an appropriate choice of the strain-displacement relationships, neither the geometric nonlinearity, nor physical nonlinearity of the piezoelectric material could be consistently accounted for in the linear constitutive equations.

Constitutive equations containing nonlinear products of the components of the strain tensor (geometrically nonlinear effect) and physically nonlinear terms have been obtained by Maugin et al. [START_REF] Maugin | Nonlinear Electromechanical Couplings[END_REF]. However, the complexity of these equations, as well as the lack of experimental data for the coefficients in the nonlinear terms, prevented wide acceptance of nonlinear constitutive equations in practical design efforts and in practical applications. This circumstance makes it important to elucidate a relative contribution of geometrically and physically nonlinear terms to the dynamic response of piezoelectric materials subject to relatively strong electric fields. There is an obvious incentive to assess the need in considering nonlinear terms in the dynamic analysis of structural elements made of such materials. While the comprehensive review of studies of the effect of physical nonlinearity is outside the scope of this paper, the reader is referred to the forthcoming papers [START_REF] Birman | Effect of Physical Nonlinearity on Vibrations of Piezoelectric Rods[END_REF], [START_REF] Birman | Effect of Material Nonlinearity on the Mechanical Response of Some Piezo-electric and Photonic Systems[END_REF] for a reference on this subject.

The problem of the mechanical behavior of piezoelectric transducers driven by voltage, which is applied along the direction of the wave propagation, has been considered since the mid-sixties (see, for example, the monograph by Mason,[START_REF] Mason | Crystal Physics of Interaction Processes[END_REF]). Physically nonlinear effects in piezoelectric transducers were considered in recent papers of Tan and Tong [START_REF] Tan | A Onedimensional Model for Non-linear Behavior of Piezoelectric Composite Materials[END_REF], [START_REF] Tan | Micromechanics Models for Non-linear Behavior of Piezoelectric Fiber Reinforced Composite Materials[END_REF], Wagner [START_REF] Wagner | Nonlinear Longitudinal Vibrations of Piezoceramics Excited by Weak Electric Fields[END_REF], [START_REF] Wagner | Non-linear Longitudinal Vibrations of Non-slender Piezoelectric Rods[END_REF] and Wagner and Hagedorn [START_REF] Wagner | Nonlinear Effects of Piezoceramics Excited by Weak Electric Fields[END_REF]. The recent papers by the authors [START_REF] Birman | Effect of Physical Nonlinearity on Vibrations of Piezoelectric Rods[END_REF], [START_REF] Birman | Effect of Material Nonlinearity on the Mechanical Response of Some Piezo-electric and Photonic Systems[END_REF] also address this problem. However, thermal effects were not included in these formulations.

The present paper contains a comprehensive formulation of equations of motion of a piezoelectric rod driven by a dynamic voltage in the direction of the wave propagation (axial direction) as shown in Fig. 1. The emphasis is on the assessment of the effect of temperature.

ANALYSIS

The energy density is obtained without considering the geometric nonlinearities as an extension of the density expression obtained by Maugin et al. [START_REF] Maugin | Nonlinear Electromechanical Couplings[END_REF] for the case when the effect of temperature was not accounted for: are other coefficients accounting for the interaction of the tensor of strains, vector of electric field and temperature. Note that the terms, which are dependent on temperature only, are excluded from (1), since these terms do not affect the dynamic problem under consideration. Additional higher-order terms could be included in [START_REF] Maugin | Nonlinear Electromechanical Couplings[END_REF], but this is not advisable, because their effect is relatively small and the information about some of the coefficients in (1) is not available, even for this "truncated" formulation.

T E T E E T E T T E p T T E E l E E E E E E e E e C a m ma n m mn m m m m n m mn p n m mnp n m mn m m m m ε κ χ η ε ε ρ ε δ ε λ ε ε ε ε ε ε ε ε β α αβ α α α α α α β α αβ α α β α αβ - - - - - - - - - - - - = Ψ 2 1 2 1 2 1 6 1 2 1 2 1 2 1 2 2 (1)
The constitutive relationships for the stresses and electric displacements are obtained as

m m E D ∂ Ψ ∂ - = ∂ Ψ ∂ = , α α ε σ (2)
This results in the following expressions for the components of the stress tensor and the electric displacements:

T E T T T E E l E e E e C m m n m mn m m m m α β αβ α α α β αβ α β αβ α κ ε ρ δ λ ε ε σ - - - - - - - = 2 2 1 (3) (4)
The terms underlined in the right side of these equations are retained in the linear formulation. As indicated above, numerous authors implicitly account for the nonlinearity by using variable coefficients of these terms, dependent on stress, electric field or temperature (see discussion in [START_REF] Birman | Effect of Physical Nonlinearity on Vibrations of Piezoelectric Rods[END_REF]). It is observed that by retaining linear terms in strain, electric field and temperature one could reduce the eqns. [START_REF] Birman | Effect of Material Nonlinearity on the Mechanical Response of Some Piezo-electric and Photonic Systems[END_REF][START_REF] Mason | Crystal Physics of Interaction Processes[END_REF] to the equations derived using the thermodynamic electric Gibbs function [START_REF] Mason | Crystal Physics of Interaction Processes[END_REF].

Notably, the terms quadratic in strain and/or electric field components and multiplied by temperature, i.e.,

T E E T E T n m m , , α β α ε ε ε
, are excluded from ( 3) and ( 4) as a result of eliminating higher-order temperature dependent terms in [START_REF] Maugin | Nonlinear Electromechanical Couplings[END_REF].

Eqns. [START_REF] Birman | Effect of Material Nonlinearity on the Mechanical Response of Some Piezo-electric and Photonic Systems[END_REF][START_REF] Mason | Crystal Physics of Interaction Processes[END_REF] are further simplified reflecting the assumption that temperature is independent of time. Accordingly, retaining dynamic terms only yields

T E T E E l E e E e C m m n m mn m m m m α β αβ α β αβ α β αβ α κ ε ρ ε ε σ - - - - - = 2 1
(5)

T T E E l E E E e e D m m mn n mn p n mnp n mn m m m α α α α β α αβ α α ε κ χ ε ε ε ε ε ε + + + + + + = 2 1 2 1 ( 6
)
In the subsequent analysis of forced vibrations generated by an alternating electric field we use the expression (5) for the stress. It is noteworthy that "material constants" in (5) are functions of the electric field and temperature:

E T mn mn T m m T m m T E l l e e e e C C αβ αβ α α αβ αβ α α αβ αβ ρ ρ = = = = = , , , , , (7)
Here the subscripts indicate that the corresponding term was measured at a constant value of the vector of the electric field E or temperatureT . Equations ( 5) and ( 6) represent a weak physically nonlinear formulation for a piezoelectric material, where the higher-order terms discussed above were not considered.

In the problem addressed in this paper, we are concerned with a dynamic response of a piezoelectric rod to the electric field applied in the axial direction, i.e.

( ) (

)

t E t E z = 3 .
Accordingly, the nonzero elements of the tensor of stress given by ( 5) are ( 8)

T T E T T p E l E E E e e D m n mn m m n mn p n mnp n mn m m m α α α α β α αβ α α ε κ χ η ε ε ε ε ε ε + + + + + + + + = 2 2
          -           - - - -                     - - - - - - - - - - - - - - - - - - =           κ κ κ ε ε ε ρ ρ ρ ρ ρ ρ ρ ρ ρ σ σ σ θ θ
The equation of motion of a polarized piezoelectric rod with arbitrary boundary conditions driven by an electric field in the axial direction can be derived from the Hamilton principle [START_REF] Birman | Effect of Physical Nonlinearity on Vibrations of Piezoelectric Rods[END_REF]. In this paper, such an equation was obtained from the generalized Galerkin method that enables one to account for the static (stress) boundary conditions. In the present problem it is assumed that the rod is clamped at one end (z=0), while its other end is free. In such a case, the generalized Galerkin procedure implies that

( ) ( ) ( ) ( ) 0 , 0 = = - - ∫ h Z h z dz z Z w m i z h i z z σ σ & & (9)
In ( 9), h is the height of the rod and ( ) z Z i are generalized coordinates that are used to characterize axial vibrations in terms of series that represent a sum of products of functions of time and coordinate:

( ) ( ) ( ) ∑ = i i i t T z Z t z w , (10)
The axial stress that appears in ( 9) can be derived from [START_REF] Wagner | Non-linear Longitudinal Vibrations of Non-slender Piezoelectric Rods[END_REF]. Particularly, in a onedimensional problem, the expression for this stress is ( ) ( )

2 3 333 3 33 33 33 3 333 33 2 1 E l E T e T E e C z z - - - - - = κ ε ρ σ (11)
where

z z w, = ε .
The axial displacements are chosen in the form of series satisfying the corresponding boundary conditions. However, it appears impossible to choose such series, while satisfying the conditions both at the "anchored" cross section 0 = z (where the axial displacement is equal to zero) and at the free end of the rod h z = (where the axial stress is equal to zero). Therefore, the solution is sought in the form that satisfies only the kinematic condition at the clamped end (z=0), while the boundary condition at z=h is accounted for through the generalized Galerkin formulation (9):

( ) h z i h t W w i i 2 sin π ∑ = (12)
The substitution of ( 11) and ( 12) into [START_REF] Wagner | Nonlinear Effects of Piezoceramics Excited by Weak Electric Fields[END_REF], where

h z i Z i 2 sin π =
, leads to a system of uncoupled differential equations. In particular, the i-th equation of this system is ( )

2 sin 2 1 ) ( 8 2 2 333 33 33 2 2 33 333 33 2 π κ π ρ i E l E T e W i T E e C W mh z z i z i       + - = - - + & & (13)
It is assumed that the rod is driven by a periodic harmonic electric field, ( )

t
π ω ω κ π ω ρ i t E l t E T e E l W i t E e T C W mh i i       + - + = - - + & & (14)
The coefficient of thermal expansion does not enter the equation of motion, since the thermal expansion does not affect the motion of the rod, which is unconstrained in the axial direction in a geometrically linear formulation.

The coefficients in the constitutive equations and in the equation of motion, i.e. and 333 l , depend on temperature. In addition, the explicit effect of temperature is present through the terms proportional toT . Accordingly, available information on the effects of temperature on the stiffness and "material constants" based on a comparison of empirical data with linear constitutive equations reflects two basic phenomena: (a) effect of temperature on the material constants, and (b) terms explicitly dependent on temperature that are not considered in a linear formulation.

The present formulation can be extended to the problem of axisymmetric axial-radial vibrations. In this case, extending the generalized Galerkin formulation implies [START_REF] Birman | Effect of Physical Nonlinearity on Vibrations of Piezoelectric Rods[END_REF], [START_REF] Birman | Effect of Material Nonlinearity on the Mechanical Response of Some Piezo-electric and Photonic Systems[END_REF]:

( ) 0 2 sin 2 sin 2 0 0 0 2 0 0 = = -       - - + ∂ ∂ ∫ ∫ ∫ ∫ ∫ θ π σ θ π σ σ σ π π θ dzd i a a r drdzd a r i r u m r r h a h r r r & & ( ) 0 2 sin 2 sin 2 0 0 0 2 0 0 = = -       - ∂ ∂ ∫ ∫ ∫ ∫ ∫ θ π σ θ π σ π π drd i r h z drdzd h z i r w m z h a a z z & & ( 15
)
Here a is the radius of the rod, the axial motion is given by (12) and the radial motion is represented by

( ) ∑ = = 1 2 sin i i a r i t U u π (16)
The substitution of the stresses (8) and the expression for the harmonic electric field in the axial direction as that considered above yield the system of two coupled nonhomogeneous Mathieu equations for each value of an integer i that is omitted here for brevity.

NUMERICAL EXAMPLES

The numerical analysis is confined to the case of axial vibrations, i.e., the effect of radial motion is disregarded. Accordingly, the analysis relies on the solution of the nonhomogeneous Mathieu equation ( 14). This solution is discussed in detail in the paper by Birman [START_REF] Birman | Effect of Physical Nonlinearity on Vibrations of Piezoelectric Rods[END_REF]. The height of the rod does not affect the following results. The amplitude of the electric field was chosen to reflect reported highest allowable fields that can be applied to currently used piezoelectric materials (less or equal to 2MV/m).

The following examples illustrate the effect of temperature on two representative materials, i.e. PZT-5H and PZN-4.5%PT. The properties of these materials at room temperature are shown in Table 1 [START_REF] Birman | Effect of Physical Nonlinearity on Vibrations of Piezoelectric Rods[END_REF]. The information on the effect of temperature on material properties of piezoelectrics is limited, but it is known that this effect on the stiffness within the range of operating temperatures is small. Therefore, in the subsequent examples it is assumed that (

)

33 33 33 C T C ≈ -ρ
.

The effect of temperature on the coefficient 333 l is not known, but numerical analysis [START_REF] Birman | Effect of Physical Nonlinearity on Vibrations of Piezoelectric Rods[END_REF] showed that the effect of this coefficient is negligible. Although data on the effect of temperature on the piezoelectric coefficient 33 e are not available, it is known that the piezoelectric constant 33 d increases with temperature [START_REF] Mukherjee | Nonlinear Properties of Piezoelectric Ceramics[END_REF], [11]. Therefore, the former constant should also increase with temperature reflecting the relationship [START_REF] Mason | Crystal Physics of Interaction Processes[END_REF] where the stiffness is measured for a constant electric field. As reported, the variations of the coefficient 33 d with temperature can be as high as hundreds percent, even if temperature varies between zero and 100 o C (see for example, [11]). Accordingly, if the stiffness remains little affected by such variations of temperature, the coefficient T e 33 33 κ may significantly increase within such limited range of temperature variations. This is reflected in the following examples (a direct correlation between temperature and the value of the above-mentioned coefficient is impossible due to the lack of experimental data).

Numerical results were obtained for a physically "reasonable" range of variations of the temperature adjusted piezoelectric coefficient T e e 33 33 33 κ -= that accounts for both the effect of temperature on the piezoelectric constant, as well as the presence of the term explicitly dependent on temperature. The examples focus on the analysis of the effect of the adjusted piezoelectric coefficient 33 e on the response of the rod, assuming that other thermal effects, while explicitly shown in the theoretical solution, are less important. The effect of physical nonlinearity in the constitutive equations has been investigated by Birman [START_REF] Birman | Effect of Physical Nonlinearity on Vibrations of Piezoelectric Rods[END_REF] and Birman and Suhir [START_REF] Birman | Effect of Material Nonlinearity on the Mechanical Response of Some Piezo-electric and Photonic Systems[END_REF]. It was shown that this nonlinearity has a significant effect on the amplitude of forced axial vibrations of the rod, particularly in the vicinity of resonant frequencies. For example, the ratio of amplitude obtained by the linear theory to its nonlinear counterpart ( ) R is shown as a function of the ratio of the frequency of the electric field to the fundamental frequency of axial vibrations ( ) F in Fig. 2. As follows from this figure, physical nonlinearity cannot be neglected in the vicinity of the resonance conditions. Moreover, the width of the spectrum of driving frequencies, where neglecting physically nonlinear effects is unacceptable, becomes broader for higher amplitude of the electric field.

The effect of the adjusted piezoelectric coefficient (as explained above, this coefficient reflects thermal contributions) is shown in Figs. 345for PZT-5H rods and three different amplitudes of the electric field. The results for PZN-4.5% rods are also shown in Fig. 6. In these figures, the amplitude of axial vibrations, i.e.W , is depicted as a function of the nondimensional driving frequency F . The results shown in Figs. 3456clearly illustrate that a higher temperature results in larger amplitude of forced axial vibrations of a rod (transducer). Predictably, the amplitude of motion increases with the amplitude of the electric field. However, a higher electric field does not change the qualitative effect of temperature.

CONCLUSIONS

The following conclusions can be drawn from the performed analysis.

• The paper presents a physically nonlinear formulation for a piezoelectric material, accounting for thermal effects. The formulation is applied to the solution of the problem of forced vibrations of a piezoelectric rod (transducer) driven by an alternating electric field applied in the axial direction. The solution is obtained by using the generalized Galerkin procedure, satisfying all the boundary conditions. As follows from numerical examples, physical nonlinearity cannot be disregarded when predicting the response of a rod in the vicinity of resonant frequency. An elevated temperature can result in a noticeable increase of the amplitude of motion.

Figures

Fig. 1. Piezoelectric rod with the coordinate systems (r-θ-z and 1-2-3) used in the analysis. The end z=0 is assumed clamped, while the end z=h is free.

θ (2) r (1) 2a h z (3)
.

Fig. 2 .

 2 Fig.2. Effect of voltage amplitude on the ratio of linear-to-nonlinear amplitudes of vibrations for a PZT-5H rod. Curves correspond to the amplitude voltage of 0.5MV/m ( ), 1.0MV/m ( ), and 1.5MV/m ( ), respectively.

Fig. 3 .

 3 Fig. 3. Effect of the piezoelectric coefficient 33 e on the amplitude of physically nonlinear vibrations of a PZT-5H rod. The amplitude of the electric field is equal to 0.5MV/m. Curves correspond to 33 e =24.5 C/m 2 ( ), 44.5 C/m 2 ( ), and 64.5C/m 2 ( ), respectively.

Fig. 4 .

 4 Fig. 4. Effect of the piezoelectric coefficient 33 e on the amplitude of physically nonlinear vibrations of a PZT-5H rod. The amplitude of the electric field is equal to 1.0MV/m. Curves correspond to 33 e =24.5 C/m 2 ( ), 44.5 C/m 2 ( ), and 64.5C/m 2 ( ), respectively.

Fig. 5 .

 5 Fig. 5. Effect of the piezoelectric coefficient 33 e on the amplitude of physically nonlinear vibrations of a PZT-5H rod. The amplitude of the electric field is equal to 2.0MV/m. Curves correspond to 33 e =24.5 C/m 2 ( ), 44.5 C/m 2 ( ), and 64.5C/m 2 ( ), respectively.

Fig. 6 .

 6 Fig. 6. Effect of the piezoelectric coefficient 33 e on the amplitude of physically nonlinear vibrations of a PZN-4.5% PT rod. The amplitude of the electric field is equal to 2.0MV/m. Curves correspond to 33 e =12.1 ( ), 22.1 ( ), and 32.1 C/m 2 ( ), respectively.

Table 1 .

 1 Material constant of representative piezoelectrics

	Material	33 e	( ) 2 m C	333 e	( ) 2 m C	C 33	( ) GPa	l 333	( ) m F
	PZT-5H	24.54 (or 31.02) 5.7*10 4	108.0	-21.21*10 -6
	PZN-4.5%PT	12.1 (or 13.02)	4.2*10 4	89.0		-1.61*10 -6

 TIMA EDITIONS / THERMINIC 2005 ISBN 2-916187-01-4

Acknowledgement: This research was partially supported by the Army Research Office Contract W911NF-04-1-0192. The support and encouragement of the Program Director, Dr. Gary L. Anderson are greatly appreciated.


Robert Darveaux

Effect of Assembly Stiffness and Solder Properties on Thermal Cycle Acceleration Factors

The impact of assembly stiffness and solder properties on thermal cycle acceleration factors was investigated using a simplified model of an electronic assembly. It was seen that stress -strain hysteresis loops change dramatically with modest changes in assembly stiffness or imposed strain. The impact on strain range or strain energy density saturates as the assembly stiffness significantly exceeds that of the solder.

The acceleration factors between various temperature cycle conditions vary with assembly stiffness and solder properties. Hence, standard qualification tests do not accelerate all types of surface mount components equally. Strain energy density appears to be a more robust damage indicator than strain range, based on its lower sensitivity to specific constitutive relations.

INTRODUCTION

Reliability prediction models have several important uses: 1) optimization of a product design before it is manufactured, 2) optimization of a product after it has already failed a qualification test, 3) estimation of field use reliability, and 4) establishing appropriate accelerated tests based on the field use conditions. In the case of solders used in electronic packaging, fatigue life under cyclic loading conditions is typically correlated to strain range per cycle [1- [START_REF] Perkins | Thermo-Mechanical Failure Comparison and Evaluation of CCGA and CBGA Elecronic Packages[END_REF][START_REF] Clech | An Extension of the Omega Method to Primary and Tertiary Creep of Lead-Free Solders[END_REF] or strain energy density per cycle [START_REF] Clech | Engineering Analysis of Thermal Cycling Accelerated Tests for Surface-Mount Attachment Reliability Evaluation[END_REF][START_REF] Perkins | Thermo-Mechanical Failure Comparison and Evaluation of CCGA and CBGA Elecronic Packages[END_REF][START_REF] Syed | Accumulated Creep Strain and Energy Density Based Thermal Fatigue Life Prediction Models for SnAgCu Solder Joints[END_REF][12][START_REF] Clech | Surface Mount Assembly Failure Statistics and Failure Free Time[END_REF][START_REF] Clech | A Comprehensive Surface Mount Reliability (CSMR) Model Covering Several Generations of Assembly Technology[END_REF][START_REF] Sarihan | Energy-Based Methodology for Damage and Life Prediction of Solder Joints under Thermal Cycling[END_REF][START_REF] Darveaux | Reliability of Plastic Ball Grid Array Assembly[END_REF][START_REF] Darveaux | Optimizing the Reliability of Thin Small Outline Package (TSOP) Solder Joints[END_REF][START_REF] Darveaux | Effect of Simulation Methodology on Solder Joint Crack Growth Correlation and Fatigue Life Prediction[END_REF][START_REF] Clech | Acceleration Factors and Thermal Cycling Test Efficiency for Lead-Free Sn-Ag-Cu Assemblies[END_REF][START_REF] Pang | Lead-Free 95.5Sn-3.8Ag-0.7Cu Solder Joint Reliability Analysis for Micro-BGA Assembly[END_REF][START_REF] Clech | An Extension of the Omega Method to Primary and Tertiary Creep of Lead-Free Solders[END_REF].

In all predictive models, accurate calculation of the stresses and strains in the solder joint is critical. In order to perform a structural calculation, constitutive models are needed which describe the relationship between stress, strain, strain rate, temperature, and microstructure. There are several forms of constitutive relations for solders in use today (see for example refs [START_REF] Sarihan | Energy-Based Methodology for Damage and Life Prediction of Solder Joints under Thermal Cycling[END_REF][START_REF] Pang | Lead-Free 95.5Sn-3.8Ag-0.7Cu Solder Joint Reliability Analysis for Micro-BGA Assembly[END_REF][START_REF] Lau | Acceleration Models, Constitutive Equations, and Reliability of Lead-Free Solders and Joints[END_REF][START_REF] Clech | An Obstacle-Controlled Creep Model for Sn-Pb and Sn-Based Lead-Free Solders[END_REF][START_REF] Clech | An Extension of the Omega Method to Primary and Tertiary Creep of Lead-Free Solders[END_REF][START_REF] Pang | Creep and Fatigue Characterization of Lead Free 95.5Sn-3.8Ag-0.7Cu Solder[END_REF][START_REF] Wiese | Creep of Eutectic SnAgCu in Thermally Treated Solder Joints[END_REF]).

Product qualification generally takes a different approach. The industry has established a series of qualification tests over the years, and each new generation of product is expected to pass the same number of cycles as the previous generation. In some cases, companies will relax certain criteria over time, but they are always reluctant due to liability concerns over field returns.

For example, temperature cycle qualification for cellular phone applications has traditionally been 1000 cycles -40C 125C, with 15min ramps, 15min dwells. Over time, this has been reduced in some cases to 300 to 600 cycles (although the "target" is still 1000 cycles). In order to make such a reduction, it is necessary to have a good understanding of the acceleration factors between the qualification test conditions and the field use conditions. In general, these acceleration factors will depend on both the assembly stiffness and the solder properties.

The present paper explores the interaction of assembly stiffness and solder properties with respect to thermal cycle acceleration factors. A 1-dimensional mechanical model is used for the analysis. Three solder alloys are compared: 62Sn36Pb2Ag, 97.5Pb2.5Sn, and 96.5Sn3.5Ag.

The properties of these alloys are representative of those being used in various flip chip, BGA, and SMT applications. In all simulations, both the inelastic strain range and the total strain energy density are used as damage indicators to predict fatigue life.

SIMPLIFIED MODEL OF AN ELECTRONIC ASSEMBLY

The mechanical model used in the present analysis is shown schematically in Figure 1. There are 4 elements which represent (1) the elastic response of the assembly, (2) the elastic response of the solder, (3) the time-independent plastic flow response of the solder, and (4) the time and temperature dependent creep response of the solder. When the assembly undergoes a temperature change, it is represented by an imposed strain on the system. Methods to estimate the assembly stiffness and imposed strain are discussed in refs [START_REF] Sauber | Predicting Thermal Fatigue Lifetimes for SMT Solder Joints[END_REF][START_REF] Darveaux | Constitutive Relations for Tin Based Solder Joints[END_REF] For example, in moving from the "initial" condition to the "ramp" condition in Figure 1, the assembly and solder springs stretch, and the solder undergoes both plastic flow and creep. The net strain of these 4 elements is equal to the imposed strain on the system. The relative magnitude of these effects depends on the stiffness of the assembly and solder elements, and on the plastic flow and creep properties for the solder.

During the "dwell" condition in Figure 1, the assembly and solder springs contract and the creep deformation of the solder increases. The overall strain is held fixed, and the stress relaxes. The magnitude of this stress relaxation depends on the stiffness of the assembly and solder elements, and on the creep properties of the solder.

After the dwell, the imposed strain is applied in the opposite direction. It is assumed that the solder recovers during the dwell, and the plastic flow and primary creep behavior start again from zero, i.e., there is no cyclic work hardening.

Evidence of this effect is given in Refs [START_REF] Darveaux | Effect of Simulation Methodology on Solder Joint Crack Growth Correlation[END_REF][START_REF] Darveaux | Shear Deformation of Indium Solder Joints[END_REF].

Figure 1. Simplified Model of an Electronic Assembly

CONSTITUTIVE RELATIONS FOR SOLDER

The constitutive relations for the solder used in the present study were taken from refs [START_REF] Sarihan | Energy-Based Methodology for Damage and Life Prediction of Solder Joints under Thermal Cycling[END_REF][START_REF] Wiese | Creep of Eutectic SnAgCu in Thermally Treated Solder Joints[END_REF]. The inelastic solder properties were measured on BGA joints loaded in shear. The values were then converted to equivalent tensile constants.

The elastic modulus is given by

E = E o -E 1 (T-273) (1)
where E o is the modulus at 273K, T is the absolute temperature, and E 1 gives the temperature dependence.

The steady state (or minimum) creep rate is given by

dε s dt = C ss [sinh(ασ)] n exp( -Q a kT ) (2)
where dε s /dt is the steady state strain rate, k is Boltzmann's constant, σ is the applied stress, Qa is the apparent activation energy, n is the stress exponent, α prescribes the stress level at which the power law dependence breaks down, and C ss is a constant.

Steady state creep is not generally achieved immediately when stress is applied. A certain amount of primary creep occurs before attaining steady state.

For normal decelerating primary creep, the strain rate starts high and decreases to a minimum value as the material work hardens. The instantaneous creep rate at constant stress and temperature is given by

dε c dt = t)) dt s d (-B exp B T + 1 ( dt s d (3)
where dε c /dt is the instantaneous creep rate, ε T is the transient creep strain, and B is the transient creep coefficient.

In addition to time dependent creep, some time independent plastic flow can occur in solders at high stresses. The following strain hardening law was used to describe time independent plastic flow

ε p = C p σ m ( 4
)
where ε p is the time-independent plastic strain, and C p and m are constants. The elastic strain is calculated by Hook's law

ε e = σ/Ε (5)
The total strain is given by the sum of elastic strain, creep strain, and plastic strain

ε t = ε e + ε c + ε p (6)
where ε t is the total strain.

The constitutive constants for the 3 solder alloys used in the present analysis are given in Table 1.

THERMAL CYCLE ACCELERATION FACTORS

An acceleration factor is defined as the number of cycles to failure under field use conditions relative to accelerated test conditions. AF = N field / N acc test [START_REF] Uegai | A Method of Fatigue Life Prediction for Surface-Mount Solder Joints of Electronic Devices By Mechanical Fatigue Test[END_REF] For a given joint size, solder fatigue is generally found to be inversely proportional to strain energy density to the power j. Hence when comparing 2 different thermal cycle conditions, the acceleration factor is given by

AF 1,2 = (∆W 2 / ∆W 1 ) j (8)
where AF 1,2 is the ratio of fatigue life for temperature cycle condition 1 versus condition 2, and ∆W is the strain energy density accumulated per cycle. The exponent, j, is typically found to be close to 1 [START_REF] Clech | Engineering Analysis of Thermal Cycling Accelerated Tests for Surface-Mount Attachment Reliability Evaluation[END_REF][START_REF] Perkins | Thermo-Mechanical Failure Comparison and Evaluation of CCGA and CBGA Elecronic Packages[END_REF][START_REF] Syed | Accumulated Creep Strain and Energy Density Based Thermal Fatigue Life Prediction Models for SnAgCu Solder Joints[END_REF][12][START_REF] Clech | Surface Mount Assembly Failure Statistics and Failure Free Time[END_REF][START_REF] Clech | A Comprehensive Surface Mount Reliability (CSMR) Model Covering Several Generations of Assembly Technology[END_REF][START_REF] Sarihan | Energy-Based Methodology for Damage and Life Prediction of Solder Joints under Thermal Cycling[END_REF][START_REF] Darveaux | Reliability of Plastic Ball Grid Array Assembly[END_REF][START_REF] Darveaux | Optimizing the Reliability of Thin Small Outline Package (TSOP) Solder Joints[END_REF][START_REF] Darveaux | Effect of Simulation Methodology on Solder Joint Crack Growth Correlation and Fatigue Life Prediction[END_REF][START_REF] Clech | Acceleration Factors and Thermal Cycling Test Efficiency for Lead-Free Sn-Ag-Cu Assemblies[END_REF][START_REF] Pang | Lead-Free 95.5Sn-3.8Ag-0.7Cu Solder Joint Reliability Analysis for Micro-BGA Assembly[END_REF][START_REF] Clech | An Extension of the Omega Method to Primary and Tertiary Creep of Lead-Free Solders[END_REF].

When using inelastic strain range as the failure indicator, it is typically found that the fatigue life is inversely proportional to the strain range to the k power. Hence the acceleration factor between two temperature cycle conditions is given by

AF 1,2 = (∆ε 2 / ∆ε 1 ) k ( 9
)
where ∆ε is the inelastic strain range. The exponent k is found to be in the range of 1 to 2 [1-10,23].

In the present study, five different temperature cycle conditions were considered, as shown in 2 is a typical stress -strain hysteresis loop generated by stepping through time and solving for the stress and strain at each time increment. The stress is considered to represent equivalent normal stress, but the positive and negative signs have no real meaning in the current context. In general, solder will be under shear + compressive stress on one side of the temperature cycle and shear + tensile stress on the other side. Whether tensile occurs on the hot or cold side of the chamber depends on the location in a joint and on the specifics of the assembly structure and expansion mismatches.

The inelastic strain range and strain energy density are defined in Figure 2. Notice that the inelastic strain range is not equal in both directions. This is called thermal ratcheting, and it results in a net change in joint shape with cycling (see for example refs [START_REF] Schubert | Fatigue Life Models for SnAgCu and SnPb Solder Joints Evaluated by Experiments and Simulation[END_REF][START_REF] Darveaux | Effect of Deformation Behavior on Solder Joint Reliability Prediction[END_REF]. The inelastic strain energy density is the area within the hysteresis loop.

In the present study, the total strain energy density (which includes the elastic component) will be used to calculate acceleration factors. It was found in ref [START_REF] Park | Reliability of Lead-Free Copper Columns in Comparison with Tin-Lead Solder Column Interconnects[END_REF] that total strain energy density gave a better correlation than inelastic strain energy density for cyclic bending tests that lasted between 10,000 and 100,000 cycles.

Subsequent work has shown that total strain energy density per cycle also correlates well to thermal cycle fatigue data that is in the range of 100 to 10,000 cycles.

Figure 2. Stress -inelastic strain hysteresis loop for 96.5Sn3.5Ag, Temperature cycle condition 3 (-40C 125C, 900 sec ramps, 900 sec dwells), 5.E5 psi assembly stiffness, 1.E-4/C imposed strain.

EFFECT OF ASSEMBLY STIFFNESS AND IMPOSED STRAIN

The impact of increasing assembly stiffness on stressstrain hysteresis loops is shown in Figure 3. With all other factors held constant, it is seen that the peak stresses, the strain range, and the strain energy density all increase dramatically with increases in assembly stiffness.

In BGA assemblies, things that increase stiffness are double sided assemblies and increasing motherboard thickness, die thickness, and mold cap thickness.

Reducing the pad area and the number of solder joints also increases the effective assembly stiffness compared to the joint array. In leaded assemblies, the lead design is a key parameter, in addition to the factors discussed above for BGAs. See refs [START_REF] Darveaux | Reliability of Plastic Ball Grid Array Assembly[END_REF][START_REF] Park | Reliability of Lead-Free Copper Columns in Comparison with Tin-Lead Solder Column Interconnects[END_REF][START_REF] Darveaux | Reliability of Area Array Solder Joints in Bending[END_REF][START_REF] Darveaux | Effect of Test Board Design on the 2 nd Level Reliability of a Fine Pitch BGA Package[END_REF] for examples on how measured fatigue life was reduced by such increases in assembly stiffness.

The effect of increasing imposed strain is shown in Figure 4. With all other factors held constant, it is seen that the peak stresses, the strain range, and the strain energy density all increase dramatically with increases in imposed strain. A 4X increase in imposed strain results in an 8X increase in inelastic strain range due to the highly nonlinear nature of the solder constitutive relations.

Factors that affect the imposed strain in a structure are die size, package size, joint height, and expansion mismatches between the materials. There are numerous examples in the literature showing these effects. A contour map showing how strain energy density per cycle is impacted by combinations of assembly stiffness and imposed strain is shown in Figure 5. For a given level of imposed strain, the strain energy density increases rapidly with stiffness, then starts to saturate at a maximum value. Saturation occurs as the stiffness of the assembly becomes significantly greater than the stiffness of the solder.

A contour plot for inelastic strain range is shown in Figure 6, which demonstrates similar trends to Figure 5. Figure 6. Inelastic strain range for 96.5Sn3.5Ag under temperature cycle condition 1 (-55C 125C, 180 sec ramps, 720 sec dwells).

INTERACTION OF TEMPERATURE CYCLE CONDITIONS AND SOLDER PROPERTIES

The impact of assembly stiffness on strain energy density and inelastic strain range is shown in Figures 789101112for five temperature cycle conditions and three alloys. The imposed strain was 1.E-4/C for all simulations. It is seen that there is a rapidly increasing portion to the curves, followed by a saturation portion as the assembly stiffness becomes greater than the solder stiffness.

The vertical separation between curves in these figures varies with stiffness. This indicates that the acceleration factor between two different thermal cycle conditions will also vary with assembly stiffness. In some regimes the acceleration factor is increasing with stiffness, and in other regimes it is decreasing.

This has important ramifications for product qualification. Typically, the exact same tests are used to do board level qualification of different components in a product. Since the effective assembly stiffness for the various components will be different, the equivalent number of field use years will also be different.

In comparing the strain energy density plots (Figures 7,9,11) to the strain range plots (Figures 8,10,12) it is seen that different acceleration factors will be predicted depending on the damage indicator. For example, the simulations for 62Sn36Pb2Ag show that for assembly stiffness less than 5.E5 psi, TC #5 is more damaging than TC #4 if inelastic strain range is considered (Figure 8). On the other hand, the opposite trend is predicted if strain energy density is used as the damage indicator (Figure 7).

Furthermore, a measurable difference in fatigue life is expected between temperature cycle profiles #1, #2, and #3 when considering strain energy density. However, curves for these conditions are virtually convergent when considering inelastic strain range.

Plots comparing the three alloys with respect to strain energy density per cycle are shown in Figures 13 and14.

It is seen that the trend with assembly stiffness is slightly different for each alloy under TC #1 conditions. However, the differences are much greater between alloys under TC #5 conditions. Hence, the acceleration factors between qualification and field use conditions will be alloy dependent. This effect is shown more clearly in Figures 15 and16, where the fatigue life ratio is estimated using strain energy density and inelastic strain range, respectively. An imposed strain of 1.E-4/C was used in the simulations. Again, the acceleration factor depends on alloy, stiffness, and damage indicator. The largest discrepancies between alloys are in the stiffness range of 5.E5 psi and below. Also, inelastic strain range will predict a larger variation between alloys than strain energy density. A comparison in fatigue life ratio between TC#1 and TC#2 is shown in Figures 17 and18. These conditions have the same temperature range and cyclic frequency, but vary in the ramp times and dwell times. A traditional Norris-Landsberg approach would predict the exact same life for both of these conditions (ratio = 1). The fatigue life ratio based on strain energy density is alloy and stiffness dependent and it varies from 0.95 to 1.25. If one considers inelastic strain range, the predicted life ratio is generally closer to 1, but still depends on alloy and stiffness.

EFFECT OF CONSTITUTIVE RELATIONS

The inelastic constitutive relations used in this study include time independent plastic flow, primary creep, and steady state creep. As discussed previously, there are several versions of constitutive relations being used in the industry to predict solder fatigue. A common practice is to only use steady state creep, since it is most conveniently available in finite element software. Also, there is rather large variation in the measured material constants depending on the investigators sample configuration and test method (see ref [START_REF] Norris | Reliability of Controlled Collapse Interconnections[END_REF] for more discussion).

Shown in Figure 19 are stress-strain hysteresis loops using progressively more complex constitutive relations. It is seen that only using steady state creep results in the largest stress range, but the smallest strain range. Adding primary creep, and then plastic flow progressively reduces the predicted stress range and increases the strain range.

The effect of constitutive relations on strain energy density and inelastic strain range is shown for the three alloys in Figures 202122232425. In all cases the imposed strain was 1.E-4/C.

In the stiffness range greater than 1.E6 psi, using only steady state creep tends to over-predict the strain energy density per cycle. At lower assembly stiffness, the strain energy density is under predicted. Using only steady state creep always under-predicts the inelastic strain range.

Figure 19.

Stress -Strain hysteresis loops for 96.5Sn3.5Ag solder with various constitutive models, temperature cycle condition 3 (-40C 125C, 900 sec ramps, 900 sec dwells), 4.0E5 psi assembly stiffness, 1.5E-4/C imposed strain.

Adding plastic flow becomes important in the low stiffness regime for 97.5Pb2.5Ag and 96.5Sn3.5Ag solders. In general, the calculated strain energy density is less sensitive to constitutive relations than inelastic strain range is.

DISCUSSION

The present analysis has shown that thermal cycle acceleration factors depend on several factors: 1) imposed strain, 2) assembly stiffness and 3) solder properties. This has ramifications with respect to product qualification because the industry uses a common set of tests to qualify product. Unfortunately, it is probably not practical to tailor qualification tests to assembly stiffness or solder alloy type.

The magnitude of calculated acceleration factors also depends on the damage indicator that is chosen, i.e., strain range vs. strain energy density. The analysis community is somewhat polarized on this issue. People have their preference for using either strain range or strain energy density. The preference is usually based on historical practice within a company, believed theoretical "correctness", or academic lineage.

Based on the current analysis, it is clear that changing temperature cycling conditions is a good way to delineate which method gives the most accurate prediction. For example, it is expected that increasing ramp rate will reduce fatigue life based on strain energy density considerations, but have little effect on fatigue life based on strain range considerations. In reviewing the measured data in refs [START_REF] Darveaux | Shear Deformation of Lead Free Solder Joints[END_REF][START_REF] Fan | Effects of Dwell Time and Ramp Rate on Lead-Free Solder Joints in FCBGA Packages[END_REF], fatigue life did decrease when ramp rate was increased significantly. This would tend to support strain energy density as a better damage indicator.

Solder alloy properties also have an impact on thermal cycle acceleration factors. Life predictions using strain range are more sensitive to variations in constitutive relations than those using strain energy density. Hence, it would follow that strain energy density is a more robust damage indicator.

CONCLUSIONS

1) The impact of assembly stiffness and solder properties on thermal cycle acceleration factors was investigated using a simplified model of an electronic assembly.

2) Stress -strain hysteresis loops change dramatically with modest changes in the assembly stiffness or the imposed strain.

3) The impact on strain range and strain energy density saturates as the assembly stiffness significantly exceeds that of the solder.

4) The acceleration factors between various temperature cycle conditions vary with assembly stiffness and solder properties. Hence, standard qualification tests do not accelerate all types of surface mount components equally. 5) Strain energy density appears to be a more robust damage indicator than strain range, based on its lower sensitivity to specific constitutive relations.

Figure 3 .

 3 Figure 3. Stress -inelastic strain hysteresis loops for 96.5Sn3.5Ag, Temperature cycle condition 3 (-40C 125C, 900 sec ramps, 900 sec dwells), 1.E-4/C imposed strain.

Figure 4 .

 4 Figure 4. Stress -inelastic strain hysteresis loops for 96.5Sn3.5Ag, Temperature cycle condition 3 (-40C 125C, 900 sec ramps, 900 sec dwells), 5.E5 psi assembly stiffness.

4 Figure 5 .

 45 Figure 5. Strain energy density per cycle for 96.5Sn3.5Ag under temperature cycle condition 1 (-55C 125C, 180 sec ramps, 720 sec dwells).

Figure 7 .

 7 Figure 7.Strain energy density per cycle for 62Sn36Pb2Ag under various temperature cycle conditions, 1.E-4/C imposed strain.

Figure 8 .

 8 Figure 8. Inelastic strain range for 62Sn36Pb2Ag under various temperature cycle conditions, 1.E-4/C imposed strain.

Figure 9 .

 9 Figure 9. Strain energy density per cycle for 97.5Pb2.5Sn under various temperature cycle conditions, 1.E-4/C imposed strain.

Figure 10 .Figure 11 .Figure 12 .Figure 13 .

 10111213 Figure 10. Inelastic strain range for 97.5Pb2.5Sn under various temperature cycle conditions, 1.E-4/C imposed strain.

Figure 14 .

 14 Figure 14. Comparison of strain energy density per cycle for three alloys, temperature cycle condition 5 (25C 75C, 600sec ramps, 13800 sec dwells), 1.E-4/C imposed strain.

Figure 15 .Figure 16 .Figure 17 .Figure 18 .

 15161718 Figure 15. Fatigue life ratio based on (Strain Energy Density) -1 for temperature cycle condition 1 (-55C 125C, 180 sec ramps, 720 sec dwells) vs. condition 5 (25C 75C, 600 sec ramps, 13800 sec dwells), 1.E-4/C imposed strain.

Figure 25 .

 25 Figure 25. Effect of constitutive relations on Inelastic Strain Range for 96.5Sn3.5Ag solder under temperature cycle condition 4 (0C 100C, 600 sec ramps, 300 sec dwells) and 1.0E-4/C imposed strain.

Table 1 .

 1 Constitutive constants for 3 solder alloys.

	Parameter 96.5Sn3.5Ag 97.5Pb2.5Sn 62Sn36Pb2Ag
	E o (Mpsi)	7.6	3.5	5.1
	E 1 (Kpsi/K)	27	4.1	22
	C ss (1/sec)	1.42e5	2.66e11	8.03E4
	α (1/psi)	3.64e-4	4.63e-4	4.62e-4
	n	5.5	7.0	3.3
	Qa (eV)	0.75	1.15	0.70
	ε T	0.096	0.066	0.023
	Β	227	237	263
	C p (psi -m )	1.20e-18	8.16e-15	1.19E-23
	m	4.39	3.10	5.53

Table 2 .

 2 Profile #1 is a typical package level qualification test. Profiles #3 and #4 are typically used in board level qualification studies. Profile #5 is meant to represent a field use condition for a computing application.

			Table 2
		Simulated Thermal Cycle Profiles
	Profile	Low	High	Ramp	Dwell
		Temp	Temp	Time	Time
		(C)	(C)	(sec)	(sec)
	1	-55	125	180	720
	2	-55	125	720	180
	3	-40	125	900	900
	4	0	100	600	300
	5	25	75	600	13800
	Shown in Figure

psi) Inelastic Strain Range Steady State + Prim ary Creep + Plastic Flow Steady State + Prim ary Creep Steady State Creep

			Steady State + Prim ary Creep + Plastic Flow
			Steady State + Prim ary Creep
			Steady State Creep
		100
	Strain Energy Density per	Cycle (psi)	10
			1
			1.E+05	1.E+06	1.E+07
				Assembly Stiffness (psi)
	Figure 20. Effect of constitutive relations on Strain
	Energy Density per cycle for 62Sn36Pb2Ag solder under
	temperature cycle condition 4 (0C 100C, 600 sec
	ramps, 300 sec dwells) and 1.0E-4/C imposed strain.
	4000
	2000	1.E-02
	0
	-2000
	-4000	1.E-03
	-6000
	-8000
	-10000
	-0.	1.E-04
			1.E+05	1.E+06	1.E+07
	Assembly Stiffness (Figure 21. Effect of constitutive relations on Inelastic
	Strain Range for 62Sn36Pb2Ag solder under temperature
	cycle condition 4 (0C 100C, 600 sec ramps, 300 sec
	dwells) and 1.0E-4/C imposed strain.

006 -0.004 -0.002 0.000 0.002 0.004 0.006 Inelastic Strain Stress (psi) Steady State Creep Primary Creep + Steady State Creep Plastic Flow + Primary Creep + Steady State Creep

			Steady State + Prim ary Creep + Plastic Flow
			Steady State + Prim ary Creep
			Steady State Creep
		100
	Strain Energy Density per	Cycle (psi)	10
			1
			1.E+05	1.E+06	1.E+07
				Assembly Stiffness (psi)
	Figure 22. Effect of constitutive relations on Strain
	Energy Density per cycle for 97.5Pb2.5Sn solder under
	temperature cycle condition 4 (0C 100C, 600 sec
	ramps, 300 sec dwells) and 1.0E-4/C imposed strain.
		1.E-02
		1.E-03
		1.E-04
			1.E+05	1.E+06	1.E+07
				Assembly Stiffness (

psi) Inelastic Strain Range Steady State + Prim ary Creep + Plastic Flow Steady State + Prim ary Creep Steady State Creep

			Steady State + Prim ary Creep + Plastic Flow
			Steady State + Prim ary Creep
			Steady State Creep
		100
	Strain Energy Density per	Cycle (psi)	10
			1
			1.E+05	1.E+06	1.E+07
				Assembly Stiffness (psi)
	Figure 24. Effect of constitutive relations on Strain
	Energy Density per cycle for 96.5Sn3.5Ag solder under
	temperature cycle condition 4 (0C 100C, 600 sec
	ramps, 300 sec dwells) and 1.0E-4/C imposed strain.
		1.E-02
		1.E-03
		1.E-04
			1.E+05	1.E+06	1.E+07
				Assembly Stiffness (
	Figure 23. Effect of constitutive relations on Inelastic
	Strain Range for 97.5Pb2.5Sn solder under temperature
	cycle condition 4 (0C 100C, 600 sec ramps, 300 sec
	dwells) and 1.0E-4/C imposed strain.

psi) Inelastic Strain Range Steady State + Prim ary Creep + Plastic Flow Steady State + Prim ary Creep Steady State Creep

 TIMA EDITIONS / THERMINIC 2005 ISBN 2-916187-01-4


Raúl Antón

Miguel Castiella

Hans Jonsson

Bahram Moshfegh

SMOKE AND CFD VISUALIZATION OF THE FLOW AFTER AN EMC SCREEN IN A SUBRACK MODEL

Keywords: RNG k-ε, perforated plate, porosity, EMC screen, subrack, flow pattern, inclined flow, smoke-wire, flow visualization

An experimental parametric study was performed in order to study the flow pattern inside a radio base station subrack model. The parameters that were analysed are the velocity, the EMC screen porosity and the subrack geometry. The technique used to visualize the flow pattern is the smoke-wire technique.

The objective is not only to gain insight about the influence of those parameters in the flow pattern of the subrack but also to validate the flow pattern predictions of a detailed CFD model. A further objective is to show the importance of using an accurate CFD model in order to visualize the flow pattern in a situation in which there is a 90 degree turn followed by a perforated plate. The fact that the airflow impinges on the screen in an inclined way means that when the flow collides with the interior face of the hole, it suffers a sudden change in direction, even if the perforated plate is very thin. This change in direction is very important in order to get the right flow pattern after the screen.

INTRODUCTION

Enclosing electronics in sealed metal boxes provides good electromagnetic shielding, but evidently restricts the air movement necessary for adequate cooling. In order to meet electromagnetic compatibility requirements, a perforated plate (an EMC screen) must be used. The design of this screen must provide a sufficient free area ratio for the adequate airflow, but at the same time the holes must be small enough to block electromagnetic radiation. According to Ott [START_REF] Ott | Noise Reduction Techniques in Electronics Systems[END_REF], the maximum aperture size should be less than 1/20 th of the source wavelength.

Nowadays, thermal management engineers use Computational Fluid Dynamics (CFD) as a design tool to predict the thermal and hydraulic pattern of the airflow at the component level (see Rodgers et al. [START_REF] Rodgers | An Experimental Assessment of Numerical Predictive Accuracy for Electronic Component Heat Transfer in Forced Convection Part 1: Experimental Methods and Numerical Modeling[END_REF]), as well as at a system level (Lee et al. [START_REF] Lee | Application of CFD Technology to Electronic Thermal Management[END_REF] and Lee & Manhalingam [START_REF] Lee | Application of a CFD Tool for System-Level Thermal Simulation[END_REF]). There are still gaps in system level modelling, such as the uncertainty in input parameters like the screen loss coefficient (Joshi et al. [START_REF] Joshi | Challenges in Thermal Modeling of Electronics at the System Level: Summary of Panel Held at the Therminic 2000[END_REF]). Many turbulence models have been developed and introduced to the market. It is still required to validate those turbulence models experimentally. Smoke visualization has always had an extended use for visualizations at not high velocity levels. Lohan et al. [START_REF] Lohan | Visualization of forced air flows over a populated printed circuit board and their impact on convective heat transfer[END_REF] studied the flow characteristics around a electronic component; among other techniques, they used the smoke-wire method for flow visualization.

Several modelling techniques have been developed and are used at different stages of the design cycle of the product, as is explained by Minichielo and Belady [START_REF] Minichiello | Thermal Design Methodology for Electronic Systems[END_REF]. One possibility is to model several parts in the system level by means of a hydraulic volume or surface resistance as in Nevelsteen [START_REF] Nevelsteens | Screen Characterization under Fan Induced Swirl Conditions[END_REF], where the screens close to axial fans were modelled by a compact model with three directional loss coefficients. Baelmans et al. [START_REF] Baelmans | Flow Modeling in Air-Cooled Electronic Enclosures[END_REF] analyse experimentally the distance of influence by the screen, remarking on the difficulty of predicting the flow beside components placed close to the screen when the screen is modelled by means of hydraulic surface resistance.

In the models that use surface impedance to model the pressure drop, it may not be clear which pressure loss coefficient to use and which experimental coefficients may be needed. In a typical Radio Base Station (RBS) case, the flow impinges on the screen at different velocities and angles, and it is difficult to choose the proper loss coefficients from the literature. Further, it is not clear which porosity characterizes the model. In most cases the airflow goes through the inlet opening, then makes a 90º turn (at different velocities) and finally goes through the EMC screen. Quite a lot of complex flow phenomena (such as recirculation) may arise due to the geometry of the subrack. Certainly, the result of modelling the screen could in turn be used to develop correlations for pressure loss coefficients for different geometries.

The objective of this paper is to perform a parametric experiment in order to study the flow pattern inside the subracks of a radio base station model. The parameters that are analysed are the velocity, the screen porosity and the subrack geometry. The technique was used to visualize the flow pattern is the smoke-wire technique. Further, the objective is not only to gain insight into the influence of those parameters in the flow pattern of the subrack but also to compare with the flow pattern of a detailed CFD model of a subrack slot. The detailed CFD model used in this paper has been already validated by the authors in [START_REF] Antón | Modelling of EMC Screens for Radio Base Stations Part 1: Experimental Parametric Study[END_REF], [START_REF] Antón | Modelling of EMC Screens for Radio Base Stations Part 2: Evaluation of Turbulence Models[END_REF] and [START_REF] Antón | Detailed CFD Modelling of EMC Screens for Radio Base Stations: a Benchmark Study[END_REF]. In [START_REF] Antón | Detailed CFD Modelling of EMC Screens for Radio Base Stations: a Benchmark Study[END_REF], the authors found a good agreement between predicted and observes values of the static pressure drop and local velocities at several locations before and after the screen. But the overall flow pattern was not studied. It is important to use a correct CFD model in order to visualize the flow pattern in a situation in which there is a 90 degree turn and a perforated plate. The fact that the airflow impinges on the screen in an inclined way means that when the flow collides with the interior face of the hole, it suffers a sudden change in direction, even if the perforated plate is very thin. This change in direction is very important in order to get the right flow pattern after the screen. This sudden change in direction is badly predicted if hydraulic surface impedance is used (Anton et. al. [START_REF] Antón | Modelling of EMC Screens for Radio Base Stations Part 2: Evaluation of Turbulence Models[END_REF]).

EXPERIMENTAL SET-UP

The experimental validation was performed in a subrack model inside a wind tunnel. The general layout of the wind tunnel is shown in Figure 1.

The wind tunnel walls are made of Plexiglass. The air velocity is controlled by frequency regulation of the electric power driving the fan. A subrack was modelled by means of Plexiglass plates (4 mm thick), as is shown in Figure 2 and Figure 3 (all the dimensions in mm). There are a total of 14 PCB dummies. No heat loads were mounted on the PCB dummies for simplification of the experimental setup and it can be assumed that this will not greatly affect the flow pattern. An entrance duct was placed before the inlet of the subrack in order to get a fully developed turbulent flow at the inlet of the subrack and thus low turbulence intensity levels. Further a perforated plate at the inlet of the entrance duct was used as a turbulence generator in order to decrease the required length needed to create a fully developed turbulent flow at the inlet of the subrack. Stainless steel screens (1 mm thickness) with 60° staggered circular holes were used. The location of the screen is indicated in Figure 2. Three different porosities, ε, were used for the screen. The porosity in the table below is an approximate value that is calculated using

wd D n 4 / 2 π ε = .
Where n is the number of holes in a PCB slot, D is hole diameter, d is the subrack depth (mm) and w is the distance between two PCBs. Two models were built with the same dimensions except for the depth, d (see Figure 2, top view). In Table 3 the magnitude of the average inlet velocity is shown for the two models. The average is an approximate value obtained by doing a weighted average from the hot wire measurements along the inlet width of the subrack [START_REF] Antón | Modelling of EMC Screens for Radio Base Stations Part 1: Experimental Parametric Study[END_REF]. Five wires were placed in the middle between two PCB dummies that are located in the middle of the subrack model. In Figure 4 and Table 4 the distance of the wires from the screen is shown: The smoke-wire technique consists of coating the stretched wires with oil and heating them up. After heating them, the oil evaporates and immediately after that the oil condenses, forming smoke-like microdroplets. The smoke was illuminated by xenon light. The visualized flow pattern was recorded on video.

The objective of the videos was to identify the frontier line between the bulk flow with a high velocity level and the reversed flow (recirculation) with a low velocity level (see Figure 4). This frontier line is a function of several parameters. Some of those parameters are analysed: the velocity level, the ratio between the inlet and the depth of the subrack, and the screen thickness. The frontier line is identified by 5 points (see w1-w5 in Figure 4). These points are defined for having the zcomponent of the velocity equal to zero. In order to making it easier to see the frontier point, a black paper was placed on the PCB dummy as in Figure 5. First it was studied whether the position of the camera influenced the measured frontier points. After 60 videos it was concluded that this factor has small importance due to the fact that the wires are close to the PCB. Actually, the uncertainty of the camera position was estimated to be up to about 0.5 cm (twice the standard deviation of the values observed in the videos).

Furthermore, in order to decrease the uncertainty each wire was recorded between 2 and 5 times.

GEOMETRY AND BOUNDARY CONDITIONS OF THE DETAILED CFD MODEL

Due to symmetry conditions, only half of the PCB slot was modelled with the same dimensions as in the experimental model. The dimensions of half of the slot are shown in Figure 6. The flow is considered steady, isothermal and turbulent.

At the inlet of the entrance duct a constant velocity was set and constant turbulence intensity equal to 10%. The result of using constant velocity and turbulence intensity profiles (always 10%) at the entrance duct inlet (see Figure 2) was validated by a comparison between CFD and experimental values of the velocity and turbulence intensity profiles at the inlet of the subrack in a previous paper by the authors [START_REF] Antón | Detailed CFD Modelling of EMC Screens for Radio Base Stations: a Benchmark Study[END_REF]. The air density is always equal to 1.2 kg/m 3 . For the high velocity level (V5), a small pressure gradient was measured at the outlet. However, at the low velocity level (V1) there was no pressure gradient and a constant pressure outlet boundary condition was used. In Table 5, the pressure difference (in Pa) measured and used as a boundary condition in the simulations with the high velocity level is shown. It is necessary to supply boundary conditions for the possible reversed flow at the outlet. Those boundary conditions were always set using a turbulence intensity level equal to 10% and the hydraulic diameter of the outlet.

Symmetry

MESH DENSITY DISTRIBUTION AND SENSITIVITY

A three-dimensional perforated plate was modelled with circular holes (the holes were done with polygons of 12 sides with an area equal to that of the real hole); the arrangement was 60° staggered and the screen thickness 1 mm (see Figures 2 and3). The Gambit 2.1 [START_REF]Fluent 6.1[END_REF] grid generation package was used to generate the mesh. Partial views of the detailed model are shown in Figure 7. The mesh density distribution and sensitivity as well as the overall convergence of the CFD runs were analyze by the authors in [START_REF] Antón | Detailed CFD Modelling of EMC Screens for Radio Base Stations: a Benchmark Study[END_REF].

Figure 7. Mesh of the model, partial views.

NUMERICAL ACCURACY

FLUENT 6.1 [START_REF]Gambit 1.2. Gambit Manuals[END_REF] was used for the simulations. The governing equations were solved with a segregated scheme. The momentum, turbulent kinetic energy and turbulent energy dissipation equations were solved with second order upwind schemes. The pressure-velocitycoupling algorithm SIMPLE was used to solve the continuity equation. The following under-relaxation factors were used for the pressure, momentum, turbulent kinetic energy and the dissipation rate of turbulence kinetic energy: 0.3, 0.7, 0.8 and 0.8. The simulation was conducted on a workstation (2.5 GHz). The turbulence model employed was the RNG kε model by Yakhot et al. [START_REF] Yakhot | Development of Turbulence Models for Shear Flows by a Double Expansion Technique[END_REF]. In [START_REF] Antón | Detailed CFD Modelling of EMC Screens for Radio Base Stations: a Benchmark Study[END_REF] an evaluation of the turbulence model was performed, which concluded that the RNG kε model was the one that best predicted the pressure drop and the velocity magnitude at several locations.

PARAMETRIC STUDY EFFECT OF THE SCREEN POROSITY IN THE FLOW PATTERN

This section is based on the smoke visualizations. In Figure 8, it can be seen how the frontier line moves to the left when the screen porosity decreases, as expected. An interesting aspect is that the smoke visualization gives us qualitative and quantitative information (from the frontier points, see Figure 4). Three cases are shown: a) minimum velocity and short case, b) maximum velocity and short case and c) maximum velocity and large case. In the three cases, it is shown that the bulk flow increases by 2 cm when the screen porosity decreases from porosity 1 to porosity 3 (see also Table 1). Also it is apparent that the ratio between the inlet and the outlet is an important parameter, as it is seen when comparing cases (b) and (c) in Figure 8.

EFFECT OF THE VELOCITY IN THE FLOW PATTERN

This section presents the results from the smoke visualizations. In Figure 9, it can be seen how the frontier line moves to the left when the velocity decreases. There is a large difference between the lowest velocity and the rest of the velocities. Probably the flow for the minimum velocity is in transition, not fully turbulent and that is the reason why the bulk flow (see also Figure 4) is wider at the outlet. The reason why the bulk flow does not distribute across the outlet area is believed to be due to the large inertia of the flow. In the case of the smallest velocity, the inertia forces are not large enough to maintain this uneven distribution and the flow is hence distributed more evenly. The last three velocity levels produce similar frontier lines and thus the velocity has negligible influence on the frontier line.

VALIDATION OF THE FLOW PATTERN PREDICTED BY THE DETAILED CFD MODEL

In the following figures some smoke pictures can been seen of the frontier points for each of the wires and the streak lines that are generated by the smoke. Also, the path lines predicted by the detailed CFD model are shown and the quantitative comparison between the frontier points observed in the smoke visualizations (the average of 5 videos) and the values predicted by CFD. The uncertainty in the observed values is around 0.5 cm (twice the standard deviation of the observed values of the five videos that were recorded per wire), for the wire 5; at the outlet and for the large case, the estimated uncertainty increases up to 0.5-1 cm. In the smoke streak lines and in the CFD path lines, only the first 13.5 cm along the depth of the PCB is shown (see Figure 10). In Figure 11 a similar validation for the same case as before is shown, but now at the maximum velocity level (V5). In Figure 13 shows one case for the large model and a low velocity level (V2). In this case the disagreement between the predicted and observed values is a bit larger. One reason for this is that for the last three wires the flow is quite inclined (not at all perpendicular to the wire around the frontier point) which increases the uncertainty in the smoke visualization. It can be concluded from Figures 10 through 13 that the detailed CFD model is able to predict the right change in direction when the inclined flow collides with the interior face of the holes and thus it is a model that could be used to predict how large the bulk flow area is (see Figure 4), which is an important parameter when placing electronic components on a PCB.

CONCLUSIONS AND FUTURE WORK

The smokewire technique was applied in a subrack model. The frontier line between the bulk flow (high velocity level) and the reversed flow (low velocity level) has been identified. A parametric study that gives an insight into the flow pattern after the screen was performed with the following parameters: velocity, screen porosity and ratio between inlet and depth of the subrack.

Further, a comparison of the overall flow pattern after the screen for a detailed CFD model was showed. It is apparent that the detailed CFD model has the capacity to predict accurately the flow pattern after the screen. This model will be used in a larger parametric study in which the bulk flow area (or wetted area) will be investigated.

Figure 1 .

 1 Figure 1. Wind tunnel layout.

Figure 2 .

 2 Figure 2. Subrack model.

Figure 3 .

 3 Figure 3. Model of the screen and PCB slot.

Figure 4 .

 4 Figure 4. Location of the wires after the screen

Figure 5 .

 5 Figure 5. Detail of the PCB dummy. A black paper background was used in order to get a good contrast with the smoke and to distinguish the frontier points.

Figure 6 .

 6 Figure 6. Boundary conditions and geometry of the model (half of the slot between two PCBs is presented).

Figure 8 .

 8 = 200 mm, Velocity = 1.60 m/s b) Depth = 200 mm, Velocity = 7.50 m/s c) Depth = 260 mm, Velocity = 7.7 m/s Effect of the variation of the porosity on the flow pattern based on the smoke visualizations (x and z coordinate in Figure 4).

  75 100 125 150 175 200 PCB depth (mm), x coordinate PCB height (mm)75 100 125 150 175 200 PCB depth (mm), x coordinate PCB height (mm)mm, Porosity 1= 62.2% Depth = 200 mm, Porosity 3= 46.9 %Depth = 260mm, Porosity 1= 62.2% Figure9. Effect of the variation of the velocity on the flow pattern, based on the smoke visualizations (x and z coordinate in Figure4).

Figure 10 .

 10 Figure 10. Validation for the following case: depth=200 mm, V1= 1.60m/s, porosity 1= 62.2%. See x and z coordinate in Figure 4.

  75 100 125 150 175 200 PCB depth (mm), x coordinate PCB height (mm), z coordinate CFD SmokeSmoke visualizationPath lines predicted by CFD CFD vs. Smoke Figure11. Validation for the following case: depth = 200 mm, V5 =7.50 m/s, porosity 1= 62.2% (x and z coordinate in Figure4).

  the same case as before is shown, but now for the minimum porosity (porosity 3) depth (mm), x coordinate PCB height (mm), z coordinate CFD SmokeSmoke visualizationPath lines predicted by the CFD model CFD vs. Smoke Figure12. Validation for the following case: depth = 200 mm, V5 =7.50 m/s, porosity 3= 46.9 %. See x and z coordinate in Figure4.

Figure 13 .

 13 Figure 13. Validation for the following case: depth = 260 mm, V2 = 3.2 m/s, porosity 1= 62.2%. See x and z coordinate in Figure 4.

Table 1 .

 1 Screen characteristics (see Figure3)

		D (mm)	S (mm) D+S (mm)	ε (%)
	Porosity 1	3.75	0.75	4.5	62.2
	Porosity 2	3.50	1.00	4.5	54.3
	Porosity 3	3.25	1.25	4.5	46.9

Table 2

 2

. Subrack depth for the two models Subrack depth, d (mm) Short case 200 Large case 260

Table 3 .

 3 Average inlet velocities for the two sizes

		Short model (m/s)	Large model (m/s)
	V1	1.60	1.65
	V2	3.05	3.22
	V3	4.43	4.68
	V4	5.86	6.14
	V5	7.50	7.72

Table 4 .

 4 Distance of the wires after the EMC screen (in mm)

	Wire 1 Wire 2 Wire 3 Wire4 Wire 5
	16.5	81	146.5	210	276

Table 5

 5

	Pressure difference at the outlet and maximum velocity
		(V5)
	Subrack depth (mm)	200	200	200	260
	Porosity (%)	46.9	54.3	62.2	62.2
	Pressure	1.06	0.915	1.5	1.78

 TIMA EDITIONS / THERMINIC 2005 ISBN 2-916187-01-4

ACKNOWLEDGMENTS

The authors gratefully acknowledge the funding received from KK-Foundation, Ericsson, Nokia, Emerson Energy Systems and the Antonio Aranzabal Foundation.


Y Ezzahri

S Grauby

S Dilhaire

J M Rampnoux

W Claeys

) Y Zhang

A Shakouri

Determination of Thermophysical Properties of Si/SiGe Superlattices with a Pump-Probe Technique

A Time Domain Thermo-Reflectance technique (TDTR) is used for the non-destructive evaluation of thermal and acoustic properties of Si/SiGe superlattices (SL) at room temperature. In particular, this technique allows the determination of the SL cross-plan thermal conductivity, the thermal boundary resistance at the interface between the metallic film and the SL, and the longitudinal sound velocity inside the SL. We discuss the effect of both transducer thickness and pump power on the measurements and then on values of the extracted thermophysical properties of the superlattice.

I. Introduction

The control and the stabilization of temperature are essential in microelectronics and optoelectronics, because quality and reliability of IC chips are heavily dependent on chip temperature. The temperature control of these components is typically accomplished with thermoelectric (TE) coolers. TE coolers have become essential in modern optical telecommunications to control the characteristics of laser sources, switching/routing elements, and detectors used in wavelength division multiplexed systems. Cooling requirements in microelectronic applications have also risen dramatically in recent years due to the increase in clock speed and the reduction in feature size. Generally, as these devices have become smaller, faster, and denser, the power density has greatly increased. Conventional Bi 2 Te 3 TE coolers are incompatible with integrated circuits fabrication technologies, essentially silicon technology, and therefore limited in how small they can be manufactured [START_REF] Labounty | Design and Characterization of thin film microcoolers[END_REF]. Thus it will be very interesting if we could use silicon as efficient thermoelectric material to realize on-chip cooling solution. The existing barrier is that silicon is a very poor thermoelectric material with very low figure of merit (Z). Fortunately, recent studies [START_REF] Hicks | Effect of quantum-well structures on the thermoelectric figure of merit[END_REF][START_REF] Sofo | Thermoelectric figure of merit of superlattices[END_REF][START_REF] Broido | Effect of superlattice structure on the thermoelectric figure of merit[END_REF] showed the possibility to enhance this parameter by introducing low dimensional structures. This has been supported by development of epitaxial growth techniques such as molecular beam epitaxy and special kinds of chemical vapor phase epitaxy that allow the control of individual semiconductor layers on an atomic scale. This achievement provides a strong basis for designing semiconductor materials with tailored band structures and consequently tailored electronic and optical properties. The combination of different composites in the form of ultra thin layers have opened up exciting possibilities for novel devices, and it has also increased basic understanding of electronic behavior in reduced dimensions. Thermoelectricity was one of the first domains to be interested in these structures. Several measurements [START_REF] Lee | Thermal conductivity of Si-Ge superlattices[END_REF][START_REF] Yang | Measurements of anisotropic thermoelectric properties in superlattices[END_REF][START_REF] Huxtable | Thermal conductivity of Si/SiGe and SiGe/SiGe superlattices[END_REF][START_REF] Touzelbaev | Thermal characterization of Bi 2 Te 3 /Sb 2 Te 3 superlattices[END_REF] have been done on different SLs in order to study their thermal and electric properties, and all agreed on the fact that these structures have a high thermoelectric power and a low thermal conductivity, compared to bulk or solid solution alloys.

Si doped Ge must receive a particular attention considering its physical properties and its very easy monolithic integration with microelectronics and optoelectronics components, this kind of SL could hence be a potential candidate to replace the conventional Bi 2 Te 3 thermoelectric materials used until now. In fact SiGe has become one of the most attractive thermoelectric materials especially for high temperature applications. In this paper, we use TDTR technique to study the thermal conductivity of two Si/SiGe superlattices at room temperature. This technique has a temporal resolution <1ps, which is necessary to the study of fast physical phenomena, in particular thermal phenomena. Moreover, it has a nanometer-scale depth resolution [START_REF] Cahill | Nanoscale thermal transport[END_REF], and therefore is well suited for measurements on very thin layers, such as SLs. More specifically, this method can isolate the effects of interface conductance from the thermal conductivity of a thin layer [START_REF] Capinsky | Thermalconductivity measurements of GaAs/AlAs superlattices using a picosecond optical pump-andprobe technique[END_REF].

The main objective of our work is to determine simultaneously at room temperature, the effective cross-plan thermal conductivity of Si/SiGe SL, Kapitza resistance at the interface transducer/SL, and the longitudinal sound velocity inside the SL. Besides, we discuss the effect of the transducer thickness and the effect of the pump power on the temperature decay into the structure after absorption of the pump energy. These effects directly influence the quality of the fit, and then the value of the extracted thermophysical properties.

II. Sample description and Experimental set-up

Two samples of 1µm Si/Si 0.7 Ge 0.3 and Si/Si 0.75 Ge 0.25 SL are grown on about 500µm silicon substrate, and coated by 100nm and 12nm Al film respectively which serves as a transducer. The term "transducer" is used to specify the role of the metallic film: it converts light energy into heat which creates an acoustic wave propagating deeply in the active structure. By this way, the transducer creates the acoustic vector necessary for the analysis of the acoustic contribution in the pumpprobe signal. Between the SL and the substrate is grown a 2µm buffer layer to reduce mechanical stress between both materials due to lattice mismatch between silicon and germanium which is of about 4.2% [START_REF] Douglas | Si/SiGe heterostructures: from material and physics to devices and circuits[END_REF]. Experimental set-up is shown in figure 1. A laser beam from a 80MHz Ti: Sapphire laser with a Full Width at Half Maximum (FWHM) of 100 fs is separated into two beams with an intensity ratio of 10:1 by a nonpolarizing beam splitter. The intense pump beam is used to heat the film, while the weaker probe beam is used to monitor the reflectivity variations of the film surface which can be related to its temperature variations. The pump beam passes through an Acousto-Optic Modulator (AOM) that creates a pulse train at 574 kHz. The pump beam is focused to about 20µm. The probe beam is reflected off a retroreflector mounted on a mechanical translation stage to increase the optical path length of the probe beam and hence the time delay between the pump and the probe pulses. The probe, which is centered in the heated area created by the pump pulse, is focused to about 6µm at near normal incidence. A polarizer, oriented such that only the probe beam can cross it, is positioned before the photodiode which monitors the probe beam reflection off the sample. To improve the signal to noise ratio, the detector response at 574kHz is then monitored by a lock-in amplifier, locked on 574kHz, and analyzed by a computer.

III. Heat transport theoretical modeling

In a TDTR technique, the determination of the thermal conductivity and interfacial thermal resistance is made by comparing the experimental cooling curve to the theoretical model and optimizing the free parameters. For all thin metal films, the temperature is homogeneous in the zone delimited by the optical penetration depth. The surface temperature variation ∆T f (t) can be calculated using one-dimensional heat flow since the thermal diffusion length is usually much smaller than the radius of the pump beam spot. Figure 2 shows a schematic diagram (not at scale) of the studied structure with a spatial coordinate in the cross plan direction normal to the surface. Because there are two different Al film thicknesses, theoretical modeling of heat will be different from one sample to the other.

Thin metallic film on SL

In this case, the Al film is only 12 nm thick, which is about twice the optical penetration depth ξ=7nm at the wavelength of the experiment λ=780nm. After absorption of the laser pulse, all the film will be heated at the same time, and the temperature distribution will be uniform. Fourier classical heat diffusion formalism gives the following expression of the cooling of the transducer in Laplace domain:

( ) ( ) 0 0 T = T (1) τ τ θ τ τ τ τ τ τ τ τ τ τ + = ∆ + + + ∆ + + sf sf f sf sf sf sf R sf sf sf R sf bi y p p bi p bi y p p p p p where ( ) ( ) 2 2
, , , The application of Laplace inversion theorem [START_REF] Carslaw | Conduction of Heat in Solids[END_REF] gets a compact formula of the Al film transducer cooling in the time space:

f f sf s sf R s k s f c d d y y bi c bi R y ρ τ τ τ ρ β α ⊥ ⊥ = = = = and ( ) ( ) 0 1 ρ - ∆ = f f R Q T c
( ) ( ) ( ) 2 2 2 2 0 2 2 1 u t sf f R sf T e T t du u u τ π τ τ ∞ - ∆ ∆ = - + ∫ 2

. Thick metallic film on SL

In this case, the SL is covered by a 100 nm Al film and diffusion of heat into the transducer can not be neglected. Such as in the 12 nm Al case, theoretical modeling of heat propagation into the structure after the laser pump pulse is absorbed, is based on solving heat Fourier equations in Laplace domain. We have used Thermal Quadrupoles Method (TQM) [START_REF] Maillet | Thermal Quadripoles: Solving of Heat Equation through Integral Transforms[END_REF], where we have done these three assumptions: (i) The characteristic time of macroscopic thermal phenomena is much greater than the laser pulse duration (about 100 fs), then the heat flux due to the laser pulse is considered to be a Dirac function of amplitude unity : ( ) ( )

t t ϕ δ = , (ii)
Due to the fact that the optical penetration depth in Aluminum (ξ=7nm) is small compared to the thickness of the film, the latter is supposed opaque, and we neglect the penetration of the light energy into it, (iii) The hypothesis of small thermal diffusivity of the SL allows us to consider the latter as a semi-infinite medium: neither buffer-layer, nor substrate will influence the temperature variation on Al film top surface during the time scale of the experiment. After some algebra, we find this expression of the temperature variation in Laplace domain:

( ) θ β - ⊥ -   + -     =   --     2 2 2 1 1 1 3 2 1 1 f f qd f qd f e E Sq e E
, where 1 , and

f k f s f p E R q q ε β ε α ⊥ ⊥ ⊥ ⊥ = + + = .
, ,

α β ε ⊥ ⊥ ⊥ f f
f are respectively the cross-plan thermal diffusivity, conductivity, and effusivity of the Al film, s ε ⊥ is the thermal effusivity of the SL.

In the approximation of short time scale 0 t → , (3) can be expanded with respect to

2 2 1 f qd e E -   Γ = -    
, the inverse Laplace transform is then applied term by term, and we get these two expressions for the first order ( )

, FO f T d t -
and the second order ( )

, SO f T d t -
of this expansion respectively in the time domain:

( ) ( ) ( ) 2 2 2 1 , 1 2
4 2 f f f f d t f FO f f hd th f f f f T d t e t S d e Erfc h t t α α α π β ψα α α ⊥ ⊥ ⊥ - ⊥ ⊥ + ⊥ ⊥     - = + -           +     and ( ) ( ) ( ) ( ) 2 2 4 4 2 2 , 4 2
2 5 α α ψα α π β β ψα β α β α α ⊥ ⊥ ⊥ ⊥ - ⊥ ⊥ ⊥ ⊥ ⊥ ⊥ + ⊥ ⊥   - = + -       + + ×     +     f f f f d t f f SO f f k f f k f f f k f hd th f f f t T d t e t S R R d h t R d e Erfc h t t where 1 f s k f h R ε ε β ⊥ ⊥ ⊥ + =
, and

2 8 k f R S ψ β ⊥ =
, Erfc is the complementary error function. Then, the Al film top surface temperature variation, corrected to the second order in Γ , is given by:

( ) ( ) ( ) ( ) , , , 6
FO SO f f f f T d t T d t T d t ∆ - = - + - .
We must note here that this expansion is valid for the case Γ<<1, which is equivalent in Laplace space We have studied the sensitivity of the temperature variation on Al film top surface ∆T f to both Kapitza resistance R k at the interface Al film/Si/SiGe SL, and thermal diffusivity s α ⊥ of the SL, using numerical Laplace inversion [START_REF] Maillet | Thermal Quadripoles: Solving of Heat Equation through Integral Transforms[END_REF]. This study has shown that for the 100 nm Al film, the temperature variation is more sensitive to R k than to s α ⊥ . But for the 12 nm Al film, the cooling is more sensitive to s α ⊥ than to R k . In the latter case, numerical calculation has shown a very small value of R k . We can then consider a perfect thermal interface between the Al film and the Si/SiGe SL, in this case 0 R τ → , and equation ( 2) becomes:

( ) ( )

exp 7 f sf sf t t T t T Erfc τ τ       ∆ = ∆           .
Since we measure the reflectivity change instead of the temperature variation, we must now link these two quantities, but it is well known [START_REF] Capinsky | Thermalconductivity measurements of GaAs/AlAs superlattices using a picosecond optical pump-andprobe technique[END_REF] that for small temperature range, the change in reflectivity on the surface of most metals is proportional to the temperature variation:

( ) ( ) ( ) , , 8
f f R d t T d t R κ ∆ - = ∆ -
where κ is the thermoreflectance coefficient.

IV. Results and discussion

Figure 3 shows the relative optical reflectivity variation on 100nm, and 12nm Al films as a function of time. We should note here that we have used normalized curves to be able to compare experimental and theoretical curves. Experimental curves depend on experimental conditions via some parameters that we cannot include in the modeling, for instance, the thermoreflectance coefficient between temperature variations and relative reflectivity variations. The absolute change in the relative reflectivity depends on the pump-power, but in our case; it was in the range 10 -4 -6×10 -3 . The curves show a thermal decay on which are superimposed some spikes. These spikes arise from acoustic echoes bouncing back and forth inside the structure. The propagation of the elastic strain inside a layer alters the optical properties of this layer, in particular those of the transducer; these changes appear by consequence in the measured reflectivity signal when the elastic strain reaches the free surface and gives an echo. For the 100 nm Al curve, there is a fast thermal decay characteristic of the diffusion of heat inside the Al film, followed by a slow decay, characteristic of the diffusion of heat inside the SL. In the opposite, the 12 nm Al curve shows only one thermal decay which describes diffusion of heat inside the SL. In addition to this echo, we can also see an other echo, at t=327ps, which is the signature of the buffer layer.

Figure 4 shows comparison between experimental thermoreflectance signal and theoretical modeling of the 12 nm Al film case, where we have used equation [START_REF] Huxtable | Thermal conductivity of Si/SiGe and SiGe/SiGe superlattices[END_REF] to fit measurements using the value of s α ⊥ obtained by Least square Method (LSM). We founded [START_REF] Ezzahri | Study of Thermomechanical Properties of Si/SiGe Superlattices using Femtosecond Transient Thermoreflectance Technique[END_REF]

β α ρ ⊥ ⊥ - - = ×
, the specific heat s c of the SL can be found using the same formula [START_REF] Capinsky | Thermalconductivity measurements of GaAs/AlAs superlattices using a picosecond optical pump-andprobe technique[END_REF] by replacing ρ j by c j . The value of the crossplan thermal conductivity found is in good agreement with that found by Huxtable et al [START_REF] Huxtable | The effect of defects and acoustic impedance mismatch on heat conduction SiGe based superlattices[END_REF] using 3ω method [START_REF] Cahill | Thermal conductivity measurement from 30 to 750K: the 3ω method[END_REF].

We should note here that the small discrepancy between theoretical curve and experimental signal can be explained by the fact that a quantity of light energy crosses the Al film and reaches directly the SL [START_REF] Ezzahri | Study of Thermomechanical Properties of Si/SiGe Superlattices using Femtosecond Transient Thermoreflectance Technique[END_REF], which is not taken into account in the model discussed above.

Figure 4: Comparison between experimental thermoreflectance signal (solid-line), and theoretical modeling (dashed-line) for the best fit in the case of a 12nm Al transducer, on a Si/Si 0.75 Ge 0.25 superlattice

Thick metallic film on SL

Similarly, we start analyzing the acoustic contribution. In this case, from the first echo arrival time (34 ps) and the knowledge of the bulk sound velocity in Al film (v Al ≈6400 m/s), we can extract the Al film thickness. We found d f =108 nm at the impact point which is in good agreement with the assumed value of 100 nm. From the second echo arrival time and the thickness of the SL, we can extract the sound velocity inside it, we find Figure 5 shows comparison between experimental thermoreflectance signal and theoretical modeling of the 100 nm Al film case, where we have used both numerical Laplace inversion transform (dashed line), and analytical inversion based on the approximation at short time scale (6) (dotted line). As we can see the two curves are quite superimposed (curves 2 and 3) at short time, and then the analytical solution diverges. This is the reason why we have used numerical inversion to fit experimental result using LSM. We see that we have a good agreement for long time greater than 100 ps. The parameters used for the best fit are: 108

f d nm =
, and Let us note that the value of cross-plan thermal diffusivity of Al film is very high compared to that of bulk material. In fact we should have a small one due to the small thickness of the metallic film, but this value could be understood if we take into account the low specific heat of electrons, and their large thermal penetration at the first tens picoseconds [START_REF] Tas | Electron diffusion in metals studied by picosecond ultrasonics[END_REF]. Due to the fact that the Al film surface temperature variation is more sensitive to R k than to s α ⊥ , we have used for s α ⊥ the value obtained from the case of 12nm Al film above.

Although the compositions of the two SL are different, we make the assumption that the slight difference between the SL does not cause a large difference between their thermal conductivities. Using all these parameters in conjunction with the LSM allows us to get the following value for Kapitza resistance: R k ≈2.25×10 -8 m²K/W. This value is in reasonable agreement with that found in literature [START_REF] Smith | Measurement of Thermal Boundary Resistance between Thin Metal Films and Dielectric Substrates[END_REF][START_REF] Stoner | Kapitza conductance and heat flow between solids at temperatures from 50 to 300K[END_REF]. We should note here that the Kapitza resistance is highly dependent on several parameters, as it was intensively discussed by R. J. Stoner and H. J. Maris [START_REF] Stoner | Kapitza conductance and heat flow between solids at temperatures from 50 to 300K[END_REF]. It depends on: the difference of Debye temperatures between the two solids, the nature of the heat transfer process across the interface, in such a way that it can involve only phonons or both phonons and electrons, the anharmonicity of lattice potential, and the nature of the diffusion process at the interface, in such a way that it can be elastic or inelastic. In The discrepancy observed for the first 100ps, could be explained by the large electron thermal penetration. When the pump light is absorbed, its energy is at first communicated to electrons near the front surface of Al film which is excited to higher energy states. These hot electrons quickly diffuse away from the Al surface, but are confined to the Al film by the Schottky barrier at the interface Al/SL. Within several picoseconds, the hot electrons transfer their energy to the lattice by emission of phonons or electron-phonon collisions, slightly raising the temperature of the Al film. The depth over which electrons are able to diffuse before completely transferring their energy to phonons bath has been estimated by G. Tas and H. J. Maris to be approximately 100nm [START_REF] Tas | Electron diffusion in metals studied by picosecond ultrasonics[END_REF]. By contrast other works [START_REF] Qiu | Short-Pulse laser heating on metals[END_REF][START_REF] Qiu | Heat Transfer Mechanisms During Short-Pulse Laser Heating of Metals[END_REF] have shown that time over which the electron gas will be completely thermalized within the phonons bath is about 5ps. After this time the Al film can be characterized by one temperature which is that of the lattice. Then Fourier heat classical equation can properly be applied.

Another assumption to explain the discrepancy is the effect of the pump power. In fact, in the majority of the models used in TDTR to extract thermophysical properties, we assume that the properties of the structure are temperature independent. The high value of the pump peak power (about 10 kW) in conjunction with the slow cycle rate due to the modulation of the pump beam could modify the film properties via an increase of its temperature. To check this effect, we have done some measurements on the same structure following the same experimental process and varying the pump power. The results are shown in figure 6. As we can see, the influence of the pump power does not seem significant as the thermal decay does not vary much. We are making new experiments to be sure that the thermal decay does not depend on pump power.

We are now working to make more experimental measurements to explain the discrepancy on the first 100ps and to develop an accurate modeling of the energy transport in thick metallic film at short time scale which must include microscopic interactions between electrons themselves and between the phonons bath, as well as the temperature dependence of the physical properties. In this case, a possibility to access absolute temperature measurement on the sample will be more difficult since the relation between reflectivity change and temperature change will be more complicated.

V. Conclusion

In this paper, two Si/SiGe superlattices structure, with slight different compositions, covered by different thickness Al transducers were studied using TDTR technique. Taken into account the short time scale of the experiment, we have shown the influence of the transducer thickness on the sensitivity of the measurement to the thermophysical properties of the underlying superlattices. By this way we were able to measure both cross-plan thermal conductivity of the SL and the Kapitza resistance at the interface transducer/SL. We were also able to estimate the sound velocity inside the Si/SiGe SL from the arrival time of acoustic echoes. This technique therefore proved to be a useful tool to characterize both thermophysical and mechanical properties of thin film multilayers. The preliminary results show that the discrepancy between experimental and theoretical curves observed in the case of thicker transducer, could not be explained by the effect of the pump power. Diffusion of electrons and then the penetration of the heat source seems to be responsible for this behavior. A more precise model which takes into account this effect should be the convenient solution.

Figure 1 :

 1 Figure 1: Experimental set up of the TDTR technique.

Figure 2 :

 2 Figure 2: Schematic diagram of the whole studied structure.

  Sd where R is the reflection coefficient of Al film, Q is the pump pulse energy, d f is the Al film thickness, and S is the illuminated area of the Al film surface. (ρC) S , and (ρC) f are the specific heat per unit volume of Si/SiGe SL and Al film respectively. s β ⊥ and R K are the cross-plan thermal conductivity of the SL and Kapitza resistance at the interface Al/SL respectively.

Figure 3 :

 3 Figure 3: Experimental thermoreflectance signal of 100nm, and 12 nm Al film on Si/SiGe superlattice.

  Then we can compare this measured value to the one estimated using harmonic average[START_REF] Rossignol | Elastic properties of ultrathin permalloy/alumina multilayer films using picosecond ultrasonics and Brillouin light scattering[END_REF] this case the estimated value is about 13% less than the measured value.

  the analytical solution diverges after ~140 ps as the validity condition (

Figure 5 :

 5 Figure 5: Comparison between experimental thermoreflectance signal and theoretical modeling for the best fit in the case of a 100nm Al transducer, on a Si/Si 0.7 Ge 0.3 superlattice.

Figure 6 :

 6 Figure 6: Experimental thermoreflectance signal of 100nm Al film on Si/SiGe superlattice with different pump powers.

table 1 ,

 1 are listed all identified properties.

	Properties	100nm	12nm
		thick	thick
		transducer	transducer
	Transducer thickness	108	12
	(nm)
	Thermal conductivity	8.5*	8.5
	(W/m/K)
	Thermal boundary	2.25×10 -8	//
	resistance (m².K/W)
	SL sound velocity	~9000	~7782
	(m/s)

Table 1 :

 1 Different

properties of the two samples identified by comparison of experimental and theoretical curves, * not identified.