Lucas Duperrex
email: lucas.duperrex@lynred.com

Raphaël Pesci

Pascal Le Boterf

Olivier Mailliart

Simulation and measurement of residual stress and warpage in a HgCdTe-based infrared detector at 100 K

published or not. The documents may come

Simulation and measurement of residual stress and warpage in a HgCdTe-based infrared detector at 100 K Lucas Duperrex, Raphaël Pesci, Pascal Le Boterf, Olivier

Introduction

Hg 1-x Cd x Te is the most used material to manufacture photodiode arrays for high performance InfraRed (IR) detection. The semiconductor properties lead to the best photodiodes [START_REF] Rogalski | Type-II superlattice photodetectors versus HgCdTe photodiodes[END_REF][START_REF] Kinch | The future of infrared; III-Vs or HgCdTe ?[END_REF] and it can be adapted to every IR wavelength by adjusting the Hg/Cd ratio. It is part of II-VI materials, nicknamed "soft brittle" [START_REF] Irwan | Mechanical properties and material removal characteristics of soft-brittle HgCdTe single crystals[END_REF][START_REF] Schlesinger | Cadmium zinc telluride and its use as a nuclear radiation detector material[END_REF] with low hardness and high brittleness compared to classical semiconductor materials such as silicon. However, the poor mechanical properties can lead to defects in some photodiodes and the emergence of dead pixels. Despite decades of research and significant progress, this issue still hasn't been entirely solved.

The elaboration of HgCdTe IR detectors can be divided into four main steps: epitaxy of thin Hg 1-x Cd x Te layer on lattice-matched Cd 1- y Zn y Te, manufacturing of photodiode array, hybridization (flip-chip bonding) to silicon-based Read Out Integrated Circuit (ROIC) and integration in cryostat (after CdZnTe removal). This work focuses on the last two steps' contribution to the level of residual stress in HgCdTe and to warpage, in view of the optical planarity requirements. The key features are the thermal expansion mismatch between HgCdTe and silicon and the operating temperature of the detector around 100 K. They lead to thermal bending and biaxial tension in the HgCdTe layer at low temperature. Since the 90s, many patents have described solutions to reduce the level of stress, based on compensation layers added to the backside of the silicon substrate of the ROIC in order to obtain a coefficient of thermal expansion (CTE) close to that of HgCdTe. However, the CTE mismatch between HgCdTe and silicon has proven to be acceptable if the epitaxial substrate (CdZnTe) is removed (including for very large formats such as 4 k*4 k, 15 μm pitch array [START_REF] Blank | Development and production of the H4RG-15 focal plane array[END_REF]). Nevertheless, thermomechanical studies of IR FPAs are still ongoing. Among the expected improvements are the enlargement of detectors format, the reduction of the pixel pitch and of the stress-induced defects in HgCdTe in order to reach higher operating temperatures. In other words, the limits of HgCdTe are far from being reached, as confirmed by recent results [START_REF] Lee | High-operating temperature HgCdTe: a vision for the near future[END_REF].

Some XRD stress measurements in HgCdTe have already been published in literature, considering a simple epitaxial layer on a CdZnTe substrate [START_REF] Colin | Elastic and plastic deformation in low mismatched CdxHg1-xTe/Cd1-yZnyTe[END_REF][START_REF] Ballet | Strain determination in quasi-lattice-matched LWIR HgCdTe/CdZnTe layers[END_REF] and an achieved IR detector [START_REF] Lebaudy | X-ray diffraction residual stress measurement at room temperature and 77 K in a microelectronic multi-layered single-crystal structure used for infrared detection[END_REF][START_REF] Lebaudy | Multilayer CdHgTe-based infrared detector: 2D/3D microtomography, synchrotron emission and finite element modelling with stress distribution at room temperature and 100 K[END_REF]. For the epitaxy process, biaxial stress values were measured, mainly induced by imperfect lattice match between Hg 1-x Cd x Te and Cd 1-y Zn y Te. Considering the IR detector, biaxial stresses (compressive at room temperature and tensile at low temperature) were also highlighted and the corresponding values were in good agreement with the above-mentioned difference in CTE between HgCdTe and silicon. A gradient between the center and the edges of the samples was reported [START_REF] Lebaudy | X-ray diffraction residual stress measurement at room temperature and 77 K in a microelectronic multi-layered single-crystal structure used for infrared detection[END_REF], but the origin remains to be clearly identified, even with Finite Element Method (FEM) [START_REF] Lebaudy | Multilayer CdHgTe-based infrared detector: 2D/3D microtomography, synchrotron emission and finite element modelling with stress distribution at room temperature and 100 K[END_REF]. The warpage evolution with temperature for such detectors does not seem to have ever been studied before, but numerical simulations have been performed on InSb-based detectors [START_REF] Zhang | Modeling and deformation analyzing of InSb focal plane arrays detector under thermal shock[END_REF][START_REF] Zhang | Assessment of InSb infrared detector arrays assembly procedure employing ANSYS[END_REF] without any measurements at low temperature. Thus, the literature only gives an incomplete description of the thermomechanical behavior of HgCdTe-based detectors when cooling to operating temperature.

Thermomechanical issues in such detectors also exist at the micrometric scale, such as stress gradient in HgCdTe at the pixel scale [START_REF] Lebaudy | Multilayer CdHgTe-based infrared detector: 2D/3D microtomography, synchrotron emission and finite element modelling with stress distribution at room temperature and 100 K[END_REF] and other phenomena which impact reliability of interconnections in all flip-chip devices. From what we know, both scales are important, but neither experimental nor simulation methods enable to study both of them at the same time. It is necessary to study, quantify, and clearly understand them separately, before being able to take them into account simultaneously (and maybe finally identify the most critical one). This work deals with strain and stress variations at the millimeter scale and is thus a step towards a complete (multi-scale) understanding of the thermomechanical phenomena that take place in such detectors. It is based on a multi-method thermomechanical analysis of a 320 × 256, 30 μm pitch, middle wave IR FPA, for both the detector on its own and after assembled to the hosting substrate. It includes planarity measurements, XRD stress analyses in HgCdTe ((111) single crystal), an analytical approach based on Stoney [START_REF] Stoney | The tension of metallic films deposited by electrolysis[END_REF] and Timoshenko's works [START_REF] Timoshenko | Analysis of bi-metal thermostats[END_REF] as well as FEM, at both room temperature and 100 K. Finally, after confirming the relevance of the FEM in accordance with the experimental measurements, larger detector formats are extrapolated in order to predict the evolution of warpage and residual stress in HgCdTe.

Experimental procedures and results

Presentation of the detector architecture

The production of the IR detectors studied in this work starts with a 5-to-10 μm thick HgCdTe layer obtained by liquid phase epitaxy on a quasi-lattice-matched CdZnTe substrate (both are single crystals). For low mismatch (under elastic assumption), it can be considered that the lattice parameter a of the substrate is imposed to the epitaxial layer in the in-plane directions [START_REF] Colin | Elastic and plastic deformation in low mismatched CdxHg1-xTe/Cd1-yZnyTe[END_REF][START_REF] Ballet | Strain determination in quasi-lattice-matched LWIR HgCdTe/CdZnTe layers[END_REF]. Then, the gap between lattice parameters gives directly biaxial strain values:

ε ‖ = ( a substrate -a epitaxial layer ) / a epitaxial layer (1)
Significant lattice mismatches have deliberately been chosen in this work to show that this well-known phenomenon still has an impact on residual stress after substrate removal.

Photodiodes arrays are manufactured collectively at the free surface of HgCdTe, which is then sliced into single chips and hybridized (indium solder bump flip-chip bonding) to the ROIC. The latter, about 7 μm thick, consists in metal circuits in a SiO 2 -based matrix above the silicon bulk and is generally called "CMOS" in reference to the name of the transistor fabrication process (Complementary Metal-Oxide-Semiconductor).

Afterwards, an epoxy resin is introduced by capillary underfill process. The curing temperature corresponding to the samples studied in this work is 363 K. The resulting interconnection layer (indium bumps and epoxy) is about 8.5 μm thick (see zoomed view in Fig. 1). From there, HgCdTe is mechanically assembled to the ROIC and the CdZnTe substrate can be removed by both mechanical and chemical thinning. A little oversupply of epoxy leads to an underfill filet which is in major part thinned with the substrate (see Fig. 6 in the FEM section).

Excluding the sides, the assembly can be seen as a multilayer plate:

HgCdTe 5-10 μm/interconnection layer 8.5 μm/CMOS 7 μm/silicon 500 μm. The "central part" corresponds to the 320 × 256, 30 μm pitch HgCdTe photodiodes array, while silicon and CMOS complete size is 10,450 μm*12,950 μm (Fig. 1). This is the "hybridized detector state", which is studied before assembly to the hosting substrate.

After performance tests, this hybridized circuit finally needs to be bonded to a hosting substrate, which has a role of thermal and mechanical link with the cryostat. In this work, in order to limit stress and warpage, aluminum nitride has been chosen for its good thermal expansion match with silicon (other possibilities exist such as Invar or SiC). Glue with a low rigidity was used for this assembly and cured at 353 K; its thickness is approximately 20 μm.

XRD stress analyses were performed on three of these detectors (samples A, B and C: see section 2.2), warpage measurements on two (samples D and E: see section 2.3). Sample A, B and C has been chosen with HgCdTe layers coming from three different CdZnTe substrates with increasing lattice parameter to show the influence of lattice mismatch on stress results. Moreover, sample C includes a thicker silicon substrate and a higher pixel density (15 μm pitch) than others. Thus, the properties of the CMOS and interconnection layers are a little different, but, as will be shown in the modelling section, they do not affect significantly the stress level in HgCdTe; neither does the substrate thickness.

X-ray diffraction at room temperature and 100 K

The experimental procedure and the XRD device are described in detail in Ref. [START_REF] Lebaudy | X-ray diffraction residual stress measurement at room temperature and 77 K in a microelectronic multi-layered single-crystal structure used for infrared detection[END_REF]. Ortner's method [15] is used to quantify the level of residual stress in the HgCdTe single crystal at room temperature and 100 K. This method is based on the measurements of the inter-reticular distances d hkl considering Bragg's equation: λ = 2. d hkl .sin(θ), where λ is the X-ray wavelength and θ the angle between the specimen and the incident beam. All the measurements are made on {533} planes using a FeK α = 0.1936 nm radiation (corresponding diffraction peaks at 2θ = 158 • ); a collimator 1 mm in diameter is used to obtain sufficient intensity and great peak definition. First, the single crystal orientation is determined through a pole figure based on Bragg's equation. Then, the position of each diffracting plane can be simulated through 2 angles (azimuth φ and tilt ψ). All the peaks are then acquired in a 2θ range in order to determine precisely their inter-reticular distances with Bragg's law (at least six are necessary). The strain tensor is finally calculated and Hooke's law leads to the stress tensor considering the elastic constants of HgCdTe (C 11 = 53.5 GPa, C 12 = 36.8 GPa, C 44 = 19.9 GPa [START_REF] Capper | Mercury Cadmium Telluride -Growth, Properties and Applications[END_REF]).

The results are given for three samples at hybridized detector state; for one of them, the residual stress was measured before and after gluing to the hosting substrate (Table 1). Sample A was analyzed both at the center and in a corner. At room temperature, the stress is slightly compressive and close to biaxial: σ xx is a few MPa higher than σ yy, and shear stress values are very low. There is no significant difference between center and corner. At low temperature, HgCdTe is submitted to biaxial tension around 35 MPa.

For sample B, whose HgCdTe was obtained from a CdZnTe substrate with a lower Zn composition, stress was only measured at the center, showing values close to zero at room temperature and biaxial tension of 46 MPa at 100 K.

Sample C has been obtained from a CdZnTe with an even lower Zn composition. As expected with low lattice mismatch [START_REF] Colin | Elastic and plastic deformation in low mismatched CdxHg1-xTe/Cd1-yZnyTe[END_REF][START_REF] Ballet | Strain determination in quasi-lattice-matched LWIR HgCdTe/CdZnTe layers[END_REF], the lower the Zn composition, the higher the tension values. This sample is all the more interesting that these high values at 293 K (12 MPa) lead to 57 MPa at 100 K. This will be of great importance for the upcoming discussion about the yield stress of HgCdTe.

For these three hybridized detectors with no hosting substrate, considering both σ xx and σ yy , the stress difference between 293 K and 100 K is almost the same: 40-45 MPa. The significant biaxial tension values induced while cooling to 100 K are consistent with a previous study [START_REF] Lebaudy | X-ray diffraction residual stress measurement at room temperature and 77 K in a microelectronic multi-layered single-crystal structure used for infrared detection[END_REF]. Moreover, gluing then sample A to a hosting substrate ("S. A*" in Table 1) leads to an increase of biaxial compression values at room temperature and biaxial tension values of the same order when cooling down to 100 K.

Warpage evolution between room temperature and 100 K

Warpage measurements were performed with laser scanning equipment OGP FLASH CNC 300 through the window of a specially-designed cryostat. Raw data consisted in (x,y,z) maps of approximately three thousand points. Random planarity defects are induced by manufacturing process and are not negligible in comparison to variations with temperature so that a variability of the initial warpage is expected. However, the aim is to characterize the warpage evolution with temperature. Consequently, perturbed surfaces have been approximated with paraboloids: z = f(x,y) = c 1 + c 2 x + c 3 y + c 4 (x 2 + y 2 ) (approximation by least square method, an example of a result is illustrated in Fig. 2b). Measured warpage Wm can then be determined using two methods leading to the same results. One of them consists in correcting tilt and considering gap between maximum and minimum values, the other in applying the above-mentioned function f(x,y) in the following equation:

Wm = (d/2) 2 / 2R with R = ( 1 + (∂f /∂x) 2 ) 3/2 /( ∂ 2 f / ∂x 2 ) (2)
The measured surface is defined manually: it is a bit smaller than that considered with FEM because in practice, to consider only the HgCdTe surface, it is necessary to leave a 0.5 mm gap at the edge. In order to be able to compare measured warpage to FEM results, the measured values have been adjusted assuming constant curvature radius. This implies that warpage is proportional to squared diagonal. With dm the diagonal

Table 1

XRD residual stresses with corresponding uncertainties in the HgCdTe layer of IR detectors at 293 K and 100 K (all values in MPa). (*) indicates data obtained after gluing to hosting substrate.

Sample

Temp.

σ xx σ yy σ zz σ xy σ xz σ yz Uncertainty S. A center 293 K -6 -3 0 -1 -1 0 <3 S. A center 100 K 31 38 0 3 1 -2 <3 S. A corner 293 K -9 -5 0 0 0 0 <4 S. A corner 100 K 34 34 0 0 3 0 <5 S. A* center 293 K -12 -14 0 1 0 1 <4 S. A* center 100 K 35 39 0 -1 1 0 <5 S. B center 293 K 1 1 0 1 0 0 <3 S. B center 100 K 45 47 0 -2 1 -2 <4 S. C center 293 K 12 12 0 0 -1 1 <3 S. C center 100 K 56 57 0 -1 -1 0 <5
of the measured map and dn the nominal diagonal (used in FEM), the finally considered warpage W is calculated as follows:

W = Wm*dn 2 / dm 2 (3)
In this work, positive and negative values are attributed to concave and convex samples, respectively (Fig. 2c). Note that when studying thermal bending, it would be more relevant to talk about curvature radius, which is independent of the considered detector format, instead of warpage (deflection). However, this notion is often used in microelectronics, probably because it can be easily linked to pick-to-valley and flatness requirements (either regarding flip-chip assembly feasibility or optical domain).

The cryostat used to reach low temperatures consists in a nitrogen bath in contact with a cold table in a chamber under vacuum (Fig. 2a). First, the upper part of the chamber is opened and the sample is glued in a corner, "locally", in order to ensure thermal transfer without disturbing self-warping. A PT1000 temperature sensor is also glued on the holder right next to the sample, and the chamber is closed. Vacuum is then applied through a dynamic pump for a few hours in order to degas all materials (glue, polymers of connection cables); it is maintained by a static getter to avoid vibrations during measurement. Temperature is measured but not controlled; the slow warming up to room temperature (associated to progressive evaporation of liquid nitrogen) is exploited to achieve intermediary measures between 85 and 90 K and 293 K. Hot water can replace liquid nitrogen to reach high temperatures until around 323 K.

The thermal homogeneity and the reliability of the measurements were verified through a second PT1000 sensor glued on the sample: the maximum gap between the two sensors (at low temperature) was only 3 K, with an additional measurement noise of about ± 1 K. Consequently, the uncertainty on temperature values can be considered less than 5 K.

For hybridized circuits, the warpage of both samples was about 0.5 μm at room temperature and around 2.5 μm at 100 K (more concave) (Fig. 3). After gluing to hosting substrate, absolute values are shifted by about -1 μm (less concave): the evolution with temperature is then significantly lower.

Analytic and finite-element modelling

The aim was to model the detector at the sample scale with reasonable computing performances which means that the geometry of each indium bump and CMOS metal circuit could not be considered. It could be possible with 2D models, but the assumptions associated lead to errors for flip-chip packages as it has been illustrated with an example by Yao et al. [START_REF] Yao | Three-dimensional versus two-dimensional finite element modeling of flip-chip packages[END_REF]. Therefore the use of Equivalent Homogeneous Materials (EHM) is required [START_REF] Kpobie | Three-Dimensional thermomechanical simulation of fine-pitch high-count ball grid array flip-chip assemblies[END_REF] for the CMOS and the interconnection layer considering a one-factor-at-a-time method to take into account the uncertainties due to local variations of the chemical composition.

In this work, the main part of the hybridized circuit (excluding sides) can be assimilated to a multi-layer stack, including one silicon layer much thicker than the others. HgCdTe and silicon layers are single crystals (111) and (100)-oriented, respectively. For these orientations and thermal bending, it has been shown that single crystals with cubic structure can be considered as isotropic [START_REF] Brantley | Calculated elastic constants for stress problems associated with semiconductor devices[END_REF][START_REF] Wortman | Young's modulus, shear modulus, and Poisson's ratio in silicon and germanium[END_REF] (invariance of biaxial modulus in these planes). Therefore only isotropic properties are presented here.

Material properties

The properties of silicon, HgCdTe and AlN are given in Table 2 with an overview of thermal expansion values whose variations with temperature are not linear (discrete values every 10 K were considered). EHM of CMOS was calculated with a simple mixture law:

-classic Voigt model [START_REF] Voigt | Ueber die Beziehung zwischen den beiden Elasticitätsconstanten isotroper Körper[END_REF] for elastic properties.

-for the CTE, starting from the same geometric assumption as Voigt, the calculation consists in an average which is weighted with both volume and Young's modulus.

Results with this EHM assumption are then compared to assimilation of this layer to pure silicon or to SiO 2 in the local sensitivity analysis section.

Silicon, AlN and CMOS are considered linear elastic (their yield stress values are not reached). For HgCdTe, which is considered elasticperfectly plastic, based on the higher value obtained by XRD analysis including uncertainty, the yield stress value considered in the model is 65 MPa and will be discussed in the last part.

The properties of the underfilling epoxy and the glue used for assembly are presented in Fig. 4. The yield stress of epoxy is about 60 MPa at room temperature and the evolution with temperature presented in Ref. [START_REF] Fiedler | Influence of stress state and temperature on the strength of epoxy resins[END_REF] was taken into account. Concerning the glue, the yield stress is about 5 MPa for temperatures higher than 170 K (glass transition temperature) and reaches 50 MPa at 100 K.

The properties of the interconnection layer (indium bumps in epoxy matrix) are obtained considering an EHM: they vary significantly with temperature. Multi-linear stress-strain curves were calculated at several temperatures (Voigt assumption) and are presented in Fig. 5. Each of them includes a first segment for which both indium and epoxy have elastic behavior. The yield stress of indium being about 1 MPa at room temperature [START_REF] Reed | Tensile strength and ductility of indium[END_REF][START_REF] Luo | Unusual tensile behaviour of fibre-reinforced indium matrix composite and its in-situ TEM straining observation[END_REF], this behavior at low strain needs a subplot to be visible. Then there is a second segment for which epoxy is still elastic and indium has plastified, and a last horizontal segment for which both are plastic.

Epoxy resins are generally initially strained in tension by the shrinkage due to curing (not the one due to thermal expansion), the main part of which occurs in liquid phase, but it has been shown that it can lead to significant warpage [START_REF] Lim | Impact of mold compound cure shrinkage on substrate block warpage simulation[END_REF] and residual stress [START_REF] Schindler-Saefkow | Measurements of the mechanical stress induced in flip chip dies by the underfill and simulation of the underlying phenomena of thermalmechanical and chemical reactions[END_REF] in the case of flip-chip assembly. The "mechanical shrinkage" and the resulting residual stress values are dependent on the curing conditions (temperature, pressure, moisture and containment): they are therefore very difficult to estimate. They are neglected in nominal simulations but the impact of mechanical shrinkage values of 0.22% [START_REF] Lim | Impact of mold compound cure shrinkage on substrate block warpage simulation[END_REF] and 1% (which seems to be the maximum expectable value for a "low shrinkage" epoxy) will be studied in the local sensitivity analysis section. These ΔL/L values are directly added to the CTE in the case of underfill filet residue and prorated to epoxy volume fraction in the case of interconnection layer.

Note that the FEM software used, Ansys, operates calculation with linear assumption between all discrete values detailed in this section.

Analytical models

The possibility to consider HgCdTe and silicon single crystals as isotropic for thermal bending allows the use of analytical models based on Stoney formula [START_REF] Stoney | The tension of metallic films deposited by electrolysis[END_REF] or Timoshenko's work [START_REF] Timoshenko | Analysis of bi-metal thermostats[END_REF] for the calculation of the warpage of the hybridized detector. The considered geometry is a simple multilayer plate whose in-plane dimensions are limited to those of HgCdTe. Elastic assumption is also required. The loading consists only in cooling down from 363 K to 100 K. The corresponding results are interesting because they give a first estimation of the expected values and are very helpful to understand the phenomena involved.

The Stoney formula was initially established on strips, but it can be adapted to plates using the biaxial modulus rather than the Young modulus to determine the warpage W:

W = 3L 2 s 4 ⎛ ⎜ ⎜ ⎝ Ef (1-νf ) t f ΔαΔT Es (1-νs) t 2 s ⎞ ⎟ ⎟ ⎠ (4)
with L s the diagonal of the plate, E s , ν s and t s the Young's modulus, Poisson's ratio and the thickness of the substrate, E f , ν f and t f the Young's modulus, Poisson's ratio and the thickness of the thin layer, Δα = α s -α f the CTE difference between the film and the substrate, and ΔT = T final -T initial the temperature variation. Note that warpage is directly proportional to the stiffness and the thickness of the thin layer, and to the difference between its CTE and that of the substrate. At cryogenic temperatures, the use of average values for CTE does not make sense any longer because they show great variations with the considered temperature range. Therefore, the use of instantaneous CTE is more suitable, which implies in our case using an integral of (α s -α f ) over the temperature range rather than a simple Δα.ΔT mathematical product. Moreover, as shown by Stoney, the contribution of several thin layers to the warpage can be added (as long as the thickness of cumulated layers remains very thin compared to the substrate, their contribution to flexural rigidity is not significant). Taking into account these two evolutions, the following equation can be established:

W = ∑ n i=1 3L 2 s 4 ⎛ ⎜ ⎜ ⎜ ⎝ Efi (1-νfi) t fi ∫ Tfinal Tinitial ( α s -α f ) dT Es (1-νs) t 2 s ⎞ ⎟ ⎟ ⎟ ⎠ ( 5
)
where n is the number of layers.

The same adaptations can be made starting from Timoshenko's work (equations ( 6) and ( 7)), where R is the radius of curvature) and lead to equation ( 8):

W = L 2 s 8R (6) 1 R = Δα.ΔT h 2 + 2 ( Es t 3 s 12 +Ef t 3 f 12 ) h ( 1 Ests + 1 Ef tf ) (7)
W = ∑ n i=1 L 2 s 8 ∫ Tfinal Tinitial ( α s -α fi ) dT h 2 + 2 ( Es (1-νs ) t 3 s 12 + E fi ( 1-ν fi ) t 3 fi 12 ) h ⎛ ⎜ ⎜ ⎜ ⎝ 1 Es (1-νs ) ts + 1 E f ( 1-ν fi ) tfi ⎞ ⎟ ⎟ ⎟ ⎠ (8)
The two formulas lead to the same results: a warpage of 6.4 μm at 100 K including 5 μm induced between room temperature and 100 K. This value, which is significantly different from those measured, will be compared to FEM in Fig. 7. Regarding stress in thin layers, if relaxation due to bending is neglected, a first approach consists in considering that the substrate imposes its thermal expansion to the thin layer (in-plane directions), the deformation of which ε f can be calculated directly from the CTEs over the temperature range.

ε f = ∫ Tfinal Tinitial ( α s -α f ) dT (9)
The stress value is finally obtained by multiplying by the biaxial modulus, which leads in our case to:

σ f = E f ( 1 -ν f ) ∫ Tfinal Tinitial ( α s -α f ) dT ( 10
)
The obtained value for the stress evolution in HgCdTe between 363 K and 100 K is 53 MPa, with 40 MPa induced between room temperature and 100 K. This is in good agreement with XRD results (Table 1).

FEM

FEM model definition

As described before, it is assumed in this work that HgCdTe is mechanically assembled to the ROIC from the polymerization of epoxy underfill at 363 K: it corresponds to the initial state as hybridized circuit. To model the entire structure with hosting substrate, birth/death function is used. During the cooling down to 100 K of the hybridized circuit alone, the elements corresponding to the hosting substrate and the glue already exist but are deactivated. The hybridized detector is then heated to 353 K: it is thus already stressed, and the unstressed elements of the glue and the hosting substrate are activated. The entire structure is then cooled down to 100 K. Since transient thermal phenomena are not studied, thermal loadings are applied uniformly to every elements of the model.

Considering geometric conditions, only a quarter of the sample needs to be modeled. Orthogonal displacements equal to zero are imposed to the faces that correspond to symmetry planes. To prevent rigid body displacement in z direction, a fixed support is applied to one point which belongs to both symmetry planes.

The mesh is continuous, composed of 14,664 quadratic hexahedral elements (solid 186) for the hybridized circuit (Figs. 6) and 12,093 more for the hosting substrate and glue; a minimum of three elements in the thickness (z direction) of each thin layer is imposed. Convergence with the number of elements was verified in all directions regarding stress in HgCdTe (excluding edge effects, where stress concentrations exist) and warpage. Two models were tested by adjusting the size of elements either in the z direction or in both x and y directions: they gave the same results as the nominal model with a difference less than 1.5%.

FEM results

Two main simplifications exist between the analytical approach already presented and the FEM: elasticity assumption and geometry simplification. Therefore, a progressive modelling strategy has been adopted with two preliminary models before the nominal one:

-a first FEM model taking into account the same data as those considered for the analytical approach, with elastic assumption for all materials and with simplification of geometry to a multi-layer plate (HgCdTe/interconnection layer/CMOS/silicon with x and y dimensions limited to those of HgCdTe layer). -a second FEM model: still with elastic assumption but taking into account the complete geometry. -the nominal FEM: taking into account both plastic behaviors (for all materials concerned, see section 3.1), and the complete geometry (section 3.3.1).

The first case includes exactly the same assumptions for geometry and material properties as the analytic models, which means that the results are directly comparable. As shown in Fig. 7 (green curve) for warpage, it is in very good agreement with analytical models. Still under elastic assumption, but taking into account complete geometry, the second model (Fig. 7 -yellow curve) significantly reduces the warpage values. Finally, the consideration of plastic behaviors in the nominal FEM model (Fig. 7 -red curve) further reduces these values: they are in good agreement with the experimental measurements (Fig. 3a).

After assembly to hosting substrate, nominal FEM warpage values are reduced with an inflexion at around 170 K. This is still in good agreement with the measurements, as it was illustrated in Fig. 3a.

For the stress results in HgCdTe, the value obtained with the nominal FEM model is 49 MPa: it is thus less influenced by elastic assumption and geometry simplification than warpage (even equation [START_REF] Lebaudy | Multilayer CdHgTe-based infrared detector: 2D/3D microtomography, synchrotron emission and finite element modelling with stress distribution at room temperature and 100 K[END_REF] led to a value only 9% higher). As expected from the analytical models, (apart from edge effects and stress concentration), the stress in HgCdTe is biaxial and homogeneous in x and y directions (Fig. 8), as in all thin layers. For the hybridized detector at 100 K, σ xx and σ yy values are approximately 70 MPa in the interconnection layer and 30 MPa in CMOS. Stress in the silicon substrate is also biaxial (x and y directions) but the level depends on z position, with -9 MPa at the interface with the layer above and 5 MPa on the back (this z distribution is in good agreement with the analytical calculation detailed in Timoshenko's work [START_REF] Timoshenko | Analysis of bi-metal thermostats[END_REF]). The underfill filet residue is highly stressed in tension and impacts values in HgCdTe at the edges.

After gluing to hosting substrate and cooling back to 100 K, σ xx and σ yy values did not change significantly in thin layers, with 51 MPa in HgCdTe, 70 MPa in the interconnection layer and 31 MPa in CMOS. The stress in silicon is still biaxial with values depending on z position: -4 MPa at the interface with CMOS and -6 MPa at the glued interface with hosting substrate. The latter also shows a low biaxial stress with a distribution along z position: from 4 MPa at the interface with the glue to -6 MPa at the backside. Finally, the biaxial stress in the glue is about 49 MPa, close to the yield stress.

Local sensitivity analysis

In order to identify the most sensitive parameters, i.e. those which can be responsible for the difference between experimental measurements and FEM results, a local parametric study was performed considering nominal FEM through a one-factor-at-a-time method.

The studied parameters were uncertainties on CMOS layer (thermomechanical properties, initial stress and alternative assumptions: CMOS = Si and CMOS = SiO 2 ), interconnection layer (thermomechanical properties, thickness and shrinkage due to epoxy curing) and HgCdTe (thermomechanical properties and initial strain). Detailed results are presented only for the main parameters in Table 3.

One of the main tendencies highlighted by this sensitivity analysis is that the stress in HgCdTe is barely impacted by the uncertainties, with the exception of those attached to HgCdTe itself. As expected from XRD results, different initial stress states only lead to different offsets: it does not change the evolution with temperature (in Table 3, for "σ inital HgCdTe (biaxial)", stress values are presented to ease the reading, but the parameter which has been modified is the initial strain: it has been combined with the CTE).

Moreover, many parameters have a strong influence on warpage. The thickness of each thin layer impacts warpage linearly, and so do the stiffness and thermal expansion until the yield stress is reached.

Extrapolation to large formats

After complete analysis of the studied IR detector, FEM simulation has been applied to larger formats to predict the evolution of stress and warpage: 500 × 500, 1000 × 1000, 1500 × 1500 and 2000 × 2000 pixels (30 μm pixel pitch). In both cases (detector alone in Fig. 9a and detector with hosting substrate in Fig. 10a), the FEM results show that the warpage grows dramatically with the size of detectors. The analytical models indicate proportionality between warpage and squared diagonal of the sample (eqs. ( 6) and ( 7)). Figs. 9b and 10b confirm that this relation is still relevant, even without considering the elastic assumption and the geometrical simplification associated with analytical models.

On the other hand, stress in HgCdTe is not significantly impacted by the increase of the detector format (FEM). The stress increase of +40 Fig. 7. Comparison of warpage results considering the nominal FEM, the elastic FEM with complete geometry, the elastic FEM with simplified geometry and the analytical models based on Stoney's [START_REF] Stoney | The tension of metallic films deposited by electrolysis[END_REF] and Timoshenko's [START_REF] Timoshenko | Analysis of bi-metal thermostats[END_REF] works.

MPa from room temperature to 100 K remains constant whatever the format. This result is in good agreement with the analytical models (eq. ( 8)) wherein the thin layer stress calculation does not depend on format.

Discussion

Stress in HgCdTe

For the hybridized detector only, if there are some differences for the three samples (A, B and C) due to variations of the lattice-mismatchinduced stress from the epitaxial growth, all the values measured by XRD shown in Table 1 are biaxial in the xy plane and consistent with the analytical and finite element models (Fig. 8). The variations as a function of temperature are similar: approximately +40-45 MPa from room temperature up to 100 K. The biaxiality can be explained as a first approximation by the fact that the substrate imposes its CTE on thin layers in all the in-plane directions (eq. ( 7)), while the orthogonal direction (z) is stress-free (free surface). As regards the lattice parameter variations coming from epitaxy, the expected results are confirmed (eq. ( 1)), with the most compressive stress for the sample A whose substrate CZT had the lowest lattice parameter, and the highest tensile stress for the sample C whose CZT substrate had the highest lattice parameter.

Sample C shows other differences

-it first has a slightly thicker silicon substrate, but as illustrated by equation 10 and the work of Stoney, as long as this thickness is very large compared to that of thin films, the residual stresses are not significantly impacted. -the CMOS and interconnection layers have higher pixel densities, and therefore slightly different thermomechanical properties from those of samples A and B. But the sensitivity analysis showed that slight variations in the properties of these thin layers do not significantly impact the stresses in HgCdTe.

This sample is also particularly interesting when discussing the yield stress of HgCdTe, since it is the one which leads to the highest tensile values: 57 MPa at 100 K. Since the behavior of HgCdTe is expected to be elastic-perfectly plastic, the yield stress is therefore at least equal to this value. It is significantly higher than those mentioned in the literature (around 15 MPa at 293 K [START_REF] Ballet | Strain determination in quasi-lattice-matched LWIR HgCdTe/CdZnTe layers[END_REF][START_REF] Capper | Mercury Cadmium Telluride -Growth, Properties and Applications[END_REF]), but this can be explained by the cryogenic temperature of 100 K as well as by the possible diffusion of zinc before removal of the CZT substrate.

After assembly to the hosting substrate, the level of stress does not change significantly, whether considering XRD measurements or FEM simulations. This confirms the relevance of the choice of AlN, the CTE of which is sufficiently close to that of silicon not to modify the stress in HgCdTe.

Regarding the evolution with the format of the detector, the FEM simulations are consistent with the analytical models and do not show a significant impact on the stress in HgCdTe (whether with or without a hosting substrate). It is thus a parameter whose modelling is quite reliable (uncertainty of only a few MPa), whatever the initial value. It is therefore simply necessary to keep the well-known recommendations: favor the lowest lattice parameter variations from the epitaxial growth to minimize the density of dislocations and be careful not to increase the stress level in HgCdTe with an unsuitable hosting substrate.

Warpage

For the hybridized detector only, a 0.5 μm difference exists at the initial state between the measurements and the FEM simulations carried out for the two samples D and E, but it is not significant since a variability in the initial state of flatness due to imperfection of the manufacturing process is expected. The most important is that the evolution with temperature is well predicted, with a very good correlation between experiments and simulations. The analytical models predict an evolution of warpage approximately twice bigger, but as it was shown by the FEM simulations, this difference is due to the geometrical simplification and the elastic assumption (FEM gives equivalent results with the analytical models only in the case of a multilayer plate under elastic assumption). However, the sensitivity analysis shows that part of the relationships that can be deduced from the analytical models remains relevant for the hybridized detector: the contribution of the layers is proportional to their thickness, as well as to their stiffness and their CTE as long as the yield stress is not reached.

The sensitivity study (Table 3) makes it possible to identify two highly influential parameters: the initial stress value in CMOS and the cure shrinkage of underfill epoxy (these two parameters lead together to an uncertainty of about 80% on the warpage value at 100 K). The consistency between simulations and measured values suggests that there were no big errors in the input data, especially when making the assumption of negligible cure shrinkage. Regarding the evolution with format, whatever the temperature, and whether with or without hosting substrate, the FEM simulations (Figs. 9b and10b) show a proportional relationship between warpage and squared diagonal of the detector: this implies a constant radius of curvature, in accordance with the analytical models.

After assembly to the hosting substrate, the warpage is reduced whatever the format as shown in Fig. 10, with a minimum around 110 K, very close to the operating temperature as expected, and a maximum around 200 K. The latter corresponds both to the glass transition temperature of the glue and to the temperature below which the CTE of AlN and Si diverge.

The gap between the two largest samples extrapolated (86 mm diagonal) and measured (13.4 mm diagonal) shows the limit of the models; the uncertainty in the simulated values for large formats is indeed very high (about 80%), with no experimental measurement for validation. For the hybridized detector only, considering the 2000 × 2000 format, the value reached at 100 K (140 μm) is such that it could not meet the flatness specifications around 20 μm for this type of detectors (optical requirements). For the hybridized detector on the hosting substrate, the predicted value (15 μm) may be sufficient to meet specifications, but this should be verified by direct measurements.

Unlike the stress in HgCdTe that does not change with the format, flatness is therefore an important parameter to take into account when enlarging the size of the detectors. If there does not seem to be a technological barrier for the moment, since there are many possibilities to limit variations in warpage with temperature (increase in the thickness of the hosting substrate, deposit of thin compensation layers at the rear …), a warpage measurement will be systematically necessary to ensure the performances at operating temperature.

Conclusions

The aim of this work was to estimate the thermomechanical behavior of large format infrared detectors by FEM based on a multi-method analysis of a 320 × 256, 30 μm pixel pitch IR detector, focusing on stress in the HgCdTe layer and warpage. The hybridized detector was studied at room temperature and 100 K, both alone and after being glued to an AlN hosting substrate. Then, the behavior of larger formats was extrapolated by FEM and compared to predictions from analytical models. The main results are:

-from XRD results, the yield stress of HgCdTe at 100 K is at least 57 MPa.

-the stress values in HgCdTe are all biaxial, with either compression or tension at room temperature mainly due to the lattice mismatch during epitaxy from CdZnTe, and with an increase of about 40-45 MPa when cooling down to 100 K (consistent between FEM, XRD and analytical models). The initial value is variable from one sample to another mainly depending on lattice mismatch during epitaxy. Gluing to AlN hosting substrate has a very low impact on the results.

-warpage increases by around 2,5 μm when cooling from room temperature to 100 K for the hybridized circuit alone and is less than 1 μm after being glued to the hosting substrate. The measurements on two samples fit well with FEM results. -the analytical models and the FEM local analysis show that the prediction of stress in HgCdTe is reliable, with only a few uncertainties coming from the CTE values and the level of initial stress due to the lattice mismatch during HgCdTe epitaxy on CdZnTe. The warpage values depend on more parameters, mainly the properties of the CMOS layer and the cure shrinkage of epoxy. -enlarging the format has no significant impact on stress in HgCdTe, but warpage is proportional to the squared diagonal of the sample: it is therefore a key parameter for the thermomechanical design and reliability of large format detectors in relation to the optical planarity requirements.

Credit author statement

Fig. 1 .

 1 Fig. 1. a) Architecture of the IR detector, b) typical cross-section of a CMOS layer of a ROIC.

Fig. 2 .Fig. 3 .

 23 Fig. 2. a) Cryostat schematic view, b) example of a measured surface (in red) and the paraboloid approximation (in black), and c) warpage convention used in this work. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fig. 4 .

 4 Fig. 4. Coefficient of thermal expansion α (instantaneous) and Young's modulus of the underfilling epoxy, the glue used for assembly to hosting substrate, indium and interconnection layer (EHM).

Fig. 5 .

 5 Fig. 5. Stress-strain curves of the interconnection layer at different temperatures with a subplot for values at low strain.

Fig. 6 .

 6 Fig. 6. Meshed geometry of the hybridized detector with HgCdTe in red, epoxy underfill residue and interconnection layer in green, CMOS layer in light blue and silicon in dark blue. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fig. 8 .

 8 Fig. 8. FEM results at 100 K (z displacement Uz, σ xx and σ yy ) for the hybridized detector, before (a, b, c) and after (d, e, f) gluing to hosting substrate.

Fig. 9 .

 9 Fig. 9. Evolution of warpage with the format of the hybridized detector as a function of a) temperature for several formats and b) squared diagonal for several temperatures (FEM results).

Fig. 10 .

 10 Fig. 10. Evolution of warpage with the format of the detector after gluing to the hosting substrate as a function of a) temperature for several formats and b) squared diagonal for several temperatures (FEM results).

  Lucas DUPERREX: Conceptualization, Methodology, Software, Validation, Formal analysis, Investigation, Resources, Writing -Original Draft, Writing -Review & Editing, Visualization, Project administration. Raphaël PESCI: Conceptualization, Methodology, Validation, Investigation, Resources, Writing -Original Draft, Writing -Review & Editing, Supervision, Project administration, Funding acquisition. Pascal LE BOTERF: Conceptualization, Methodology, Validation, Resources, Writing -Review & Editing, Supervision, Project administration, Funding acquisition. Olivier MAILLIART: Conceptualization, Validation, Resources, Writing -Review & Editing, Supervision, Project administration, Funding acquisition.

Table 2

 2 Coefficient of thermal expansion α (instantaneous) and Young's modulus of silicon, HgCdTe, AlN and CMOS (EHM).

	Material	Instantaneous CTE α (10 -6 .K -1 )	Young's modulus (GPa)
		293	100	Reference	293	100	Reference
		K	K		K	K
	Silicon	2.5	-0.3	[22] + isotropy	130	132	[24,25]
	(100)			[23]
	HgCdTe	5.03	1.9	[16,26]	50	52	[27,28]
	(111)
	AlN	2.6	0.5	[22,29]	310	309	[30]
	CMOS	3.5	1.2	EHM	84	83	EHM

Table 3

 3 FEM local sensitivity analysis on the main parameters (one-factor-at-a-time method): "low" indicates variations under 5%.

	Parameter	Nominal or unit	Value or variation imposed	Impact on results
				σ HgCdTe 100 K (51 MPa)	σ HgCdTe RT (12 MPa)	Warpage 100 K (3 μm)	Warpage RT (0.6 μm)
	CMOS hypothesis	EHM	CMOS=Si	low	low	-21%	-24%
	ε initial CMOS (biaxial)	0	CMOS=SiO 2 -2.10 -4 +2.10 -4	low low low	low low low	-62% -13% +13%	-65% -66% +66%
	Epoxy shrinkage	ΔL/L = 0	ΔL/L = 0.22%	low	low	+11%	+19%
	CTE HgCdTe	10 -6 /K	ΔL/L = 1% -10%	low -17%	low -23%	+71% -6%	+85% -9%
			+10%	+17%	+23%	+6%	-9%
	σ initial HgCdTe (biaxial)	0	-15 MPa	-15 MPa	-15 MPa	-10%	-48%
			+15 MPa	+15 MPa	+15 MPa	+10%	+47%

Acknowledgements

This work was supported by Association Nationale de la Recherche Technologique (grant number 2017-1323). It was done in the context of DEFIR, the common laboratory of LYNRED and CEA-Leti. The authors want to thank all the people involved in the fabrication and the characterization of the samples, especially X. Derveaux for the warpage measurements.

Data availability statement

Not all the raw/processed data required to reproduce these findings can be shared but some is available from the corresponding authors on reasonable request.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.


A Ndiaye
email: ndiaye@cerfacs.fr

M Bauerheim
email: bauerheim@cerfacs.fr

S Moreau

F Nicoud

UNCERTAINTY QUANTIFICATION OF THERMOACOUSTIC INSTABILITIES IN A SWIRLED STABILIZED COMBUSTOR

Combustion instabilities can develop in modern gasturbines as large amplitude pressure oscillations coupled with heat release fluctuations. In extreme cases, they lead to irreversible damage which can destroy the combustor. Prediction and control of all acoustic modes of the configuration at the design stage are therefore required to avoid these instabilities. This is a challenging task because of the large number of parameters involved. This situation becomes even more complex when considering uncertainties of the underlying models and input parameters. The forward uncertainty quantification problem is addressed in the case of a single swirled burner combustor. First, a Helmholtz solver is used to analyze the thermoacoustic modes of the combustion chamber. The Flame Transfer Function measured experimentally is used as a flame model for the Helmholtz solver. Then, the frequency of oscillation and the growth rate of the first thermoacoustic mode are computed in 24 different operating points. Comparisons between experimental

and numerical results show good agreements except for modes which are marginally stable/unstable. The main reason is that the uncertainties can arbitrary change the nature of these modes (stable vs unstable) ; in other words, the usual mode classification stable/unstable must be replaced by a more continuous description such as the risk factor, i.e. the probability for a mode to be unstable given the uncertainties on the input parameters. To do so, a Monte Carlo analysis is performed using 4000 Helmholtz simulations of a single experimental operating point but with random perturbations on the FTF parameters. This allows the computation of the risk factor associated to this acoustic mode. Finally, the analysis of the Monte Carlo database suggests that a reduced two-step UQ strategy may be efficient to deal with thermoacoustics in such a system. First, two bilinear surrogate models are tuned from a moderate number of Helmholtz solutions (a few tens). Then, these algebraic models are used to perform a Monte Carlo analysis at reduced cost and approximate the risk factor of the mode. The accuracy and efficiency of this reduced UQ strategy are assessed by comparing the reference risk factor given by the full Monte Carlo database and the approximate risk factor obtained by the surrogate models. It shows a good agreement which proves that reduced efficient methods can be used to predict unstable modes.

NOMENCLATURE

INTRODUCTION

Combustion instabilities remain an active research domain for industries that involve combustion processes. They refer to the sustained acoustic fluctuations in combustors where unsteady combustion takes place. Usually this problem arises once the full combustor is tested, thus requiring expensive modifications on either the fuel injection system or the combustion chamber itself. Therefore, combustion instabilities have to be tackled at the design stage. However, this problem remains difficult because of the non-linearities associated with the turbulent flow or chemical reactions, the complex interactions between acoustics and unsteady heat-release, impedance boundary conditions as well as the numerous parameters involved. This issue becomes even more challenging when considering uncertainties on input parameters and underlying models. Especially, flame models obtained experimentally or numerically are known to be uncertain. Addressing the sensitivity of thermoacoustic results with respect to the input parameters is thus a necessary and important step towards reliable predictions of unstable modes in gas turbines. This is illustrated in Fig. 1 showing a typical result of a ther- moacoustic analysis, i.e. a set of modes, each with its own frequency and growth rate. When no uncertainty is present, each mode corresponds to a single point (black symbols) in the frequency plane. Here, modes 1, 4 and 5 are dangerous and should be also controlled since the growth rate ω i is positive. If uncertainties are present, each mode belongs to an admissible region of the frequency plane. Mode 2 (and maybe 6) are now dangerous and should be controlled. It suggests that taking into account uncertainties is required for reliable predictions of combustion instabilities. The objectives of this paper is therefore to (1) determine the uncertain region of the acoustic modes of a single swirled burner combustor using a brute force approach and (2) propose a methodology to deal with uncertainties at lower cost. This paper is organized as follows : the thermoacoustic framework used in this study is described in section 2 and section 3 presents the single swirled combustor configuration. The fourth section deals with the brute force Monte Carlo methodology applied to evaluate the uncertain region associated to the acoustic mode of interest. A statistical convergence study shows that 4000 Helmholtz simulations are required. It suggests that such a method can become CPU-demanding for more complex configurations using multiple uncertain parameters. In section 5, two surrogate models are tested to avoid expensive brute force methods. The purpose is to compute uncertain region of the acoustic mode with an affordable number of simulations. Finally, the prediction and the accuracy of the low-order models are validated against the initial full Monte Carlo analysis in the section5.

Thermoacoustic framework

Thermoacoustic instabilities can be studied using several approaches : theoretical models, low-order network methods, full scale Large Eddy Simulations, or experiments. Analytical approaches [?] allow the study of the underlying processes leading to combustion instabilities. However, numerous assumptions (null Mach Number, linear framework etc.) are required to simplify the problem and hence only academic or simplified configurations can be investigated. These assumptions can be relaxed using Large Eddy Simulation to predict turbulent flames interacting with the acoustic waves and reproduce self-excited combustion instabilities. However, even when LES simulations confirm that a combustor is unstable, they do not suggest how to control the instability. Moreover, LES techniques are CPU demanding : for example, performing parametric studies using hundreds of simulations is out of reach today. Faster tools such as acoustic solvers [?] are thus required to simplify the design process of the system. To do so, a linear wave equation for small pressure perturbations p 1 (x,t) can be derived from the reactive Navier-Stokes equations by neglecting turbulence and viscous effects and assuming a frozen baseline flow [?] :

∇. 1 ρ 0 ∇p 1 - 1 γ p 0 ∂ 2 p 1 ∂t 2 = - γ -1 γ p 0 ∂ q 1 ∂t (1)
where q 1 is the heat release fluctuation. Assuming harmonic pressure (p 1 = p(x)e iωt ) and heat release fluctuations (q 1 = q(x)e iωt ), Eq. ( 1) becomes :

∇. 1 ρ 0 ∇ p + ω 2 γ p 0 p = iω γ -1 γ p 0 q(x) (2)
where p and q are the complex amplitude of the pressure and heat release disturbance and ω = 2π f is the complex angular frequency of the thermoacoustic mode. The density ρ 0 and the specific heat ratio γ can depend on space x and are known quantities related to the baseline flow ; the thermodynamic pressure p 0 is constant under the zero Mach number assumption. In Eq. ( 2), the right hand side term is closed by a Flame Transfer Function using the n-τ model [?].

The expression of the global FTF reads :

F (ω) = n(ω)e iωτ(ω) ( 3
)
where n is the amplitude of the flame response and τ corresponds to a time delay. These parameters n and τ are usually obtained from experiment or LES simulations. They are known to control combustion instabilities although being uncertain. One objective of this paper is to assess the sensitivity of the computed complex angular frequency to these uncertain parameters. To do so, a Helmholtz solver called AVSP [?] is used to solve the wave equation (Eq. ( 2)) in a 3D single swirled burner geometry. The outputs of the solver are the pressure field p and the complex angular frequency ω of the acoustic mode : its real part ω r corresponds to the frequency of oscillation and ω i is the growth rate of the acoustic disturbances. AVSP is used to study how the uncertainties on n and τ propagate into uncertainties on the growth rate ω i and to determine the risk factor of the acoustic mode i.e. the probability for a mode to be unstable (ω i > 0) :

risk factor(%) = 100 ∞ 0 PDF(ω i )dω i (4)
Where PDF(ω i ) stands for the Probability density function of the growth rate of the acoustic disturbances. To fairly assess the risk factor, it is necessary to have a realistic statistical distribution of the input parameters n and τ, given here by the experiment presented in the following section.

The single swirled burner configuration

The experimental configuration targeted in this paper corresponds to a single injector combustor designed and studied by Palies et al. [?] (Fig. 2). It comprises a confined swirled flame, an upstream manifold, an injection unit equipped with a swirler and a cylindrical flame tube. The geometry of the combustor is displayed in Fig. 2. In order to ease the numerical study, the system was decomposed in three dimensional coupled tubes with variable lengths. The length l 2 which corresponds to the axisymmetric convergent tube is fixed. The length l 1 corresponds to the upstream manifold which is designed to have three different lengths and the combustion chamber l 3 may have four different lengths. This leads to twelve different set-ups, presented in mentally tested leading to flame A and B, the former with a smaller power than the latter. Finally, 24 operating points were studied for this configuration. The stability of the system has been experimentally evaluated by [?] and compared with numerical results obtained by AVSP assuming no acoustic losses. These losses can be obtained experimentally [?] and are given for the two flames : α A = 82(s -1 ) for flame A and α B = 125(s -1 ) for flame B with an uncertainty of ∆α = ±10(s -1 ). Numerically, stable regime is considered if the growth rate ω i is smaller than the damping rate α, and unstable if the growth rate is larger than the damping rate. When considering the error ∆α, this classification becomes : --Stable S :

ω i < α -∆α --Unstable U : ω i > α + ∆α --Marginal S/U : α -∆α < ω i < α + ∆α
The comparison between the experiment and the numerical results provided by [?] reveals a good agreement as illustrated in (Table . 2). However, three partial disagreements are observed and correspond to cases where the experiment gives a marginal stability (S -U) while the computation predicts an instability (U) or vice versa. In the remaining of the paper, one focuses on the operating point which corresponds to the configuration 07 for the Flame B where partial disagreement has been found between the experiment and numerical simulation. Moreover, it is hardly conceivable to perform a computation for all the configurations as one Helmholtz computation takes almost 10 minutes when using 24 processors.

Case

Uncertainty quantification study

Computing the risk factor (Eq. 4) associated to the first acoustic mode ( f 0 = ω 0 2π = 121Hz), i.e. the probability for this mode to be unstable constitutes one objective of this paper. Using uncertainty quantification is justified since repeating some tasks in this experiment may yield different results. They can be related to a statistical distribution of the FTF uncertain parameters n and τ. In the context of non-intrusive methods where the acoustic code acts as a "black box", the choice of sampling points is essential. The next section details the classical brute force Monte Carlo analysis to provide a reference database for the growth rate probability density function (PDF), a necessary quantity to compute the risk factor (Eq. ( 4)). Then, surrogate models are fitted using few Helmholtz simulations extracted from the Monte Carlo database by solving a least mean square problem. These surrogate models are then reused to perform a Monte Carlo analysis at reduced cost and estimate the risk factor. Finally, a discussion on the prediction of the risk factor is carried out by evaluating the accuracy of the low-order models against the initial full Monte Carlo method.

Classical Monte-Carlo analysis

First, a Monte Carlo sampling is used on the two dimensional parameters space n and τ (Eq. 3) to evaluate the probability of the first acoustic mode ( f 0 = ω 0 2π = 121Hz) to be unstable. The literature does not confer a clear accurate analysis on the uncertainty range of these parameters and this may have an impact on the stability of the system. We thus analysed some unpublished data from experimentalists at EM2C (Paris) and IMFT (Toulouse), two leading groups regarding the experimental measurement of flame response. Typically, data corresponding to the same configuration but gathered at different days were analysed to obtain an assessment of the variability of n and τ : from this variation, standard deviation of the FTF parameters σ n and σ τ are fixed to The uniform PDF reads :

f U X = 1 [x max -x min ] 1 ||x max -x min || for x min ≤ x ≤ x max (5)
Therefore, the mean µ U X and the variance v U X are :

µ U X = x min + x max 2 and v U X = 1 12 (R U µ U X ) 2 (6)
where R U represents the normalized range x max -x min µ U X of the uniform distribution : here R U = 10%. The β -distribution is characterized by its density function :

f β Y = B(α, β ) -1 y α-1 (1 -y) β -1 for 0 ≤ y ≤ 1 ( 7
)
where B(α, β )= Γ(α)Γ(β ) Γ(α+β ) denotes the beta function, Γ(.) is the gamma function, and α and β are two free parameters. Note that f

α = µ β Y µ β Y (1 -µ β Y ) v β Y -1 (8)
and

β = (1 -µ β Y ) µ β Y (1 -µ β Y ) v β Y -1 (9)
To close the problem, the reduced variable Y in [0, 1] is related to the desired random variable X in [x min , x max ] :

X = µ β X (1 + R β [2Y -1]) (10)
Taking the mean and variance of the previous equation leads to the following relations between characteristics of X and Y :

µ β Y = 1/2 and v β Y = ν β X 4R 2 β (µ β X ) 2 (11)
the same characteristics, i.e. µ β X = µ U X and ν β X = ν U X . Note however, that the range of the β -PDF appears in (µ β X )(Eq. ( 11)). If this range is chosen equal to the range of the previous uniform PDF (i.e. R β = R U = 10%) then the β -distribution degenerates to the previous uniform PDF. Consequently, the range R β is an additional free parameter. For this study, this range is fixed to R β = 30% leading to the characteristic values α = β = 2.87. In the following sections, the results from these two distributions are described and the statistical convergence of the risk factor is discussed.

The uniform distribution

The Monte Carlo sampling is displayed in Fig. 4 for the configuration 07 of the Flame B. Each point of Fig. 4 corres- ponds to a Helmholtz simulation in the complex domain. The horizontal solid lines denotes the acoustic losses α : 115 s -1 < α B < 135 s -1 . The stable or unstable regions are evaluated using the difference ω iα :

1. ω i -115 s -1 < 0 corresponds to a stable system (S).

2. ω i -135 s -1 > 0 corresponds to a unstable system (U). 3. 115 s -1 < ω i < 135 s -1 corresponds to a situation where the system is marginal (neither stable nor unstable) (S/U). The 4000 samples are then classified into three types : stable regime (S), unstable regime (U) and marginal regime (S/U). The histogram of the growth rate is shown in Fig. 5 and corresponds to an approximation of the growth rate Probability Density Function PDF(ω i ). Most of the thermoacoustic modes found by AVSP are in the stable regime. This leads to a risk factor close to 24 %.

The β -distribution

Following a similar methodology, 4000 runs are performed using the AVSP code when n and τ follow a β -distribution. A stability analysis is performed to classify each Helmholtz solution as stable (S), unstable (U) and marginal (S/U) (Fig. 6). As expected, the scatter plot is more diffused than in the case of the Uniform PDF because extrema values are less expected than intermediate ones when using the β distribution. The histogram of the growth rate is plotted in Fig. 7 where the medium value of the acoustic loss is reported ( ᾱB = 125 s -1 , solid vertical line). From this histogram, a risk factor close to 22 % is obtained. The risk factor obtained from the Uniform and Beta distribution are close (24 % vs 22%) showing that UQ results are weakly affected by the distributions chosen for the input parameters n and τ for this case. This suggests that assessing the risk factor of a mode without a clear knowledge of the uncertainties on the input data is relevant.

Convergence study

For both the uniform and β -distribution, a brute force Monte Carlo analysis is performed using 4000 runs. Nevertheless, the question of the convergence of the resulting risk factor is still open. Consequently, a convergence study is realized by doing the same exercise but with other dataset size : from 100 to 4000. The results of this study are displayed in Fig. 8. In both figures, the dashed line represents the risk factor of the system determined by the full Monte Carlo database with 4000 runs and the solid line corresponds to the variation of the risk factor. This study reveals that 4000 samples are by far sufficient to reach a reliable convergence of the risk factor. This method is relatively simple to implement and affordable when few parameters (in our case n and τ) are involved but becomes expensive when the dimensionality increases (this is typically the case for annular combustors fed by several burners) or when the CPU cost for each computation is large. This suggests that reduced methodologies have to be introduced to obtain risk factors in more complex configurations : this is the second objective of this paper discussed in the following.

Multiple linear regression

Because Eq. ( 2) is an eigenvalue problem which is nonlinear in the pulsation ω i , the surface response ω i = ω i (n, τ) is implicit and non-linear. To speed up the UQ analysis, it is worth investigating if this surface response can be approximated by explicit surrogate models. Two such models are tested in this paper :

1. One based on the n-τ parameters :

Y model 1 = β 0 + β 1 n + β 2 τ (12)
2. One based on the Flame Transfer Function evaluated at ω = ω 0 where ω 0 corresponds to the mode without flame coupling. FTFs incorporate here physical non-linearities into the model :

Y model 2 = β 0 + β 1 ℜ(ne jω 0 τ ) + β 2 ℑ(ne jω 0 τ ) (13)
These two algebraic models can be written in linear algebra notation as follows :

Y = Xβ + ε = Y model + ε (14)
where Xβ is the matrix-vector product, β = [β 0 , β 1 , β 2 ] T corresponds to the regression coefficients of the model. These coefficients represent the mean change in the response variable for one unit of change in the predictor variable. Y is a N × 1 dimensional vector containing the growth rate ω i determined from N Helmholtz computations, X is the matrix containing 1, n and τ for each samples and ε the N×1 vector of residuals :

Y =      Y 1 Y 2 . . . Y N      , X =      1 n 1 τ 1 1 n 2 τ 2 . . . . . . . . . 1 n N τ N      , β =   β 0 β 1 β 2   and ε =      ε 1 ε 2 . . . ε N     
Since models of Eq. ( 12) are linear, a least squares methodology is used to determine the coefficients β which minimize the error ε :

β = X t X -1 X t Y (15)
where β corresponds to the estimated parameters from the least squares, (X t X) -1 is called the "information matrix" and X t corresponds to the transpose of the X matrix. The predicted values Ỹ for the mean of Y are then determined as follows :

Ỹ = X β = X X t X -1 X t Y (16)
A Monte Carlo analysis is performed at reduced cost to provide an approximate risk factor for the mode of interest f 0 = ω 0 2π = 121Hz. A two-steps UQ strategy is investigated in the following :

1. Find the three regression coefficients associated to the models at reasonable cost using Eq (15) and only few Helmholtz computations .

2. Use the two surrogate models to perform a Monte Carlo analysis for assessing the risk factor.

The correlation coefficient is evaluated to measure the correlation between the two surrogate models (Eq (12) and Eq (13)) and the reference Monte Carlo database obtained in section 4 : The correlation is 95% for the linear model 1 (Eq (12)) and 98% correlation for the linear model 2 (Eq (13)). The corresponding growth rates are plotted against their exact counterpart in Fig. 9. The algebraic surrogate model using the FTF formulation (Eq (13)) is rather accurate in mimicking the actual surface response of the system. That is the reason why the latter is then used to plot the histogram of the growth rate based on the reference Monte Carlo database using 4000 runs for both the Uniform and the Beta distribution, see Fig. 10. The shape of both histograms is very similar to the histograms plotted previously in section 4 for the Uniform and β distributions. Also, the associated risk factor is approximated using the surrogate model 2 : (i) for the Uniform distribution the RF is close to 23% (ii)for the β -distribution the RF is close to 21% . The next step is to use Eq (13) with the aim to find a way to approximate the risk factor at low cost i.

Corr 2 = ∑ i (Y model -Ȳ ) 2 ∑ i (Y -Ȳ ) 2 (17)

CONCLUSIONS

In this paper, a UQ study applied to thermoacoustic instabilities in a single swirled combustor experiment has been carried out. An Helmholtz solver has been used to analyze the thermoacoustic modes of the combustion chamber. The FTF measured experimentally is used as a flame model for the Helmholtz solver. The frequency of oscillation as well as the growth rate of the first thermoacoustic mode are computed in 24 different operating points. The stability of the system has been evaluated by [?] and predictions show a good agreement with experimental behavior of the combustor. However, it also reveals that uncertainties on the Flame Transfer Function parameters may lead to modes which are marginally stable/unstable. Towards a more accurate mode classification instead of the usual stable/unstable one, a continuous description has been adopted based on the risk factor defined as the probability for a mode to be unstable given the uncertainties on the input parameters. The risk factor associated to the first acoustic mode of the combustor has been evaluated using a Monte Carlo analysis based on 4000 Helmholtz simula-tions of a single experimental operating point but with random perturbations on the FTF parameters. The analysis of the Monte Carlo database suggests that a two-step UQ strategy may be efficient to deal with thermoacoustics in such a system. First, two bilinear surrogate models were tuned from a moderate number of Helmholtz solutions (a few tens). Then, these algebraic models were used to perform a Monte Carlo analysis at reduced cost and approximate the risk factor of the mode. A comparison between the reference risk factor given by the full Monte Carlo database and the approximate risk factor obtained by the surrogate models shows a good agreement which proves that reduced efficient methods can be used to predict unstable modes.

q 1

 1 Heat release fluctuations p Harmonic complex pressure amplitude q Harmonic heat release fluctuations p 0 Mean pressure c 0 Mean field of sound speed n FTF amplitude n Mean FTF amplitude f U Probability density function of the continuous uniform distribution f β Probability density function of the β -distribution B The beta function Γ The gamma function R U Range of the uniform distribution R B Range of the β -distribution ℜ Real part ℑ Imaginary part Corr Correlation coefficient Subscripts X Random variable Y Random variable Ỹ The mean of Y t Transpose operator Greek letters ρ 0 Mean density ω Complex angular frequency ω 0 Complex angular frequency without flame coupling ω r Real part of ω : Pulsation of acoustic pressure wave ω i Imaginary part of ω : Growth rate of the acoustic pressure disturbances τ FTF time delay τ Mean FTF time delay γ Specific heat ratio α Shape parameter of the β -distribution β Shape parameter of the β -distribution α A Acoustic losses rate of the flame A α B Acoustic losses rate of the flame B ∆ α Uncertainty of the acoustic losses rate µ U Mean value of the uniform distribution ν U Variance of the uniform distribution µ β Mean value of the β -distribution ν β Variance of the β -distribution

FIGURE 1 .

 1 FIGURE 1. Location of the first six thermo-acoustic modes in a typical combustor.

FIGURE 2 .

 2 FIGURE 2. Numerical/experimental configuration. From ref. [?]

  σ n n = σ τ τ = 10%, where n = 1080J/m and τ = 4.73ms correspond to the nominal experimental values. From these very scarce data, we decided to keep 10 % uncertainty for both input. In absence of more information regarding the probability density functions, two different distributions have been considered : a Uniform distribution (Section 4.2) and a β -distribution (Section 4.3) with the same mean and variance. The ranges of the uniform distributions are directly deduced from the experimental values of the amplitude and time delay i.e 10% of the mean values (Fig. (3)).

FIGURE 3 .

 3 FIGURE 3. The uniform and the β -PDF of an arbitrary random variable X with similar mean (µ) and standard deviation (σ ), but with different range (R)

βY

  is only defined for a reduced random variable Y on [0, 1]. The parameters α and β which characterize the β -PDF are deduced from the desired mean µ β Y and variance ν β Y of this reduced variable Y :

FIGURE 4 .

 4 FIGURE 4. Monte Carlo results for the configuration 07 for the Flame B. Uncertain region for the first acoustic mode.

FIGURE 5 .

 5 FIGURE 5. Histogram of the growth rate of acoustic disturbance for 4000 samples using the Uniform distribution for the parameters n and τ.

FIGURE 6 .

 6 FIGURE 6. Monte Carlo results for the configuration 07 for the Flame B. Uncertain region for the first acoustic mode.

FIGURE 7 .

 7 FIGURE 7.Histogram followed by the growth rate of acoustic disturbance for 4000 samples using the Uniform distribution for the parameters n and τ.

FIGURE 8 .

 8 FIGURE 8. Convergence study of the risk factor ( in %) ; Top : Uniform distribution, Bottom : β -distribution

FIGURE 9 .

 9 FIGURE 9. Multiple linear regression computation ; Top : Using the linear model 1 Eq (12), Bottom : Using linear model 2 Eq (13).

  e. relying on much less then 4000 Helmholtz computations. To assess the number of Helmholtz computations required to tune the surrogate model, several tuning of the β coefficients are performed based on 3, 5, 10, 40, 70 and 100 Helmholtz computations (instead of 4000

FIGURE 11 .FIGURE 12 .

 1112 FIGURE 11. Risk factor from Eq (13) ; Left : 3 points calculations are used, Middle : 5 points calculations are used, Right :10 points calculations are used

FIGURE 10 .

 10 FIGURE 10. Histogram followed by the growth rate using algebraic surrogate model 2( Eq (13)) : Top : With the full Monte Carlo database using the Uniform Distribution, Bottom : With the full Monte Carlo database using the β -distribution

Table . 1

 . . Additionally, two different air flow rates were experi-

	Cases studied		l3=100 l3=150 l3=200 l3=400
	Expe./Simu.	l1=96.0	C01	C02	C03	C04
	Expe./Simu.	l1=160.0	C05	C06	C07	C08
	Expe./Simu.	l1=224.0	C09	C10	C11	C12

TABLE 1 .

 1 Twelve

different configurations explored : l 1 indicates the upstream manifold length and l 3 corresponds to the combustion chamber length. Dimensions are given in millimeters. From ref.

[?]

TABLE 2 .

 2 Linear

		Flame A				Flame B
		C01	C02 C03 C04	C01	C02 C03 C04
	Experiment	S	S	S	U	S	S	S-U	U
	Simulation	S	S	S	U	S	S	S-U	U
		C05	C06 C07 C08	C05	C06 C07 C08
	Experiment	S	S	S-U	U	S	S	S	UU
	Simulation	S	S	S-U	U	S	S	S-U	U
		C09	C10 C11 C12	C09	C10 C11 C12
	Experiment	S	S	S-U	U	S	S	S-U	U
	Simulation	S	S	U	U	S	S	S	U

stability analysis of flame A and flame B. Comparison between experimental and numerical results. (S) Stable, (S/U) Marginally stable/unstable, (U) Unstable. The geometrical configurations C01 to C12 are defined in Table.

1.

ACKNOWLEDGEMENTS

This study was performed within the UMRIDA project funded by the European commission (FP7-AAT-2013-RTD-1 -605036). The authors also thank the Center for Turbulence Research for its support during the 2014 Summer Program. Prof. G. Iaccarino (Stanford University) and P. Constantine (Umines Colorado) are acknowledged for their help regarding the UQ strategy.