L Manchon
email: louis.manchon@ias.u-psud.fr

K Belkacem

R Samadi

T Sonoi

J P C Marques

H.-G Ludwig

E Caffau

Influence of metallicity on the near-surface effect on oscillation frequencies

Keywords: asteroseismology -convection -stars, low-mass -stars, oscillations -stars, solar-type

Context. The CoRoT and Kepler missions have provided high-quality measurements of the frequency spectra of solar-like pulsators, enabling us to probe stellar interiors with a very high degree of accuracy by comparing the observed and modelled frequencies. However, the frequencies computed with 1D models suffer from systematic errors related to the poor modelling of the uppermost layers of stars. These biases are what is commonly named the near-surface effect. The dominant effect is thought to be related to the turbulent pressure that modifies the hydrostatic equilibrium and thus the frequencies. This has already been investigated using grids of 3D hydrodynamical simulations, which also were used to constrain the parameters of the empirical correction models. However, the effect of metallicity has not been considered so far. Aims. We aim to study the impact of metallicity on the surface effect, investigating its influence across the Hertzsprung-Russell diagram, and providing a method for accounting for it when using the empirical correction models. Methods. We computed a grid of patched 1D stellar models with the stellar evolution code CESTAM in which poorly modelled surface layers have been replaced by averaged stratification computed with the 3D hydrodynamical code CO 5 BOLD. It allowed us to investigate the dependence of both the surface effect and the empirical correction functions on the metallicity. Results. We found that metallicity has a strong impact on the surface effect: keeping T eff and log g constant, the frequency residuals can vary by up to a factor of two (for instance from [Fe/H] = +0.0 to [Fe/H] = +0.5). Therefore, the influence of metallicity cannot be neglected. We found that the correct way of accounting for it is to consider the surface Rosseland mean opacity. It allowed us to give a physically grounded justification as well as a scaling relation for the frequency differences at ν max as a function of T eff , log g and κ. Finally, we provide prescriptions for the fitting parameters of the most commonly used correction functions. Conclusions. We show that the impact of metallicity through the Rosseland mean opacity must be taken into account when studying and correcting the surface effect.

Introduction

Thespace-bornemissionsCoRoT (Baglin et al.2006;Michel et al. 2008;[START_REF] Auvergne | [END_REF] and Kepler (Borucki et al. 2010) have provided a rich harvest of high-quality seismic data for solar-like pulsators. This has allowed a leap forward in our understanding and modelling of low-mass stars (see the reviews by Chaplin & Miglio 2013;Hekker & Christensen-Dalsgaard 2017). However, for the last three decades (e.g. Dziembowski et al. 1988) it has been known that the comparison between modelled and observed acoustic-mode frequencies suffer from systematic discrepancies. This bias is called the surface effect and has been widely studied in the solar case (Rosenthal et al. 1995;Christensen-Dalsgaard & Thompson 1997;Rosenthal & Christensen-Dalsgaard 1999). They are attributed to our deficient modelling of the uppermost layers of stars with a convective envelope. Indeed, 1D stellar models hardly take into account the complexity of these layers that are subject Present address: Institut d'Astrophysique Spatiale, Université Paris Sud, Orsay, France to highly turbulent flows as well as a complex transition between a convective to a radiative-dominated energy flux transport (e.g. Kupka & Muthsam 2017).

More generally, these frequency residuals prevent a direct comparison between modelled and observed frequencies. Frequency combinations are commonly used to circumvent this problem (e.g. Roxburgh & Vorontsov 2003), but still, an accurate determination of frequencies is highly desirable to take advantage of the full potential of asteroseismology. To reach this goal, a handful of empirical prescriptions with adjustable free parameters have been proposed (Kjeldsen et al. 2008;[START_REF] Ball | [END_REF], 2017;Sonoi et al. 2015) and allow one to apply a posteriori corrections to the modelled frequencies. Such an approach is now widely used (e.g. Lebreton & Goupil 2014;Silva Aguirre et al. 2017) and has proven to be quite efficient in inferring a stellar model that fits the observed frequencies. However, it suffers from some fundamental drawbacks. The choice of the parameters is not physically motivated. Consequently, there is no guarantee that this optimal model is unique and accurate (i.e. that it properly reproduces the real physical structure of the observed star).

Another complementary approach then consists of investigating the physical nature of the surface effect. This motivated a number of studies to unveil and constrain the physical ingredients responsible for these biases. More precisely, surface effect has been shown to be the result of two distinct effects (e.g. Houdek et al. 2017): structural effects coming mainly from turbulent pressure in the hydrostatic equation which is usually absent in 1D stellar evolution codes, and modal effects gathering modifications of the eigenmodes, mostly due to nonadiabaticity (e.g. Balmforth 1992;Houdek et al. 2017) as well as the perturbation of turbulent pressure induced by the oscillations (Sonoi et al. 2017). Other related processes were also invoked, such as convective backwarming (Trampedach et al. 2017) or magnetic activity (Piau et al. 2014;Magic & Weiss 2016). Nonetheless, as demonstrated by the early work by Rosenthal & Christensen-Dalsgaard (1999) on the Sun using a 3D hydrodynamical simulation, the dominant physical ingredient is thought to be the turbulent pressure that modifies the hydrostatic equilibrium and subsequently introduces an elevation of the star surface. Then, the acoustic cavity is modified and therefore the frequencies are as well.

Based on a grid of 3D numerical simulations, this method was used by Sonoi et al. (2015), Ball et al. (2016), Trampedach et al. (2017) who investigated the surface effect variations across the Hertzsprung-Russell diagram. These works clearly demonstrated that surface effects sharply depend on effective temperature and surface gravity of star. In addition, Sonoi et al. (2015) presented a way to provide parameters for the empirical surface corrections by fitting them against a physically motivated scaling relation derived by Samadi et al. (2013). However, all these works considered solar metallicity models while the distribution of metallicity for observed solar-like pulsators is quite large (see e.g. Pinsonneault et al. 2014). Our goal is thus to study the influence of metallicity on the surface effects and propose a method to account for it.

The article is organized as follows: in Sect. 2 we describe the method of model patching, which is constructed by replacing the upper layers of a 1D model by horizontally averaged stratification of a 3D model atmosphere, and our set of models. Then in Sect. 3 we show that metallicity has a strong impact on the frequency residuals and therefore its influence cannot be ignored. We also study the variation of the frequency differences with effective temperature, surface gravity and opacity and give a theoretical justification. Finally, in Sect. 4 we provide constraints on the various parameters usually used in the empirical surface effect function across the T efflog glog κ space.

Model-patching method

In this section we explain the method we used to patch our models and describe our final set of models.

Grid of 3D models

We used a grid of 3D hydrodynamical models from the CIFIST grid of stellar atmosphere including the superadiabatic region to the shallowest layers of the photosphere, computed using the CO 5 BOLD code (see Ludwig et al. 2009;Freytag et al. 2012 for details). The chemical mixture is based on the solar abundances of Grevesse & Sauval (1998) apart from the CNO elements which follow [START_REF] Asplund | Cosmic Abundances as Records of Stellar Evolution and Nucleosynthesis[END_REF]. We considered 29 models with effective temperature (T eff ) ranging from 4500 K to 6800 K, a surface gravity (log g) ranging from 3.5 to 4.5, and a metallicity [Fe/H] = -1.0, -0.5, +0.0, +0.5. We note that 5000 5500 6000 6500

T eff [K] 3.50 1 summarizes the global parameters of the 3D models. The range of metallicities we considered corresponds to the metallicities of observed solar-like pulsators (Anders et al. 2017;Serenelli et al. 2017). Table 1 exhibits small groups of models (labelled with same first letter) with very similar T eff and log g. Those groups for instance in Fig. 1 at log g = 4.0, help us to investigate the influence of metallicity on the surface effect by keeping other global parameters fixed. However, we pointed out that, whereas within a group the dispersion in log g is rather small (of the order of 0.1%), the dispersion in T eff is much higher (of the order of 1%). Indeed, surface gravity is an input parameter of the hydrodynamical simulations while effective temperature is controlled by adjusting the entropy at the bottom of 3D models. It is therefore difficult to match an accurate effective temperature.

Computation of patched models

For each 3D model, both a patched (hereafter PM) and an unpatched model (hereafter UPM) have been constructed. A patched model is a model computed using a 1D stellar evolutionary code in which we replaced the surface layers with the stratification obtained by horizontally averaging a 3D model computed with a R-MHD code. The fully 1D model is called an unpatched model. The construction of PM and UPM has been widely described in Trampedach (1997), Samadi et al. (2007Samadi et al. ( , 2008)), Sonoi et al. (2015), Jørgensen et al. (2017). The 1D counterparts of 3D hydrodynamical models have been obtained using the 1D stellar evolutionary code CESTAM (Morel 1997;Marques et al. 2013) by tuning the age (or the central temperature for advanced stages), the total stellar mass M, and the mixing length parameter α MLT in order to match the effective temperature, the surface gravity and the temperature at the bottom of the 3D model, located just below the superadiabatic region. We chose to remove the first four bottom layers and the last top layer of the 3D hydrodynamical model to be sure to remove any numerically induced errors and that the patching point is deeply inside the Notes. The final three letters of the model labels correspond to [Fe/H]: (m00, m05, m10, p05) refer to [Fe/H] = (0.0, -0.5, -1.0, +0.5), respectively. T b is the mean temperature at the bottom of the 3D model and ν max is the frequency with the largest amplitude in the oscillation power spectrum (ν max = 3050.0(M/M )(R 2 /R 2 )(5777/T eff ) 1/2 , Kjeldsen & Bedding 1995) and M is the mass of the PM which differs by a fraction 10 -7 from the one of UPM. The initial helium and metal abundances are close to the ones at the surface. We recall that the metal abundance is different from the iron abundance [Fe/H] imposed in our models. The evolutionary stages PMS, MS, and SG stand for pre-main-sequence, main-sequence and sub-giant.

adiabatic region, which has been shown to be a condition for a obtaining reliable PM (Jørgensen et al. 2017). The 1D models use the equation of state, and opacities given by OPAL2005 (Rogers & Nayfonov 2002;Iglesias & Rogers 1996) and implement standard mixing-length theory (Böhm-Vitense 1958) with no overshoot. We ignore diffusion processes, rotation and turbulent pressure. The atmosphere is computed using the Eddington approximation. The helium abundance in 1D models is set to the one used in 3D models.

Finally, we note that for some 3D models one can find a degenerate solution for the corresponding 1D model: we could patch either a PMS or a sub-giant model. We opted for evolved models since they corresponds to stars in which solar-like oscillations are observed so far. However, when the evolved models are too old (older than the age of the Universe) we kept the PMS model, except if lying on the Hayashi track.

Table 1 also summarizes the stellar parameters of both UPM and PM together with relative radius differences R PM /R UPM -1. Our set of models covers a wide portion of the Hertzsprung-Russel Diagram for intermediate mass stars. We note that our patched models with metallicity [Fe/H] = +0.5 only have log g = 4.0. Indeed, 3D models from the CIFIST grid with [Fe/H] = +0.5 were only available for log g ≥ 4.0. In addition, 3D models with log g 4.5 are located below the main sequence diagonal, and therefore it is impossible to find a 1D model matching their characteristics (with the physical ingredients we used). Thus, a large portion of our initial [Fe/H] = +0.5 3D hydrodynamical models were not suitable for our purposes.

Computation of oscillation frequencies

In this work, we consider only structural effects and an adiabatic treatment of the oscillations. The frequencies are computed using the ADIPLS code (Christensen-Dalsgaard 2011) for both UPM and PM by assuming the gas Γ 1 approximation, which assumes that the relative Lagrangian perturbations of gas pressure and turbulent pressure are equal (Rosenthal & Christensen-Dalsgaard 1999;Sonoi et al. 2017). Besides this distinction in the treatment of Γ 1 entering the calculation of the model frequencies, we emphasize that the frequency differences studied in this work are only emerging from structural effects. Therefore, it should be emphasized that the frequency differences studied in this work concern only purely structural effects. We have checked that we recovered the previous results of Sonoi et al. (2015) for the solar metallicity. For the sake of simplicity, we mainly focussed on the surface effect affecting radial modes: non-radial modes exhibit a mixed behaviour that would make our analysis more complex (however, see Sect. 4.2.2 for a discussion).

Influence of metallicity

Until now, surface effects have always been studied assuming a solar metallicity. Corrections depend only on T eff and log g such as the power law proposed by Kjeldsen et al. (2008), cubic and combined inverse-cubic laws [START_REF] Ball | [END_REF], or a modified Lorentzian (Sonoi et al. 2015). This section is intended to motivate the investigation of the dependence of the surface effect on metallicity.

Qualitative influence of metallicity on frequency differences

We begin this section by quickly describing the effects of a change of metallicity on the frequency residuals. Frequency differences are induced by the surface elevation between PM and UPM due to turbulent pressure, which extends the size of the resonant cavities and therefore decreases the mode frequencies for PM, leading to negative frequency differences δν = ν PM -ν UPM .

Up to now, only the influence of effective temperature and surface gravity on surface effects have been studied. However, the abundance of heavy elements has a strong impact on opacity and hence on the convective flux imposed by a change in the radiative flux. In turn, a change in the convective flux leads to a change of convective velocity and therefore a change of turbulent pressure and finally it changes the location of the surface. We mention here that metallicity also has an effect on gas pressure, through the mean molecular weight µ, which varies in the opposite direction of the turbulent pressure and therefore counteracts its effect. Finally, while mechanisms by which a change of metallicity can act on the surface effect are known, those mechanisms are too intricate to identify the resulting effect on the variations of surface term without a deeper analysis as will be demonstrated in the following (see Sect. 4).

Figure 2 shows the (purely structural) frequency differences for three groups of models that have approximately the same effective temperature and surface gravity. The discrepancies in ν n between two models appear at relatively low frequencies and generally increase towards high frequencies. As for finding a general trend of the evolution of the surface effect against the metallicity, it seems from Fig. 2 no such trend exists: in the top panel, frequency differences, at ν max for instance, slightly decrease from [Fe/H] = -0.5 to 0.0 and then are much higher for the [Fe/H] = +0.5 model. In the middle panel, the frequency residual at ν max significantly increases from [Fe/H] = -0.5 to 0.0. Finally, in the bottom panel, very little variations at ν max can be noticed from one composition to an other. However, the variation of the frequency differences seems to follow closely the variations of the elevation of the stellar surface between UPM and PM:

∆r ≡ R PM -R UPM .
(1)

The slight disagreement between [Fe/H] = 0.0 and -0.5 in the top panel may be explained by the large dispersion in effective temperature.

Effect of the elevation on the frequency differences

To gain some insight into the influence of metallicity on surface effect, we tried to scale the normalized frequency differences at ν max for our set of models. This is a necessary step to allow an estimate of the surface effect correction parameters (see Sect. 4). Thus, let us start with the perturbative approach as adopted by (Christensen-Dalsgaard & Thompson 1997; see also Goldreich et al. 1991;Balmforth et al. 1996). The authors show that the frequency difference can be well approximated by

δν ν = R 0 Kn c 2 ,v δ m c 2 c 2 + Kn v,c 2 δ m v v dr (2) R 0 Kn v,c 2 δ m v v dr , (3)
where c is the adiabatic sound speed, the variable v is defined by v = Γ 1 /c, Kn c 2 ,v and Kn v,c 2 are the kernels that can be determined from eigenfunctions, δ m c 2 and δ m v are the Lagrangian differences of c 2 and v, respectively, at fixed mass.

Rosenthal & Christensen-Dalsgaard (1999) further approximated the frequency differences for radial modes, based on the expression of Kn v,c 2 and using a first-order asymptotic expansion for the eigen-function, by

δν ν ∆ν∆r c ph , ( 4
)
where ∆ν is the asymptotic large frequency separation, ∆r is the previously defined elevation, and c ph the photospheric sound speed (see Appendix A for a demonstration of this relation). This relation has been previously tested by Sonoi et al. (2015) at solar metallicity using surface effect derived from a grid of 3D numerical simulations. It turns out that Eq. ( 4) reproduces the overall scale of the surface effect (such as in Fig. 3 were the surface effect is considered at ν max ) for a set of models. It is thus necessary to determine whether this relation holds for models with a non-solar metallicity. To this end, we have compared frequency residuals at ν = ν max given by Eq. ( 4) as shown in Fig. 3 (top panel). There is still a good agreement between the frequency differences and the approximated expression given by Eq. ( 4). Moreover, it appears that the frequency differences are dominated by the surface elevation ∆r. To understand the link to metallicity, it is thus necessary to go a step further and to investigate the relation between surface elevation and metallicity.

Scaling law for the frequency differences

In this section, we aim to determine a relation between frequency differences at ν max and global parameters of the models. First, as shown in the previous section, there is no clear trend between the surface term and metallicity. Indeed, at constant metallicity and considering our rather large range of effective temperatures and surface gravities, the dominant opacity mechanisms are not the same from a model to an other for instance, the opacity at the surface is dominated by the negative hydrogen ions for T eff 5000 K. Therefore, the relation between δν/ν and Z is non-trivial. To overcome this problem, we directly consider the Rosseland mean opacity at the photosphere instead of the metallicity as a global parameter in addition to the effective temperature and to the surface gravity (in the following, the photosphere is defined as the radius at which T = T eff ).

Let us begin by considering the elevation in Eq. (1) which must be expressed as a function of these global parameters. Using the hydrostatic equilibrium equation, it reads

∆r = R PM 0 H PM p dp tot p tot - R UPM 0 H UPM p dp g p g , (5)
where H PM p and H UPM p are the pressure scale heights at the photosphere associated with the patched and unpatched models, p tot is the total pressure such as p tot = p turb + p g with p turb and p g the turbulent and gas pressure, respectively. Further assuming that

Finally, since the pressure scale-height scales as T eff /g, the elevation scales as ∆r ∝ (T eff p turb )/(gp g ).

To go further, we need to find an expression for p turb /p g . Near the photosphere, the turbulent pressure can be written as

p turb = ρv 2 conv , (7)
where v conv is the vertical component of the convective velocity. We now need an expression for this velocity and for the density. Assuming a standard Eddington grey atmosphere, the optical depth is approximated by τ = H p ρκ, and in the Eddington approximation, we have τ = 2 /3 at the bottom of the photosphere. Then, and accordingly:

ρ ∝ g T eff κ • (8)
As for finding an expression for v conv , we note that F tot = F rad +F conv , with F rad and F conv the radiative and convective component of the total energy flux respectively. The convective flux is proportional to the kinetic energy flux (as shown for instance within the MLT framework). Then,

ρv 3 conv ∝ T 4 eff 1 + F rad F conv -1 • (9)
The ratio F rad /F conv is assumed to remain nearly constant from one model to an other. Therefore, v conv finally reads,

v 2 conv ∝ T 8 /3 eff ρ 2 /3 • (10)
Inserting the expressions of Eqs. ( 8) and (10) into Eq. ( 7) leads to:

p turb ∝ T eff T eff 7 /3 g g 1 /3 κ κ -1 /3 , (11)
where κ = 0.415 cm 2 g -1 . From the perfect gas law for, p g ∝ ρT eff and using Eq. ( 8), we can rewrite ∆r as

∆r ∝ T eff T eff 10 /3 g g -5 /3 κ κ 2 /3 • (12)
Replacing ∆r into Eq. ( 4) one finally obtains the following estimate:

δν ν ∝ ∆ν ∆ν T eff T eff 17 /6 g g -5 /3 κ κ 2 /3 ≡ z 1 . (13)
This expression provides us with a simple relation between the frequency differences and the global parameters. The dependence on the metallicity is embedded into the Rosseland mean opacity. We note that it is possible to go further and to explicitly introduce the metallicity. For instance, in the vicinity of the solar effective temperature and gravity, the opacity is dominated by the H -so that κ ∝ ρ -1 /2 T 9 eff Z. However, given the wide range of effective temperatures and surface gravities of our grid of models, it is more relevant to keep the Rosseland mean opacity at T = T eff (surface opacity) as a global parameter. Indeed, the Rosseland mean opacity is a quantity available in any 1D stellar evolutionary code.

Then, using Eq. ( 13) as a guideline, we performed a fit where the powers of the temperature (p), gravity (q), and opacity (s) have been adjusted at ν = ν max for each model. Figure 3, bottom panel, displays the result. This figure shows a very good agreement between exponents derived in Eq. ( 13) and the one actually obtained using our simulations. Consequently this scaling can be used to provide a physically-grounded values for the parameters of the empirical correction function of the surface effect. Finally, we note that using the opacity instead of the metallicity allows us to take a detailed mixture into account.

In addition to our crude approximations, a possible source of discrepancies between values predicted by Eq. ( 13) and the one calculated can be that we did not fix the helium abundance from one model to the other when varying the metallicity. The changing helium abundances have an impact both on the evolution of the model and on its opacity at the surface. However, the helium abundances [He/H] range between -5.8 × 10 -3 and +1.2 × 10 -2 and should be a negligible source of uncertainty. A final source of error comes from the method we used to average the 3D stratifications. Indeed, since the Rosseland opacity is involved Eq. ( 13), it would be more precise to patch the models using a stratification averaged against the Rosseland optical depth instead of the actual geometrically averaged stratification, but this is beyond the scope of this paper and will be investigated in a forthcoming work.

Surface-effect corrections

A handful of empirical functions have been suggested to perform a posteriori corrections on the modelled frequencies. After having given a theoretical background that explains variations of δν/ν, we considered the most commonly used correction models to study the evolution of the related free parameters as a function of effective temperature, surface gravity, and surface opacity. This is intended to provide constraints on those parameters and thus to provide physically-grounded values for use on seismic observations.

Empirical functions for correcting modelled frequencies

4.1.1. Kjeldsen et al. (2008) power law Kjeldsen et al. (2008) proposed a power law which was found to match the frequency differences obtained between the observed and modelled solar frequencies:

δν ν max = a ν PM (n) ν max b , (14)
where a and b are the parameter to be adjusted. They found a = -4.73 and b = 4.9 for their model of the Sun by matching a subset of nine radial modes centred on ν max . Kjeldsen et al. (2008) provided a method to correct the frequency for a star similar to the Sun without having to calibrate b. Let us assume we want to model a star with near solar global parameters and we want to constrain our model using the individual frequencies. The radial mode frequencies spectrum of our best model which include a surface term are denoted ν i,best and the frequencies of solar radial modes for the same order are denoted ν i,ref . Then, Kjeldsen et al. (2008) proposed that the frequencies can be linked, to a good approximation, by ν i,best rν i,ref , using the proportionality factor r between mean densities of both models: ρbest = r 2 ρref . Using this relation and the large separations of both models, they provided a way to obtain a and b. Further assuming b constant (the value of which depends of the physical ingredients used in the model), they derived a value for a for a set of theoretical models close to the Sun.

This power law has been widely used since and many authors (e.g. Metcalfe et al. 2009;Bedding et al. 2010) have used a constant value for b (not necessarily 4.90 though) derived from solar frequency measurements. Keeping b constant is often necessary in the case for which observations do not provide enough constraints to adjust it. However, using the solar value leads to a bad correction if the modelled star is too different from the Sun (e.g. Kallinger et al. 2010). Furthermore, b depends on the input physics. Otherwise, b can be considered as a variable parameter in the modelling and therefore significantly improve the correction. Different models of the star HD 52265 have been compared by Lebreton & Goupil (2014) using various input physics and found approximatively the same predicted age models when either frequency ratios (Roxburgh & Vorontsov 2003) or individual corrected frequencies were used as constraints. The age dispersion was slightly higher with models constrained by individual corrected frequencies (∼±9.5%) and using uncorrected individual frequencies lead to ages 40% larger (Lebreton et al. 2014).

In the following, we have studied two versions of this parametric function. The first, adjusted on the whole radial mode frequency spectrum for frequency less than the acoustic cut-off frequency, will be referred to as K08. The second, adjusted on a reduced frequency interval 0 < ν/ν max < 1.05 is refered to as K08r (see Fig. 6 and Appendixes B and C).

Ball & Gizon (2014) cubic and combined inverse-cubic laws

Ball & Gizon ( 2014) suggested a new function to correct frequency differences. It is partially based on the early work by [START_REF] Gough | Progress of Seismology of the Sun and Stars[END_REF] for the cubic part). They accounted for two leading effects introducing systematic errors in the theoretical computation of the frequency spectrum: the modification of the sound speed caused by a magnetic field concentrated into a filament by convective motions, causing a frequency shift scaling as ν 3 /E [START_REF] Libbrecht | Observations of Solar Cycle Variations in Solar p-mode Frequencies and Splittings[END_REF]), E being the normalized mode inertia; and the modification of the pressure scale height caused by a poor description of convection, inducing a frequency shift scaling as ν -1 /E. This correction funtional has the advantages of being independent of a solar calibration and including a dependence on the normalized mode inertia which allows us to correct non-radial modes, without the need of re-scaling their frequency differences. Because of this, they suggested a cubic correction taking only into account the dominant effect and a combined inverse-cubic correction including the perturbation.

The cubic correction (in the following BG1) is defined by

δν ν max = a 3,BG1 E ν ν max 3 , ( 15
)
and the combined inverse-cubic correction (in the following BG2) is

δν ν max = 1 E        a -1,BG2 ν ν max -1 + a 3,BG2 ν ν max 3        , ( 16
)
where E is the normalized mode mass:

E = 4π R 0 |ξ r (r)| 2 + ( + 1)|ξ h (r)| 2 ρr 2 dr M |ξ r (R)| 2 + ( + 1)|ξ h (R)| 2 , ( 17
)
where R, M, and ρ are respectively the photospheric radius, mass and density of the star, and ξ r and ξ h are the radial and the horizontal component of the displacement of an eigenmode of degree . a 3,BG1 , a -1,BG2 , and a 3,BG2 are the parameters to be adjusted. They used the acoustic cut-off frequency ν c instead of ν max in order to normalize their fitting parameters: it only results in a modification of a -1 and a 3 and does not change the law itself.

Sonoi et al. (2015) modified Lorentzian

The final function to be introduced was a modified Lorentzian (Sonoi et al. 2015) that was found to better correct the surface effect derived from the 3D simulations at high frequency. It reads

δν ν max = α           1 - 1 1 + ν PM ν max β           , ( 18
)
where α and β parameters are to be determined. When ν PM /ν max 1 we get back to Kjeldsen et al. (2008) law. When ν PM = ν max , δν/ν max = α/2. Therefore, a and α are directly linked to δν/ν given by Eq. ( 13), which gives physical justification for its variations. In the following, we will refer to this correction law as S15. A recent comparison of the above correction laws has been performed by Ball & Gizon (2017) on six sub-and red giants from the Kepler Input Catalog. We note that since these are evolved stars, they display mixed-modes which have their frequency residuals off the general trend of radial p-modes frequency residuals. This should have consequences on the quality of the correction. They computed stellar models matching their six stars constrained by the effective temperature, the metallicity and the individual frequencies. They tested five correction relations: BG1, BG2, S15, and K08 with power b = 5.0 calibrated with a solar model computed with their own input physics and K08 with power b left free. Ball & Gizon (2017) found no correction to be clearly superior than the others for all stars. However BG2 and then BG1 performed slightly better than the others, followed by the free power law, S15 and finally K08. S15 was shown to poorly correct high frequencies and the K08 with b = 5.0 correction gave worse results than no correction for their of their stars. We present very similar conclusions in the following.

Variation of the coefficients in the ∆ν

-T eff -log g -κ space 4.2.

Prescriptions for radial modes

In order to fit the parameters of the correction functions and to determine the fitting parameters, we used a least square minimization algorithm implementing a Levenberg-Marquardt A107, page 8 of 14 method which minimizes the squared deviation defined as

D = N i = 1 ν PM,i -ν UPM,i -δν i ν max 2 , (19)
where i corresponds to the eigenmode index, N to the total number of radial modes and δν the correction computed from the considered correction relation. Tables B.1 andC.1 summarize the coefficients and their squared deviations from our computations. We also define √ D/N as the root-mean frequency differences after correction.

Coefficients a and b involved in Eq. ( 14) and a 3,BG2 from Eq. ( 16) are presented in T efflog g plane in Fig. 4. Furthermore, all coefficients are represented as a function of log z 1 in Fig. 5. a, α and a 3,BG2 (and therefore a 3,BG1 ) show similar trends. Indeed, they are related to the amplitude of the surface effect at ν max and Eq. ( 13) allows one to understand their variations. However, this theoretical justification for the variations of a, α, a 3,BG1 and a 3,BG2 does not provide a way to favour one correction law over an other. Coefficient a -1,BG2 exhibits the same behaviour as above. However, we cannot offer the same explanation for its trend because the inverse term in BG2 is a second correction to the cubic term and is not related to the amplitude of the surface effect at ν max .

The trends followed by b and β are related to the slope of the frequency differences. As shown in Sonoi et al. (2015) and in Fig. 4, the coefficients b (whatever the metallicity) increase significantly towards cooler stars, which again contradicts the assumption of a constant b. Regarding the relevance of giving a prescription for log b and log β thanks to the linear relationship with log z 1 , we can see in Fig. 5 that log b and log β are affected by a high dispersion compared to the grey line. This could mean either that we omitted a physical dependency in Eq. ( 13) that only affects the agreement with log b and log β, or a prescription based on same other physical basis should be investigated.

Table 2 shows the prescriptions for the variations in the T efflog gκ space of all coefficients c 0 studied in this article in the form log c 0 = c 1 log ∆ν/∆ν + c 2 log T eff /T eff + c 3 log g/g + c 4 log κ/κ + c 5 .

We note that the opacity has a strong impact on each of the coefficients and must, therefore be taken into account when correcting the surface effect.

The top panel of Fig. 6, top panel also shows the value of the root-mean frequency differences after correction for each model and each correction law. From this, we see that BG1 is the worst performer followed by K08. Those laws provide a correction that leaves frequency residuals comprised between 1 and 10 µHz which are still higher than the frequency resolution provided by CoRoT and Kepler. The better performance of K08 over BG1 can be explained by the fact that K08 have two degrees of freedom whereas BG1 has only one. For radial modes, the inclusion of the normalized mode mass E n in BG1 does not compensate the loss of a degree of freedom.

Then, the remaining laws K08r, S15 and BG2 provide correction almost as good as the resolution of CoRoT and Kepler. K08r and BG2 are slightly better than S15, yet K08r is applied only on the frequency range 0 < ν/ν max < 1.05.

Mixed-modes case

We also performed the same test as in Ball & Gizon (2017) on evolved models that present mixed-modes in their frequency spectrum. In Sect. 4.1.2 we see that, thanks to the dependence in the normalized mode inertia, BG1 and BG2 can be applied to non-radial modes without any change of the law. However, in order to be able to compare all empirical corrections on non-radial modes, one had to rescale the frequency differences on which K08, K08r and S15 by mean of the inertia ratio Q n for a mode of frequency ν n defined as the ratio of the inertia of this mode by the inertia of a radial mode interpolated at the frequency ν n : Rosenthal & Christensen-Dalsgaard 1999). Furthermore, we added one last empirical relation by modifying the expression given for S15 in Eq. ( 18) similarly to BG1 and BG2 in which we replaced α by α/E where E is defined in Eq. ( 17) (the new function is denoted S15E). This allows S15E to be applied directly on non-radial modes frequency differences. The empirical relations K08, K08r and S15 were then adjusted on Q n δν n , with 0 ≤ ≤ 2 and S15E, BG1 and BG2 were adjusted directly on δν n , with 0 ≤ ≤ 2. For 9 of the 16 evolved models considered, the least-square algorithm converge to a solution of K08, K08r or S15 very remote from the general trend of frequency differences, whereas for the second group of relations (S15E, BG1 and BG2), the residual root-mean frequency differences after correction are greatly improved to a value between 0.1 and 1 µHz.

Q n = E n /E n0 (ν n ) (e.g.
As a third test, we performed the same fits excluding the quadrupolar modes (i.e. we fit modes with 0 ≤ ≤ 1). This time, corrections laws accuracies are similar to the one presented in Sect. 4.2.1. As for the newly introduced S15E, it performs slightly worse than S15 but still better than K08. This third test suggests that the failure of K08, K08r and S15 in fitting Q n δν n , with 0 ≤ ≤ 2 is due to quadripolar modes. There are two reasons for it.

First, the p and g cavity are less coupled for = 2 than for = 1 mixed-modes which induces more important changes on the behaviour of a mode when the surface layers are changed between UPM and PM. Indeed, modifying the surface layers changes the frequency of pure p modes that couple with different g modes for PM and UPM (Ball & Gizon 2017). As a consequence, when computing Q n,2 δν n,2 , it so happens that we deal with mixed-modes from PM and UPM that have different properties. Second, due to the presence of mixed-modes, Q n,2 sometimes becomes higher than ten, while it is normally of the order of unity (see example of Cm10 in Fig. 7). It over-scales the corresponding quadrupolar mixed-modes and gives much weight

A m 0 0 A m 1 0 B m 0 0 B m 1 0 C m 0 0 C m 1 0 C p 0 5 D m 0 0 E m 0 5 F m 0 0 F m 0 5 F m 1 0 F p 0 5 G m 0 0 G m 0 5 G m 1 0 H m 0 0 I m 0 0 I m 0 5 J m 0 0 J m 0 5 J p 0 5 K m 0 0 K m 0 5 K m 1 0 L p 0 5 M m 0 0 M m 0 5 O p 0 5 Models 10 -2 10 -1 10 0 10 1 ν max D N [µHz] UC K08 K08r S15 BG1 BG2
Fig. 6. Root-mean frequency differences ν max √ D/N of radial modes after correction for each empirical law and for each model. It shall be noticed that the deviation for K08r is computed only on the range 0 < ν/ν max < 1.05. Black dots corresponds to uncorrected frequencies (UC), that is root mean frequency differences between PM and UPM frequencies.

to those modes, which in turns has strong impact on the quality of the fit. The peaks in the value of Q n,2 arise for mixed-modes having most of their amplitude in the g mode cavity, contrary to dipolar mixed-modes which have their Q n,1 staying close to unity. On the other hand, fitting directly δν n with S15E, BG1 and BG2 does not amplify the frequency differences affecting mixed-modes, providing a much better correction.

For all these reasons, we recommend the use of BG2 or S15E when correcting sets of radial and non-radial modes and we recommend either S15 or BG2 when correcting only radial modes (K08r can be used for frequencies ν max ). The remaining advantages of S15 over BG2 is that the coefficients α and β are better described as a function of Eq. ( 13) than coefficients of BG2 which present a bigger dispersion on Fig. 5. This being said, many of the correction laws considered in this paper gives a rootmean frequency difference of the order of 0.1 µHz (at least for few models), similar to the frequency uncertainties of CoRoT and Kepler. Furthermore, the search for the best a posteriori correction law should not set aside the need of a theoretical understanding of the surface effect.

Conclusion

We have computed a grid of 29 couples of one dimensional models using the method of patched models consisting in replacing poorly modelled surface layers of a 1D model by the stratification, averaged over the geometrical depth and time, computed from 3D hydrodynamical models. The grid includes models with effective temperature ranging from T eff = 5000 K to 6800 K, surface gravity ranging from log g = 3.5 to 4.5 and iron abundance ranging from [Fe/H] = -1.0 to +0.5.

Our aim was to estimate and understand the impact of varying metallicities on the surface effect. Our main result is that, in the considered range of metallicities (i.e. [Fe/H] = -1.0 to +0.5) the amplitude of the surface effect computed at ν max , and for models with same effective temperature and same surface gravity, can be up to a factor of three between the model with the lowest amplitude and the model with the highest one. However, it appears that studying the amplitude as a function of the metallicity does not lead to a clear trend, whereas the Rosseland mean opacity κ turned out to be the adapted quantity for understanding the variation of the surface effect. Based on relatively simple physical arguments, consolidated using the grid of 3D models, we found a scaling relation between the amplitude of the surface effect and the global parameters T eff , log g and the opacity κ computed at the photosphere.

We also tested the accuracy of existing surface effect empirical corrections of radial modes frequency differences on each model of our grid in order to obtain a prescription for the coefficients. Then, we tested those laws on radial and non-radial modes for evolved models exhibiting mixed-modes, in order to test how the empirical corrections perform when mixed-modes are involved. Overall, the combined correction law proposed by [START_REF] Ball | [END_REF] is found to give the best performer, closely followed by the law proposed by Sonoi et al. (2015). These two laws leave frequency differences that are less than 1 µHz on average, even reaching 0.1 µHz for the coolest stars of our set of model, which is of the order the frequency resolution provided by CoRoT and Kepler. We note that, on a low frequency range (0 < ν/ν max < 1.05), the Kjeldsen et al. (2008) power law (calibrated on this reduced range) gives equivalent results. Then the Kjeldsen et al. (2008) power law calibrated on the whole range of frequency and the purely cubic correction proposed by [START_REF] Ball | [END_REF] are the worst performer with remaining mean frequency differences of the order of few µHz. When applying those corrections on frequency spectra including mixed-modes, only the empirical corrections BG1 and BG2 proposed by [START_REF] Ball | [END_REF] and the modified S15E where we added a factor of 1/E improve the mean frequency dispersion. Only S15E and BG2 leave a satisfying root-mean dispersion of the order of the CoRoT and Kepler frequency resolution.

Therefore, we derived prescriptions for the fitting parameters of those radial modes correction empirical models as functions of log ∆ν, log T eff , log g and log κ which are quantities easily computed by 1D stellar evolution model. The next step will be to test our prescriptions against observed frequency spectrum in order to determine their degree of accuracy. We will focus on this in a furture work.

Finally, we only considered in this article the issue of structural effects. However, other effects such as non-adiabaticity effects may also play a non-negligible role in the propagation of acoustic waves in the surface layers. This will be studied in a forthcoming paper.

Fig. 1 .

 1 Fig. 1. Patched models represented in Kiel diagram (see Sect. 2.1). The metallicities are shape-and colour-coded and are [Fe/H] = {-1.0, -0.5, +0.0, +0.5}.

T

  H] = -0.5, ∆r = 5.81 × 10 7[Fe/H] = 0.0, ∆r = 6.43 × 10 7[Fe/H] = +0.5, ∆r = 9.49 × 10 7 H] = -0.5, ∆r = 1.84 × 10 8[Fe/H] = 0.0, ∆r = 2.53 × 10 8 eff 6500 K log g 4.00[Fe/H] = -1.0, ∆r = 1.39 × 10 8[Fe/H] = 0.0, ∆r = 1.77 × 10 8[Fe/H] = +0.5, ∆r = 2.04 × 10 8

Fig. 2 .

 2 Fig. 2. Frequency differences of PM vs. UPM for radial modes with frequencies less than the acoustic cut-off frequency ν ac = c/(4πH p ). Frequencies on the abscissa are normalized by the ν max of each models. Top panel: models J*: T eff 5900 K and log g = 4.0. Middle panel: models M*: T eff 5500 K and log g = 3.5. Bottom panel: models C*: T eff 6500 K and log g = 4.0.

Fig. 3 .

 3 Fig.3. Frequency differences scaled by the frequency, taken at ν max against a scaling relation given by Rosenthal & Christensen-Dalsgaard (1999) (top panel; with σ = 0.89) and a scaling relation given by Eq. (13) where powers are left free (bottom panel; with σ κ = 0.63).

Fig. 4 .

 4 Fig. 4. Parameters a (top panel), b panel) and a 3,BG2 (bottom panel) across the Kiel diagram from K08 and BG2. The symbols refer to [Fe/H] = -1.0 (diamond), [Fe/H] = -0.5 (square), [Fe/H] = +0.0 (circle) and [Fe/H] = +0.5 (triangle).

Fig. 5 .

 5 Fig. 5. Logarithms of K08 (first row), S15 (second row) and BG2 (third row; log |a 3,BG1 | behaves almost identically to log |a 3,BG2 |) fitting coefficients as a function of log z 1 (Eq. 13): log f = slope × log z 1 + offset. Purple dots (resp. yellow, green and red) corresponds to [Fe/H] = -1.0 (resp. -0.5, 0.0 and +0.5).

Fig. 7 .

 7 Fig. 7. Ratio Q n of a mode of degree n against frequency of the same mode normalized by ν max (only ν n /ν max > 0.75 for clarity), for the subgiant model Cm10 with T eff = 6503 K, log g = 4.0 and [Fe/H] = -1.0.Each colour corresponds to a degree . Yellow (resp. green) dots breaking from the general trend correspond to dipolar (resp. quadrupolar) mixed modes.

Table 1 .

 1 Characteristics of the 3D hydrodynamical models and of the UPM and PM.

	Model	T eff log g	T b	ν max	M	Age Y init	Z init α MLT	R PM ∆r/R PM Stage
		K		(×10 4 ) K	µHz (M )	(Myr)	%	%		(R ) (×10 -3 )
	Cp05	6443	4	3.20 1057	1.91	1341 24.2	4.14	1.61	2.29	1.3	MS
	Fp05	6188	4	2.71 1078	1.81	1820 24.2	4.14	1.68	2.23	0.89	MS
	Jp05	5877	4	2.17 1103	1.53	4129 24.2	4.16	1.74	2.05	0.67	SG
	Lp05	5431	4	2.03 1161	1.32	8361 24.2	4.14	1.74	1.89	0.45	SG
	Op05	4933	4	1.98 1205	2.12	3.808 24.2	4.14	1.75	2.41	0.16 PMS
	Am00 5776 4.44	1.55 3063	1.02	4628 24.9	1.35	1.65	1.01	0.21	MS
	Bm00 6730 4.25	7.89 1839	1.38	1306 24.9	1.35	1.68	1.46	1.4	MS
	Cm00 6490	4	2.79 1053	1.48	2397 24.9	1.35	1.65	2.01	1.3	MS
	Dm00 6434 4.25	2.76 1881	1.27	2562 24.9	1.35	1.67	1.4	0.77	MS
	Fm00	6231	4	2.30 1074	1.28	4352 24.9	1.35	1.69	1.87	0.96	SG
	Gm00 6182	3.5	3.63 341.7	1.93	1158 24.9	1.35	1.7	4.09	2.4	SG
	Hm00 6101 4.25	2.51 1923	1.13	5413 24.9	1.35	1.67	1.32	0.53	SG
	Im00	5861	4.5	2.31 3479	1.09	34.8 24.9	1.35	1.66 0.976	0.2 PMS
	Jm00	5936 3.99	2.18 1077	1.16	6753 24.9	1.35	1.74	1.8	0.51	SG
	Km00 5886	3.5	2.18 348.5	1.88	1249 24.9	1.35	1.68	4.05	1.5	SG
	Mm00 5436	3.5	2.03 364.4	2.21	1017 24.9	1.35	1.71	4.37	0.83	SG
	Em05 6227	4.5	2.22 3395	1.02	644.2 25.1	0.597	1.74 0.938	0.24	MS
	Fm05	6260	4	2.26 1079	1.06	6520 25.1	0.597	1.79	1.7	0.85	SG
	Gm05 6149	3.5	3.27 338.2	1.66	1389 25.1	0.597	1.76	3.82	2	SG
	Im05	5897	4.5	2.25 3502 0.902	4290 25.1	0.597	1.62 0.882	0.17	MS
	Jm05	5919	4	2.19 1091 0.949 10 130 25.1	0.597	1.76	1.62	0.52	SG
	Km05 5764 3.51	2.13 363.2	1.59	1623 25.1	0.597	1.74	3.66	1.5	SG
	Mm05 5463	3.5	2.01 359.4	1.63	1483 25.1	0.597	1.74	3.78	0.7	SG
	Bm10 6743 4.24	7.60 1801 0.776 12 030 25.3 0.0205	1.84	1.1	1.2	SG
	Am10 5771 4.44	1.53 3065 0.768	13.47 25.3 0.0205	1.26 0.875	0.32 PMS
	Cm10 6503	4	2.63 1052 0.765 13 770 25.3 0.0205	1.75	1.45	1.4	SG
	Fm10	6242	4	2.25 1073	1.17	3.639 25.3 0.0205	1.67	1.79	0.8 PMS
	Gm10 6213	3.5	3.14	341	1.03	5093 25.3 0.0275	1.7	2.99	3	SG
	Km10 5894	3.5	2.15 349.7 0.906	7858 25.3 0.0205	1.62	2.8	1.9	SG

Table 2 .

 2 Prescriptions for the fitting coefficients involve in the empirical relations K08, K08r, S15, BG1, BG2. The prescriptions are written: log c 0 = c 1 log ∆ν/∆ν + c 2 log T eff /T eff + c 3 log g/g + c 4 log κ/κ + c 5 . Our solar values are computed from model Am00. ∆ν = 137 [µHz], T eff = 5776 [K], g = 27511 [cm s -2 ], and κ = 0.415 [cm 2 g -1 ].

	Law	log c 0	c 1	c 2	c 3	c 4	c 5
	K08	log |a|	1.03	3.26 -1.75	0.655 -2.72
		log b	-0.185 -0.584 0.313 -0.117 0.289
	K08r log |a|	1.08	3.39 -1.82	0.683 -2.65
		log b	-0.387	-1.22 0.655 -0.246 0.647
	S15	log |α|	0.999	3.15 -1.69	0.635 -2.36
		log β	-0.477	-1.51 0.808 -0.303 0.787
	BG1 log |a 3 |	1.93	6.09 -3.26	1.22 -11.9
	BG2 log |a -1 |	2.13	6.72	-3.6	1.35	-12
		log |a 3 |	1.8	5.67 -3.04	1.14	-12
	Notes.

Acknowledgements. L.M. thanks Dr C. Pinçon for many interesting discussions. H.G.L. acknowledges financial support by the Sonderforschungsbereich SFB 881 "The Milky Way System" (sub-projects A4) of the German Research Foundation (DFG).

where ω = 2πν and ξ r are respectively the angular frequency and the radial displacement of the mode, c s is the sound speed and M and R are respectively the total mass and the total radius of the star.

The first term in parentheses in the expression of Kn v,c 2 is the inertia of a mode:

Using a first-order expansion, we can write the radial displacement as (e.g. [START_REF] Unno | Nonradial Oscillations of Stars[END_REF]):

where A is a constant and ζ is a phase factor. Further inserting Eq. (A.3) into Eq. (A.2) leads to

Averaging the cosine term gives 1 /2 and then simply

where the last line was obtained by replacing ξ 2 r by its expression. Finally, by simplifying this expression we obtain:

Eventually, inserting Eq. (A.8) into Eq. ( 3), approximating the cosine by 1 /2 as in (A.5) and using Eq. ( 19) from Rosenthal & Christensen-Dalsgaard (1999):

Appendix B: Fitting coefficients

In Table B.1 we gather the values of the fitting parameters introduced in Eqs. ( 14), ( 15), ( 16), and ( 18) used in order to perform the fits shown in Fig. 5 and toderive the coefficients for the prescriptions given in Table 2. Notes. Coefficients a, α, a 3,BG1 , a -1,BG2 and a 3,BG2 are all negative.

Appendix C: Root-mean frequency differences after correction

In Table C.1 we gather the values of the root-mean frequency differences ν max √ D/N of radial modes after correction for each empirical law shown in Fig. 6.


Julien Zollinger

B Rouat

J Guyon

S K Pillai

M Rappaz

Influence of Ir Additions and Icosahedral Short Range Order (ISRO

come

Introduction

Kurtuldu et al. [START_REF] Kurtuldu | [END_REF]2] have presented recently a new mechanism of nucleation in fcc alloys based on icosahedral quasicrystals (iQC) or clusters formation in the liquid phase and acting as pre-nuclei for the fcc phase. This so-called "iQC-mediated nucleation" mechanism was evidenced first in Al-20wt%Zn:Cr [2], and then in Au-12.5wt%Ag-12.5wt%Cu:Ir [START_REF] Kurtuldu | [END_REF]. It is still a conjecture but based on several strong experimental facts, mainly based on EBSD observations: (i) The addition of minute amounts of Cr in aluminum alloys or Ir in gold alloys drastically reduces the grain size with an abnormal fraction of twinned grain boundaries; (ii) Several assembly of nearest-neighbor (nn) grains exhibit multiple-twins (MT) relationship, with a common 110 direction; (iii) Up to 5 nn fcc grains have MT relationship compatible with the symmetry of a decahedron and up to 9 nn grains exhibit a MT symmetry compatible with an icosahedron, (iv) Some nn grains exhibit a MT relationship compatible with the facets of a 2 Materials and methods rapid solidification experiment has been described and characterized previously [16]. In this case, 4g of each alloy were melted by induction up to 1030 • C (100 • C superheat). After the temperature was stabilized, the liquid droplet was allowed to fall on a copper substrate cooled by water (about 5 cm free fall). This experiment allows to reach substantial nucleation undercooling with typical values ranging from 15 to 25 K [17]. The associated cooling rate Ṫ is in the range -10 4 K/s at the substrate/melt interface to -10 2 K/s at the top of the solidified droplet.

After standard metallographic preparation, the samples were electrolytically polished with a LectroPol-5 apparatus from Struers and prepared by ion polishing with a Gatan PECS II apparatus (4kV for 12h). Samples were examined by SEM and OM for imaging, Energy Dispersive Spectroscopy (EDS) for chemical analysis. Electron BackScatter Diffraction (EBSD) for phase identification and orientation was performed on a FEI Quanta 600 F FEG-SEM equipped with a Brucker Si-drift Quantax EDS detector and a Nordlys high-speed EBSD camera for low resolution maps and on a Zeiss Auriga 40 FEG-SEM equipped with an Oxford Instrument Symmetry camera for high resolution maps. The accuracy to identify twin grain boundaries in EBSD reconstructed microstructures was set to 5 deg. as in [START_REF] Kurtuldu | [END_REF]2]. For the samples melted and solidified in the DSC calorimeter, the whole surface of the samples was analyzed by EBSD. For the samples processed with FFD, this analysis was performed over a surface 4 times larger than that shown in the Results section (≈ 16 mm 2 ).

As will be shown in the next section, Cu-rich particles observed in the alloy with Ir have a characteristic size that is substantially lower than the size of Au-rich grains analyzed by EBSD. In order to have statistically significant microstructural information, an area of 4 mm 2 was characterized using the methodology of Jeanmaire et al. [18]. Using a 2000 magnification required to detect the smaller Cu-rich particles, the analyzed surface corresponds to about 2000 BSE-SEM images. With a software specifically developed for particle analysis, the sample stage can be driven and the sample completely analyzed automatically. For each detected Cu-rich particle, the spatial coordinates, mean equivalent diameter and surface fraction are determined by image analysis.

Thermodynamic calculations were performed with the ThermoCalc software [19]. The thermodynamic parameters for the Au-Cu, Ag-Au, Ag-Cu and Ag-Au-Cu systems were taken from Sundmann et al., Hassam et al. and Kusoffsky,respectively [20,21,22], and from Dinsdale for pure elements [23].

Results

3.1 iQC-mediated nucleation in Au-20.5Cu-4.5Ag at low Ṫ Figure 2 shows the microstructure of Au-20.5Cu-4.5Ag with and without Ir additions solidified at high cooling rate using the FFD experiment and observed using BSE-SEM. In all the figures related to FFD results, the thermal gradient goes from left to right, i.e., the FFD specimen. To detect more nn grains with such MT relationships, the tolerance has to be increased in order to account for subgrain misorientations induced by the 7.35 deg. aperture default angle associated with the decahedral configuration of 5 fcc tetrahedra.

Setting the tolerance to 7.5 deg., Figure 4(a) shows a set of 9 nn grains which re-produce the symmetry of an icosahedron when fcc tetrahedra are set on its triangular facets with the heteroepitaxy relationship found in [2], i.e., 111 fcc // 3-fold symmetry axis ico and 110 fcc ⊥ 2-fold symmetry axis ico . The grains identified by EBSD are displayed with false colors which are then used in the four 110 pole figures of Fig. 4(c)-(f). In (c), the grains 1, 2, 3, 4 and 5 have MT relationships outlined by the five {111} common planes (arcs of circle) with 110 common directions (surrounded by small circles or ellipses for near-twins). The five {111} twin planes have a common 110 direction which we label 110 c . The same occurs for grains 2, 3, 7, 8 and 9 in Fig. 4(d) with a common 110 d direction, for grains 5, 6, 7, 8 and 10 in Fig. 4(e) with a common 110 e direction, and finally for grains 1, 5, 7 and 9 (only 4 grains) in Fig. 4(f) with a common 110 f direction. This clearly shows that each of the four sets of 5 (or 4) nn grains whose pole figures are shown in (c)-(f) are in MT relationship with a common 110 direction correspond-ing to the 5-fold symmetry axis of a decahedron (assembly of 5 tetrahedra with 10 triangular facets). Moreover, these four sets of grains correspond to a single icosahedron template. If such is the case, the angles between the common 110 c-f directions should correspond to those between the 5-fold symmetry axes of the icosahedron. In an icosahedron, these angles are either 63.43 deg. (first vertex neighbors) or 116.57 deg. (second vertex neighbors). The mutual angles found between the four common 110 c-f directions of Fig. 4, as identified by the crossing of the {111} planes, are: 60.3, 66.2, 72.9, 78.2, 79.2 and 125.7 deg.. The comparison between the common 110 c-f directions of the sets of fcc grains and the 5-fold directions of the assumed icosahedron template is shown in Fig. 4(b). Considering the cumulated aperture default angle of fcc grains around a common 5-fold symmetry axis of an icosahedron (7.35 deg.), subgrain misorientations induced during growth and measurement uncertainty, this comparison is a fairly strong evidence that ISRO is again responsible of the high twinning frequency of grain boundaries in this specimen.

Nucleation, grain distribution and stability of fcc-Cu particles

Figure 5 shows a BSE-SEM image corresponding to Figure 3(b), i.e., Au-20.5Cu-4.5Ag with Ir solidified rapidly, and the corresponding density of Cu-rich particles (or dendrites when they are bigger) as a function of the distance from the copper substrate. The density of particles was measured over the surface using 0.2 mm steps. Again the vertical red dotted line corresponds to the position of the CET of the fcc-Au grain structure. The evolution of the number of Cu-rich particles decreases with increasing distance from the substrate, and thus decreasing cooling rate. At the top of the sample where the cooling rate is "only" about -100 K/s, 740 particles/mm 2 were observed.

In order to check the Cu-rich phase stability, different heat treatment were performed between 700 and 900 • C with 50 • C steps for 1h with a dynamic argon atmosphere and furnace cooled. At 700 and 750 • C, the microstructure does not evolve significantly and only a smoothing of the dendrites envelope was observed. Between 800 and 850 • C, the fraction of Cu particles start to decrease. It appears that they have partially melted with an increasing fraction of liquid, as could be seen from the pores formed near these remelted particles. Although the temperature is much lower than the melting point of pure copper, diffusion lengths on the order of 1 µm were measured. Therefore, local diffusion tends to equilibrate the compositions and thermodynamic calculations indicate that partial melting should occur in this alloy. At 850 • C, about 50% of the particles have disappeared by this local diffusion-remelting process. At 900 • C, part of the sample has melted and there is no Cu-rich-phase remaining after cooling.

Discussion

iQC-mediated nucleation of the Au-phase

It is recalled that the iQC-mediated nucleation mechanism of fcc metals proposed by Kurtuldu et al. [2] considers that 20 new fcc grains can form on the facets of an icosa-fortunately not mentioned in [START_REF] Kurtuldu | [END_REF], but is probably close to that of the FFD specimen since Au-12.5Cu-12.5Ag was cast in a copper mold. A stronger indication to validate the conjecture of iQC-mediated nucleation is that many nearest-neighbor grains are in a configuration of multiple twinning compatible with the icosahedral symmetry, as shown for example in Fig. 4.

In the columnar zone, there are very few twins, regardless whether Ir is added or not. It does not mean that the iQC-mediated mechanism did not operate at the onset of columnar growth near the copper chill in the specimen with Ir. But as growth proceeds, grain competition selects those having a hk crystallographic direction most closely aligned with the thermal gradient, kh being normally 100 for cubic metals. Because of the angles between the 3-fold facets of an icosahedron, only one of the potentially 20 fcc grains forming on an iQC can have such an orientation. Therefore, it can be assumed reasonably that columnar grains selected after some growth distance originate from separate iQC's and should not exhibit any twin, except those than can be detected within the preset accuracy among randomly oriented grains. (As a matter of fact, a perfect 100 texture cannot exhibit any twin due to the cubic symmetry).

Several questions nevertheless remain. Why does the fraction of twinned grain boundaries vary as a function of the distance from the copper chill (or cooling rate) in the FFD experiment (between 2 and 6% in Fig. 3(c)), while the average fraction measured over a larger surface is close to that of the DSC specimen (see Table 1)? And why this fraction of twins varies non-monotonously while the size of equiaxed grains remain fairly constant, with a slow gradual decrease from about 50 µm to about 35 µm? Kurtuldu et al. [START_REF] Kurtuldu | [END_REF] have shown that the grain size correlates better with the fraction of twinned grain boundaries than with the amount of Ir added to the melt. Such seems not to be the case here. Finally, why is the columnar zone more extended and the corresponding grain size larger when Ir is added to the melt, while the equiaxed zone is drastically grain refined? It is well-known in Al alloys that the addition of TiB 2 inoculant particles refines the grain size and at the same time reduces, or even eliminates the columnar zone, a mechanism very well explained by simple criteria such as derived by Hunt [25].

The answers are not simple since the final grain structure in solidification is a convolution between nucleation and growth kinetics acting concomitantly in a non-steady thermal field. Without Ir, these mechanisms are well understood because both the columnar and equiaxed grains clearly exhibit a dendritic structure (see Fig. 2(top)), for which the growth kinetics is fairly well known and about the same for both morphologies. Knowing the growth kinetics, so-called micro-macroscopic approaches allow to interpret the results and give access to the nucleation distribution [START_REF] Dantzig | Solidification[END_REF]. With Ir (Fig. 2(bottom)), not only the grain size in the equiaxed zone is drastically reduced, but the microstructure of both columnar and equiaxed grains is clearly no longer dendritic. It is in fact much coarser as compared with the dendrite arm spacing of the Ir-free specimen. In the equiaxed zone, the grains are globular-dendritic and this could be due to the small grain size. (It is to be noted that there is no growth kinetic model for such morphology besides the phase field method). For columnar grains, the microstructure is astonishing: it looks like small cells growing more or less in the direction of the thermal gradient.

It could be thought that the formation of Cu particles, which occurs first, is respon- sible for this change of morphology of the Au-rich columnar grains when Ir is added. Although these particles represent about 4% in volume, this does not influence much the composition of the liquid, and thus microsegregation when the Au-rich phase forms (i.e., about 2 wt% less Cu in the liquid). On the other hand, their density being fairly high (Fig. 5), they could possibly influence (or "drive") the growth direction of the Au-rich phase, as it has been shown for example in [START_REF] Gránásy | [END_REF]. However, Figure 6 clearly rules out this mechanism and shows that this change of morphology, and thus of growth kinetics, of columnar grains is directly induced by Ir addition.

Figure 6 shows the corresponding textures of the columnar zone of the FFD Au-20.5Cu-4.5Ag specimens, without (a) and with Ir (b). In these Inverse Pole Figure (IPF) texture representation, the thermal gradient direction, i.e., normal to the copper chill substrate, is represented with respect to the orientation of the grains composing the columnar zone. Considering the fcc symmetry, IPF are represented in a stereographic projection limited to the unit triangle 001 -101 -111 . The dendritic morphology of the Ir-free specimen has a very strong 100 texture, i.e., the columnar grains have one of their 100 directions aligned with the thermal gradient direction, as revealed by the high intensity (black) near 001 in Fig. 6(a). This is expected from the grain competition occurring within the columnar zone of "normal" fcc alloys, where 001 dendrites well aligned with the thermal gradient have the lowest velocity, and thus the smallest undercooling. With Ir (Fig. 2), the peculiar microstructure of the columnar grains exhibits a 111 texture (Fig. 6(b)). This was totally unexpected and clearly confirms that ISRO in the liquid not only contributes to nucleation via the iQC-mediated mechanism, but also modifies the growth morphology and texture via an attachment kinetics contribution. If one considers that icosahedral clusters in the liquid have to orient themselves in order to attach 3-fold symmetry facets, i.e., 3 atoms, to the growing fcc solid, this would explain why this favors the growth of (111) planes parallel to the isotherms.

This "ISRO-influence" on the growth kinetics and morphology of the fcc phase is most probably responsible of peculiar morphologies in Al-alloys previously reported.

For example, Henry et al [28] have observed the formation of "curly" 112 dendrites or the formation of fractal-type secondary arms in commercial Al-Cu-Mg alloys. Kurtuldu [29] has seen that the growth directions of trunks and arms in Al-Zn:Cr change with the solidification speed. Bedel et al [30] have observed the formation of 111 dendrites in small Al-Cu droplets rapidly solidified. At the light of our present results, all these observations start to make sense if one considers the attachment kinetics of icosahedral clusters, rather than individual atoms, at the interface of the fcc phase with the liquid.

A similar ISRO-influence on growth could also explain the peculiar arrangement of some of the MT equiaxed grains observed in Fig. 4(a). Some of them appear at several locations with an alternated sequence, e.g., grains 3 and 4 appear in a sequence 4-3-4'-3'-4"-3" around grain 2, or grains 5 and 6 in the sequence 6-5-6'-5'-6" around grains 7 and 8. This was not the case in Figure 6 of [START_REF] Kurtuldu | [END_REF]. We could first think that the complex 3D internal microstructure of the grains is such that grains appear as several patches in metallographic cross sections, as can be the case for well-developed fully dendritic specimens. But this looks quite improbable, especially with the fairly compact globular-cellular morphology observed in Fig. 2. A more plausible explanation could be again the influence of ISRO on growth, and not (or not only) on nucleation, as suggested by preliminary atomistic calculations [START_REF] Ding | Frontiers in solidification[END_REF]. If an iQC template can indeed lead to the formation of 20 MT grains, as suggested by Kurtuldu et al. [START_REF] Kurtuldu | [END_REF]2], the reverse heteroepitaxial relationship could also act during growth.

Considering that ISRO "motives" or "bricks" in the liquid phase can attach to a single growing fcc grain with the same heteroepitaxy relationship, this situation becomes similar to a local stacking fault during growth: the next (111) layer forming on the stacking fault transforms this defect into a twin. In the case of an ISRO motive attaching to the fcc grain, another (111) layer can form on another facet of the icosahedron motive leading to the formation of a twin but with a misorientation of 1.5 degree. In the present observations, this process can explain the alternance of twinned grains observed in the EBSD pictures. As far as the texture and morphology of columnar grains are concerned, it is assumed that the 111 texture and very peculiar morphology observed at high solidification rate with Ir addition is also due to this attachment of ISRO motives for two reasons: (i) the presence of ISRO decreases the mobility of atoms in the liquid phase, as measured by Kurtuldu et al. for Al-Zn:Cr [14], thus introducing an attachment kinetics contribution revealed at much lower solidification speed than for "ISRO-free" metals where this would occur at several meters per second; (ii) at the solidification speed of the droplet experiment, the attachment kinetics anisotropy already dominates the interfacial energy anisotropy, thus making the 111 direction the "easy" growth directions instead of the 110 and totally changing the overall dendritic morphology into some kind of cells.

Formation of Cu-rich particles

It has to be emphasized that the Cu-rich particles form from the liquid phase during the solidification of the droplet and not via a solid state transformation. Several arguments support this statement.

• First, the particles have a very dendritic morphology as can be seen in Figure 2 or in the detailed EBSD reconstructed map obtained after Ar-ions polishing shown in Fig. 7(a). Again, the twinned grain boundaries between the coarse Au-rich grains are white in this figure and the Cu-rich particles or dendrites are more visible in the zoom of the large brown Au grain shown in (b). The dendritic morphology of the Cu-rich phase is even more apparent in (d) and (e) where two of them have been enlarged. While dendrites can be observed during solid state transformation [START_REF] Lima | [END_REF], their formation during cooling of the solid and their large dendrite arm spacing are incompatible with the large cooling rate associated with the FFD experiment.

• Cu-rich dendrites sometimes go across a whole Au-rich grain and clearly appear as patches with the same color/orientation, e.g., in the light brown grain located at the bottom left of the figure. They sometimes even go across Au grain boundaries while keeping the same morphology and orientation, as shown for example in the neighbor pink Au grain. Finally, as revealed by the 110 pole figure of Fig. 7(c), there is no evident crystallographic orientation relationship between the Cu particles/dendrites and the surrounding Au grain. In this figure, the 110 crystallographic directions of the brown Au grain shown in (b) are identified with small black-brown circles. There are at least 8 or 9 Cu particles/dendrites in this single Au grain, with no clear orientation relationship. Therefore, this Au grain is definitely not the parent phase for the formation of these Cu particles/dendrites via a solid state transformation.

Interestingly, many of the Cu-rich dendrites exhibit a near-twin or twin relationship, as shown in Fig. 7(d) and (e), respectively. Kurtuldu et al [2] have defined a near-twin relationship as a twin, i.e., rotation of 60 deg. around a common 111 direction, followed by a rotation of about 7.5 deg. around a common 110 direction. This last rotation is associated with the default aperture angle of 7.5 deg, when setting 5 fcc tetrahedra in a decahedral symmetry around a common 110 direction. Although the level of iridium is too low to be detected by EDS, we suspect these atoms to be also responsible of this twinning mechanism via ISRO in the liquid. In Fig. 7(d), twinning seems to start from near the nucleation center of this twinned Cu-rich dendrite, while in (e) the purple and light brown dendrites look like individual entities who have "collided" during growth, an indication that twinning has occurred during growth.

The formation of copper particles is not only related to the presence of iridium in the alloy, but also on the solidification conditions associated with the FFD experiment, i.e. high undercooling and high thermal gradient. No particles were found in the DSC samples nor in alloys processed with FFD without iridium. These Cu-rich particles were also found to be metastable, since the Cu-rich phase can be eliminated by heat treatment.

While there are no known binary, ternary or quaternary phases and metastable phases in the Au-Cu-Ag-Ir system and sub-systems [START_REF] Igor L Shabalin | Ultra-High Temperature Materials I[END_REF], the three binary Au-Ir, Ag-Ir and Cu-Ir systems show a metastable miscibility gap. At high nucleation undercooling, and due to the small solubility of Ir in Au-based alloy (literature values for Ag, Au and Cu are reported in Table 2 at the invariant reaction temperature), the liquid phase can be brought to the miscibility gap region in FFD experiments. Without accurate thermodynamic description of the quaternary system, this issue can not be addressed quantitatively. It can nevertheless be stated that the liquid phase entering the miscibility gap region will undergo a spinodal decomposition process. Using an analogy to Guinier-Preston zones during precipitation, for example in Al-Cu alloys, Rappaz and Kurtuldu have also assumed a spinodal decomposition in the Al-Zn-Cr system to explain the formation of Cr-centered icosahedra in the liquid in a thermodynamic approach of nucleation [3]. In the present case, the formation of ISRO around iridium atoms would slow down diffusion, as shown by Kurtuldu et al in Al-Zn:Cr [14], allowing at high cooling rate the formation of Cu-rich and Au-rich liquid regions, both containing some amount of iridium. Unfortunately, the phase diagram of Au-Ir is totally uncertain (peritectic or eutectic according to [START_REF]Au-Ir binary phase diagram[END_REF] or [START_REF] Nielsen | [END_REF], respectively, with very different maximum solubilities), while that for Ag-Ir is eutectic-type and that for Cu-Ir is peritectic-type with the largest solubility (see Table 2). From the invariant temperature of the binary X-Ir phase diagrams, it seems natural that iQC-mediated nucleation would occur first in Cu-rich and then in Au-rich liquid regions. Cu-dendrites would then be stopped during their growth by the solidification of the Au-rich regions.

Conclusions

Processing Au-20.5wt%Cu-4.5wt%Ag alloy with 0.01 Ir addition under rapid solidification conditions has given new insight regarding ISRO and i-QC mediated nucleation. It has been shown that the mechanism previously described by Kurtuldu et al. [START_REF] Kurtuldu | [END_REF] still applies regardless of the copper and silver composition of the gold alloy. While the number of twins does not evolve significantly with increasing cooling rates, the microstructure is however strongly affected by Ir-induced ISRO of the liquid phase. At high solidification rate, additions of 100 ppm of Ir change the growth morphology from dendrite-to cellular-type as well as the preferred growth direction from standard 100 dendrites to a 111 texture during columnar growth. This correlates the orientation relationship of 111 fcc // i-3 axis observed for nucleation, indicating that ISRO also modifies the attachment kinetics.

The most surprising result presented in this paper is the nucleation of primary equiaxed copper particles/dendrites, which are even twinned, when Ir is added to the melt under high cooling rate. It appears that Ir-induced ISRO, by slowing down diffusion, favors a spinodal decomposition of the liquid, with formation of Cu-rich and Aurich liquid regions containing some iridium. Nucleation and growth in both regions is then influenced as well by ISRO since twinning of Cu-dendrites and multiple-twinning of Au-rich grains are observed.

List of Figures

Figure 6 :

 6 Figure 6: Inverse texture pole figures of the columnar zone of the FFD specimens without (a) and with 100 ppm Ir (b), shown in the unit triangle 001 -101 -111 of the fcc structure. In these stereographic projections, the grey level indicates the probability of having the thermal gradient direction aligned along a certain crystallographic direction of the Au-rich fcc grains. The texture density is scaled between 0 and 10 for both specimens.

Figure 7 :

 7 Figure 7: (a) EBSD false color map of the Au-rich grain structure and Cu-rich particles/dendrites obtained on a section of the FFD specimen with Ir after Ar-polishing. The twin grain boundaries between Au-rich grains are shown in white. A zoom of the large brown grain is shown in (b) with its internal Cu-rich particles/dendrites, while the corresponding 110 pole figure is shown in (c) with the same color code. The 110 directions of the Au-rich grain are outlined with small circles. Figures (d) and (e) show two Cu-rich dendrites at other locations of this section with their corresponding 110 pole figure showing that they are near-twinned and twinned, respectively.

6

  Inverse texture pole figures of the columnar zone of the FFD specimens without (a) and with 100 ppm Ir (b), shown in the unit triangle 001 -101 -111 of the fcc structure. In these stereographic projections, the grey level indicates the probability of having the thermal gradient direction aligned along a certain crystallographic direction of the Aurich fcc grains. The texture density is scaled between 0 and 10 for both specimens. . . . . . . . . . . . . . . . . . . . . . . . . . . .

Table 1 :

 1 . . . . . 13 7 (a) EBSD false color map of the Au-rich grain structure and Cu-rich particles/dendrites obtained on a section of the FFD specimen with Ir after Ar-polishing. The twin grain boundaries between Au-rich grains are shown in white. A zoom of the large brown grain is shown in (b) with its internal Cu-rich particles/dendrites, while the corresponding 110 pole figure is shown in (c) with the same color code. The 110 directions of the Au-rich grain are outlined with small circles. Figures (d) and (e) show two Cu-rich dendrites at other locations of this section with their corresponding 110 pole figure showing that they are neartwinned and twinned, respectively. . . . . . . . . . . . . . . . . . . . 15 Average twin boundary fraction and grain size in the equiaxed zone of the samples with/without Ir, solidified at low/high cooling rate.( only two grains were found in this sample). The error associated with the twin frequency is calculated with a 99% confidence interval, the error associated with the grain size corresponds to the standard deviation given by the EBSD measurements.

	20

Table 2 :

 2 Solubility values of Ir (at.%) in liquid and fcc solutions for the binary systems of interest at the invariant reaction temperature (eut.: eutectic, per.: peritectic). For comparison, the liquidus of Au-20.5Cu-4.5Ag is 900 • C.

Acknowledgements

The authors would like to thank Dr. S. Prades-Rödel and Dr. A. Blatter at PX Group for providing the gold alloys. J.Z. would like to thank E. Etienne for sample preparation. Y. Amal, E. Delon, D. Lacueva Oyarzabala and A. Masson are also acknowledged for their help with heat treatments and associated characterisations.


R Gantois
email: renaud.gantois@mines-albi.fr

A Cantarel
email: barthur.cantarel@iut-tarbes.fr

G Dusserre
email: cgilles.dusserre@mines-albi.frd

J.-N Félices

F Schmidt

Mold Filling Simulation of Resin Transfer Molding Combining BEM and Level Set Method

Keywords: Numerical simulation, Boundary Element Method, Level Set, Liquid Composite Molding, Darcy's law

Liquid Composite Molding (LCM) is a popular manufacturing process used in many industries. In Resin Transfer Molding (RTM), the liquid resin flows through the fibrous preform placed in a mold. Numerical simulation of the filling stage is a useful tool in mold design. In this paper the implemented method is based on coupling a Boundary Element Method (BEM) with a Level Set tracking. The present contribution is a two-dimensional approach, decoupled from kinetics, thermal analysis and reinforcement deformation occurring during the flow. Applications are presented and tested, including a flow close to industrial conditions.

Introduction

Liquid Composite Molding (LCM) is used in a wide range of industries (e.g. automotive, aeronautics, spatial, marine, defense, goods, etc.). Since the apparition of Resin Transfer Molding (RTM), a great number of variants has been developed in order to improve the forming conditions, lower production costs and reduce processing time. Among them can be mentioned Vaccuum Assisted Resin Transfer Molding (VARTM), Structural Reaction Injection Molding (SRIM) and Liquid Resin Infusion (LRI). The main feature in common is a liquid resin flow through the fibrous preform. A typical RTM setup [1] is depicted on Figure 1. The process can be divided into three major steps a) the dry plies are cut and stacked, b) the reinforcement is preformed and put in place in the mold, c) the mold is filled using the liquid resin, d) the part is cured and demolded.

In industry, a condition to meet the mechanical specifications is the entire impregnation of the preform. For years, engineers have worked using the trial and error method in order to optimize the filling conditions. An acceptable design is often reached after iterations because it involves the control of a large numbers of processing parameters. Considering the tool costs, it is clear that numerical simulation of flow is of a great interest, in particular for complex parts. Motivated by this need, numerical modeling has become an extensive field of investigation for the last three decade. Several models based on continuum mechanics and heat transfer has been studied at different scales. Simulations at fibre or tow scale are generally conducted using Stokes' flow.

Due to the complexity of the fibrous network, Representative Elementary Volumes are usually considered. Simulations at part scale are frequently adopted in industry considering an average flow modeled using Darcy's law. As the resin flows, a tracking technique is employed to follow the moving front in the mold. It can be performed using Lagrangian or Eulerian approaches. The first approach offers an interesting framework because the moving mesh explicitly describes the boundary [2][3][4]. Previous works [3] showed that the technique is highly suitable in combination with a Boundary Element Method. However, it becomes difficult to implement as soon as topological changes appears, which is frequently encountered in mold filling simulation. Eulerian approaches are clearly more consistent for that task. In such a framework, an auxiliary field is used to capture the front. In Volume Of Fluid (VOF) methods [5,6], it is defined at each cells as the fluid fraction. Front is implicitly located through connecting the partially filled cells. Some enhancements aim to approximate the front shape using reconstruction algorithms (e.g. Simple Line Interface Calculation (SLIC) or Piecewise Linear Interface Construction (PLIC)). As the field is discontinuous, a major drawback is that it requires a large number of cells to reach an acceptable level of accuracy, which tends to dramatically increase the computational time. In Level Set methods [7,8] the auxiliary field is a continuous function, usually chosen as the signed distance to the front. Thanks to the regularity of the field, an accurate shape of the front can be computed. The present paper focus on a technique combining a Boundary Element Method (BEM) [2,3,8,9] with a level set technique, following the approach developed in reference [8]. Level set theory indicates that the front motion is governed by the boundary velocity field. As a result the use of BEM together with level set techniques is straightforward to compute the filling pattern. Besides, the versatility of boundary mesh makes easier the remeshing procedure comparing to Finite Element Methods. The first part of this paper considers the governing equations. The second part describes the implemented model. The last part covers some applications.

Governing equations

Modeling Assumptions. Impregnation is modeled as a Darcy's flow. Fibrous preform is regarded as a porous medium. As a first approach, we focus on the flow front numerical modeling. Hydro-mechanical couplings are neglected, so that the medium is assumed to be not deformable and of a constant permeability. It should be noted that BEM technique allows coupling preform deformation and resin flowing but this is not addressed in this paper. The flow is considered as isothermal and resin viscosity is assumed to be constant, neglecting chemical reactions effects. Again, coupling heat balance equation including reactive source term remains possible but has not been investigated. Consequently, governing equations reduce to Darcy's law and incompressibity equation

[ ]      >= < ∇ ∇ - >= < (2) 0 . (1) ' v p K v ! ! ! ! µ where > < v ! is
     >= < ∇ ∇ - >= < (4) 0 . (3) ' v p K v e ! ! ! ! µ where e K is given as 2 1 K K K e = , 1
K and 2 K being the principal permeabilities. The isotropic equivalent transformation [3,10] is given by (5) Outline. The algorithm is implemented using Matlab environment. At the beginning, preprocessing imports a standard mesh file, performs the isotropic equivalent transformation, assigns material and processing data, and locates the injection ports. Next, the filling routine based on a Level Set formulation advances the front. At the end, a post-processing program is run in order to visualize the filling pattern once the real domain is recovered. The filling stage is divided into a finite number of quasi steady states. At each step time, governing equations are solved using a constant Boundary Element Method. The front is advanced by feeding a level set solver with the BEM-computed velocities. The process is repeated until the mold is completely filled. (10) Rewritten in a matrix form using H and G , Eq. 10 can be transformed into

1 0 0 1 2 1               =       y x K K K y

Boundary Element

∫ ∑ ∫ ∑ Γ = Γ = Γ = Γ + j j d p q d q p p c N j j N j j i i
ij 1 1 G q H p N j j ij N j j ∑ ∑ = = = (11)
where

∫ Γ Γ + = j d q H ij * ij
F AX = (12)
Where X is the vector of N unknowns, A is a 2 N matrix obtained by reordering H and G . F is the known vector computed from the boundary conditions, H and G matrix.

Level Set method. The moving front is captured using the signed-distance function φ , defined so that the zero level set corresponds to the interface [7,13]

{ } 0 ) , , ( / ) , ( ) ( 2 = ∈ = Γ t y x IR y x t φ (13)
For each point under consideration, distance is signed negative if located in the impregnated area and positive otherwise. Two types of mesh are utilized. The first one is a background grid. It is generated by a standard triangle element mesher using the CAD of the cavity. The other mesh is a computational boundary mesh. It bounds the fluid phase at each step time. Level set equations govern the evolution of the signed-distance function. To combine with BEM, we use a formulation involving a propagating interface with a velocity in its normal direction. It is given as follows

     = = = ∇ + ∂ ∂ (15) ) 0 , , ( (14)
0 0 φ φ φ φ t y x F t !
where F is the extended normal velocity at time t, built in order to coincide with the real velocity of the front (it is extrapolated outside without physical meaning). The filling stage is initialized by choosing 0 φ as the signed-distance to the inlet gates. Next, values for each point are updated using an Euler scheme

(16) ) , , ( ) , , ( t F t y x t t y x ∆ ∇ - = ∆ + φ φ φ ! Where t
∆ is the time step, adjusted to match its upper limit according to Courant-Friedrich-Levy (CFL) conditions. F and φ ∇ ! are computed at time t.

For the sake of efficiency, φ can be updated on few nodes around the front in a "narrow band" [7] instead of the entire grid. By the next, front is updated using a boundary element mesher developed during this work. It is based on interpolating the zero level set on the grid through a Marching Triangles algorithm.

It should be noted that the method directly handles topological changes, involving merging or dividing fronts. However it does not ensure that the resin remains inside the mold. Resin-mold contact is implemented using the fixed level set describing the mold walls.

Applications

Isotropic radial injection. The numerical method, was first validated by considering the academic case of an isotropic radial injection. The mold chosen for that study is a 0.25 m 2 square plate. The resin enters the part from the center and flows through the preform describing circular patterns. Four vents located on the corners ensure the complete filling. Material and processing parameters are given in Table 1. Table 1 Material and processing parameters

K [m 2 ] µ [Pa.s] ρ [kg.m -3 ] ε [-] 0 p [Pa] f p [Pa] 0 r [m]
1.10 -9 0.1 1150 0.5 2.10 5 1.10 5 2.5.10 -3

Neglecting gravity, the front motion in the free flow stage is given by the following analytical solution [3,10] (17)

) ( 4 1 1 ln 2 . 2 0 0 0 2 0 r t p p K r r r r f f f εµ - = +         -                
where f r and 0 r are the radii of the moving front and the inlet gate, 0 p and f p the inlet and outlet pressures. The mold is meshed using 2500 triangle elements. CPU time is around 20s on a 2.26 GHz / 1.93 Go of RAM laptop. The predicted filling time is 98s. Figure 2 shows the comparison between the implemented the analytical results at time 10.7s and 63.3s. The agreement is fair, with an error of 2.1%. Anisotropic injection involving complex shapes. The model was also compared with a more realistic case. A test was performed consisting of impregnating a 1x1 rib knit fabric made of glass fibers [14] using a canola oil of known viscosity. The preform was placed under a transparent flexible bag and filled using two inlet gates under vaccum conditions (infusion). Front is pulled using two vaccum lines on the upper and lower sides. The dimensions of the mold are given in Fig. 4. Some inserts were placed to assess the contact algorithm. Material and processing parameters are given in Table 2. CPU time is 200s for a model involving 3693 elements. Data was compared at 0.5s, 3s, 5s and 7s. Figure 6 indicates that agreement is fair at any time. In particular, a front merging is accurately predicted.

Conclusion

The presented method is based on combining BEM with Level Set method. Velocity computation is reduced to boundary, which tends to reduce the computational time. An important feature of the method is that the front shape is accurately computed which is more difficult to achieve in standard VOF techniques. Comparing with a Lagrangian approach, the Level Set front tracking facilitates the treatment of topological changes, typically appearing when fronts merge into each other.

Possible extensions of this work may include other physics occurring during the flow (e.g. hydromechanical couplings, kinetics). It should be noted that the implemented model includes a gravity term, and can be easily modified to take into account other body forces (inertia effect due to rotation for example). The technique remains similar in three dimensions, considering the resin flow as surfaces propagating.

  Figure 1 RTM process

(

  the physical coordinates.Finally, the resin velocity v ! is given by dividing the average velocity by the reinforcement porosity ε, which turns the system of Eqs. 3 and 4 into Neumann condition)where p (resp. q ) is the prescribed value of pressure (rep. normal derivative). In this paper, 0 = q on the mold walls (non penetration condition), p being the prescribed value of pressure at inlet and outlet.

2 N

 2 matrix. Finally Eq. 11 is reordered to take the form of a linear system as

Figure 2 Figure 3

 23 Figure 2 Filling patterns in isotropic radial injection: comparison with analytical solutionThe same simulation was conducted by activating the gravity contribution. A comparison between the implemented method and a FEM/VOF-based code (PAM-RTM™) is given Figure3. The same mesh was used for both simulations. As expected, the flow is slightly deported in the lower part of the plate. The predicted filling time is 101s for the implemented model and 89s for PAM-RTM™. The difference is less than 14%, which shows a good qualitative adequacy.

Figure 4 1 K [m 2 ] 2 KFigure 4

 4124 Figure 4 Dimensions of the mold

  and g denote respectively the resin specific mass and gravity. Combining the previous equations and transforming the coordinates in the isotropic equivalent domain e

	p	'	=	p	+	gz ρ
	where, ρ Ω of permeability

the macroscopic velocity, ] [ K the permeability tensor, µ the liquid resin viscosity and ' p the acting pressure. That pressure includes a gravity term. It is given by e K , we obtain

  Method. For clarity, the subscript e referring to the isotropic equivalent domain is omitted by the next. Laplace's equation is multiplied by the Green function and θ is the internal angle of the corner in radians. For a two-

	Where i p is the value of the pressure at a point i M located on the boundary, * q the pressure
	gradient associated with
	dimensional domain,	* p and * q are given as follows
	p	*	=	2 1 π	  	1 r	  	and	q	*	=	2 1 π	2 . n r r -	(9)
	where	r	=	M	i	P	, i M the point of application of the Dirac delta function and P the point under
	consideration. Meshing the boundary into N constant boundary elements and applying Eq. 8 leads
	to
											*								*
						1										1
					i c	i p	+	Γ ∫	pq	*	Γ d	=	∫ Γ	p	*	q	Γ d	with	i c	=	π θ 2	(8)

* p .

Integrated twice by parts over the calculation domain and using Green's theorem leads to Somigliana's equation

[3,11,12]

* p

Acknowledgments

This work was carried out in the framework of FUSCOMP project coordinated by DAHER-Socata. The authors gratefully acknowledge B. Cosson and S. Leroux for their contributions during this work.