Q Douasbin
email: quentin.douasbin@imft.fr

C Scalo

L Selle

T Poinsot

Delayed-time domain impedance boundary conditions (D-TDIBC)

Keywords:

Defining acoustically well-posed boundary conditions is one of the major numerical and theoretical challenges in compressible Navier-Stokes simulations. We present the novel Delayed-Time Domain Impedance Boundary Condition (D-TDIBC) technique developed to impose a time delay to acoustic wave reflection. Unlike previous similar TDIBC derivations (Fung and Ju, 2001-2004 [1,2], Scalo et al., 2015 [3] and Lin et al., 2016 [4]), D-TDIBC relies on the modeling of the reflection coefficient. An iterative fit is used to determine the model constants along with a low-pass filtering strategy to limit the model to the frequency range of interest. D-TDIBC can be used to truncate portions of the domain by introducing a time delay in the acoustic response of the boundary accounting for the travel time of inviscid planar acoustic waves in the truncated sections: it gives the opportunity to save computational resources and to study several geometries without the need to regenerate computational grids. The D-TDIBC method is applied here to time-delayed fully reflective conditions. D-TDIBC simulations of inviscid planar acoustic-wave propagating in truncated ducts demonstrate that the time delay is correctly reproduced, preserving wave amplitude and phase. A 2D thermoacoustically unstable combustion setup is used as a final test case: Direct Numerical Simulation (DNS) of an unstable laminar flame is performed using a reduced domain along with D-TDIBC to model the truncated portion. Results are in excellent agreement with the same calculation performed over the full domain. The unstable modes frequencies, amplitudes and shapes are accurately predicted. The results demonstrate that D-TDIBC offers a flexible and cost-effective approach for numerical investigations of problems in aeroacoustics and thermoacoustics.

Background

Direct Numerical Simulation (DNS) or Large Eddy Simulation (LES) have become standard approaches for high-fidelity simulations of unsteady compressible flows. In many cases, accurate flow predictions can be obtained only if the reflection of the acoustic waves at the boundaries is precisely defined and/or controlled [START_REF] Poinsot | Boundary conditions for direct simulations of compressible viscous flows[END_REF][START_REF] Polifke | Partially reflecting and non-reflecting boundary conditions for simulation of compressible viscous flow[END_REF][START_REF] Poinsot | Theoretical and Numerical Combustion[END_REF][START_REF] Lodato | Three-dimensional boundary conditions for direct and large-eddy simulation of compressible viscous flow[END_REF]. Impedance is a widely-used quantity to characterize the reflection of acoustic waves at boundaries; it is naturally defined in the frequency domain whereas DNS/LES are performed in the time domain. Numerical methods permitting the imposition of impedance boundary conditions (IBCs) in time domain solvers fall under the category of Time Domain Impedance Boundary Conditions (TDIBC). While existing TDIBC methods [START_REF] Fung | Broadband time-domain impedance models[END_REF][START_REF] Fung | Time-domain impedance boundary conditions for computational acoustics and aeroacoustics[END_REF][START_REF] Scalo | Compressible turbulent channel flow with impedance boundary conditions[END_REF][START_REF] Lin | High-fidelity simulation of a standing-wave thermoacoustic-piezoelectric engine[END_REF][START_REF] Christopher | Time-domain impedance boundary conditions for computational aeroacoustics[END_REF][START_REF] Özyörük | A time-domain implementation of surface acoustic impedance condition with and without flow[END_REF][START_REF] Sjoerd | Impedance models in time domain including the extended Helmholtz resonator model[END_REF][START_REF] Cotté | Time-domain impedance boundary conditions for simulations of outdoor sound propagation[END_REF][START_REF] Zhong | A controllable canonical form implementation of time domain impedance boundary conditions for broadband aeroacoustic computation[END_REF][START_REF] Jaensch | Hybrid CFD/low-order modeling of nonlinear thermoacoustic oscillations[END_REF] have demonstrated their ability to achieve such goal, one case has not yet be considered in the Computational Aeroacoustics (CAA) community: a pure time delay. Developing Delayed-Time Domain Impedance Boundary Conditions (D-TDIBC), i.e. a TDIBC allowing to impose a time delay, is the goal of the present work.

The proper implementation of boundary conditions is critical for the stability of the numerical solution. For example, directly imposing mean values of velocity or pressure at the boundary leads to complete reflection of acoustic waves, hence leading to an artificial increase in disturbance energy.

Acoustic waves are the manifestation of perturbations around a mean state: the pressure field, for instance, can be expressed as p(t) = p 0 + p ′ (t), where p 0 is the mean pressure and p ′ (t) is the acoustic perturbation. In the present work, the time harmonic behavior e +iωt is assumed, yielding the convention φ(t) = ℜ(φ(ω)e +iωt ) for a generic periodic function φ ′ (t) oscillating with angular frequency ω. The corresponding Fourier transform convention reads:

φ(ω) = +∞ -∞
φ(t)e -i ω t dt [START_REF] Fung | Broadband time-domain impedance models[END_REF] with F being the Fourier operator φ(ω) = F (φ ′ (t)), yielding the complex fluctuation amplitude φ(ω).

The specific impedance Z * (ω) is defined in the frequency domain as the non-dimensional ratio of complex pressure and the velocity component normal to the boundary [START_REF]Acoustics of Ducts and Mufflers[END_REF]:

Z * (ω) = 1 ρ 0 c 0 p(ω) u(ω) • n (2)
where ρ 0 and c 0 are the mean density and speed of sound, respectively. Z * (ω) is a frequency-dependent complex-valued function, which is commonly used to characterize the response of an acoustic element. It is unbounded, contrary to the reflection coefficient R and the wall softness S, which are defined as:

R(ω) = 1 -Z * (ω) 1 + Z * (ω)
, S(ω) = R(ω) + 1.

(

) 3
Unlike the specific impedance Z * , both the reflection coefficient R and the wall softness S are bounded in amplitude (for passive acoustic elements), making them the working variable of choice in computations. Fig. 1 illustrates the definition of the outgoing and ingoing acoustic waves, A out n and A in n , oriented in the same direction of and opposite to the outward normal n. Following this convention, the outgoing and ingoing characteristic waves are defined as:

A out n (ω) = u(ω) • n + p(ω) ρ 0 c 0 , A in n (ω) = u(ω) • n - p(ω) ρ 0 c 0 (4)
Combining Eqs. ( 2) and (3) yields:

R(ω) = A in n (ω) A out n (ω) (5)

S(ω) =

A in n (ω) + A out n (ω)

A out n (ω) (6)
Impedance, reflection coefficient and wall softness are defined in the frequency domain whereas DNS/LES are performed in the time domain. However, in order for the TDIBC to be physically admissible, the following three mathematical conditions must be fulfilled, as stated by Rienstra et al. [START_REF] Sjoerd | Impedance models in time domain including the extended Helmholtz resonator model[END_REF]:

1. Reality: The characteristic waves imposed in the time domain should be purely real. Such a property can be ensured in the frequency domain, for instance in the reflection coefficient, if:

R * (-ω) = R(ω) [START_REF] Poinsot | Theoretical and Numerical Combustion[END_REF] where the superscript * denotes the complex conjugate.

2. Causality: Any process should have a zero-response before it is triggered. Mathematically, such a condition is fulfilled if 1 :

ℜ(p i ) ≤ 0 for i ∈ [1; n poles ] (8)
where p i represents the ith complex pole of the reflection coefficient R(ω) prolonged to the Laplace plane and ℜ(•)

denotes the real part operator. 3. Passivity: The reflected wave amplitude should be lower or equal than the amplitude of the incident wave. This is fulfilled for:

|R(ω)| ≤ 1 (9)
Should one of these constrains be violated, the imposition of TDIBC is expected to contribute to numerical instabilities. It is hence imperative to bear in mind such mathematical constraints when modeling R(ω) (or S(ω)).

The time domain equivalent of R(ω) is obtained by use of the inverse Fourier transform R(t) = F -1 (R(ω)) and S(t) = F -1 (S(ω)). When dealing with the imposition of acoustic properties at a boundary condition, both the outgoing wave A out n (t) (evaluated from the computational domain) and the reflection coefficient R(ω) (respectively the softness coefficient) are known. However, the ingoing wave A in n (t) is unknown and has to be evaluated from the values of A out n (t) and R(ω)

(respectively S(ω)). This can be done by recasting Eq. ( 5) (respectively Eq. ( 6)) and taking its inverse Fourier transform. As the inverse Fourier transform of a product yields a convolution integral [START_REF] Dunford | Linear Operators[END_REF][START_REF] Rudin | Functional Analysis[END_REF], Eqs. ( 5) and ( 6) can be expressed in the time domain as follows:

A in n (t) = t 0 R(τ ) A out n (t -τ ) dτ (10) A in n (t) = -A out n (t) + t 0 S(τ ) A out n (t -τ ) dτ ( 11
)
where the bounds of the integrals in Eqs. [START_REF] Özyörük | A time-domain implementation of surface acoustic impedance condition with and without flow[END_REF] and [START_REF] Sjoerd | Impedance models in time domain including the extended Helmholtz resonator model[END_REF] These assumptions restrict the bounds of the convolution integrals as it is the case when using the causal inverse Laplace transform [START_REF] Priestley | Spectral Analysis and Time Series, Vols. 1 and 2 , Probability and Mathematical Statistics[END_REF][START_REF] Korn | Mathematical Handbook for Scientists and Engineers: Definitions, Theorems, and Formulas for Reference and Review[END_REF]. Non coincidentally, the conversion from the frequency domain to time domain as done in most TDIBC technique (see section 3) requires the formal prolongation of the support of the Fourier transform, restricted to a real-valued frequency ω, to the complex domain. This is step is hidden by an abuse of notation where the quantity iω is intended as complex and hence following the formalism of Laplace transforms (s = iω).

Although mathematically equivalent to Eqs. ( 5) and ( 6), Eqs. [START_REF] Özyörük | A time-domain implementation of surface acoustic impedance condition with and without flow[END_REF] and [START_REF] Sjoerd | Impedance models in time domain including the extended Helmholtz resonator model[END_REF] are not directly applicable in high-fidelity compressible solvers because the direct evaluation of the convolution integral requires a large amount of memory. Indeed, A out n (t) and R(t) need to be stored at every iteration and each boundary node or face. For simple geometries, Nottin [START_REF] Nottin | Overview on Research Activities for Numerical Prediction and Control of Combustion Instabilities in Ramjet Powered Missiles[END_REF][START_REF] Nottin | Développement de méthodes de prévision des instabilités de combustion dans les foyers prémélangés[END_REF] has solved Eq. ( 10) directly. In computationally intensive frameworks such as DNS/LES of turbulent flows, a more efficient compact-in-time method is required for the evaluation of A in n (t), as shown in recent papers on this topic [START_REF] Fung | Broadband time-domain impedance models[END_REF][START_REF] Fung | Time-domain impedance boundary conditions for computational acoustics and aeroacoustics[END_REF][START_REF] Scalo | Compressible turbulent channel flow with impedance boundary conditions[END_REF][START_REF] Lin | High-fidelity simulation of a standing-wave thermoacoustic-piezoelectric engine[END_REF][START_REF] Jaensch | Hybrid CFD/low-order modeling of nonlinear thermoacoustic oscillations[END_REF][START_REF] Kaess | Time-domain impedance boundary condition for compressible turbulent flows[END_REF][START_REF] Jonghoon Bin | Broadband impedance boundary conditions for the simulation of sound propagation in the time domain[END_REF][START_REF] Jaensch | On the robust, flexible and consistent implementation of time domain impedance boundary conditions for compressible flow simulations[END_REF][START_REF] Monteghetti | Design of broadband time-domain impedance boundary conditions using the oscillatorydiffusive representation of acoustical models[END_REF][START_REF] Lourier | Scale adaptive simulation of a thermoacoustic instability in a partially premixed lean swirl combustor[END_REF].

Advances in time domain impedance imposition in the electromagnetic community have driven novel acoustic boundary conditions for computational aeroacoustics [START_REF] Christopher | Time-domain impedance boundary conditions for computational aeroacoustics[END_REF][START_REF] Özyörük | A time-domain implementation of surface acoustic impedance condition with and without flow[END_REF][START_REF] Özyörük | Time-domain numerical simulation of a flow-impedance tube[END_REF][START_REF] Fung | Impedance and its time-domain extensions[END_REF]. Özyoruk et al. [START_REF] Özyörük | A time-domain implementation of surface acoustic impedance condition with and without flow[END_REF][START_REF] Özyörük | Time-domain numerical simulation of a flow-impedance tube[END_REF][START_REF] Ozyoruk | Time-domain calculation of sound propagation in lined ducts with sheared flows[END_REF] and Tam and Auriault [START_REF] Christopher | Time-domain impedance boundary conditions for computational aeroacoustics[END_REF] have first developed formulations of TDIBC based on the z-transform and a Laurent series development of the impedance in the frequency domain resulting in a time domain Ordinary Differential Equation (ODE) problem. Numerical instabilities were observed by Özyoruk et al. [START_REF] Özyörük | A time-domain implementation of surface acoustic impedance condition with and without flow[END_REF][START_REF] Özyörük | Time-domain numerical simulation of a flow-impedance tube[END_REF][START_REF] Ozyoruk | Time-domain calculation of sound propagation in lined ducts with sheared flows[END_REF] and discussed by Tam and Auriault [START_REF] Christopher | Time-domain impedance boundary conditions for computational aeroacoustics[END_REF]. Recent work on this topic suggests that these instabilities are due to the illposedness of the impedance boundary condition when coupled with Myers' slip flow velocity condition [START_REF] Sjoerd | Impedance models in time domain including the extended Helmholtz resonator model[END_REF][START_REF] Brambley | Time-domain implementation of an impedance boundary condition with boundary layer correction[END_REF].

Fung and Ju [START_REF] Fung | Broadband time-domain impedance models[END_REF][START_REF] Fung | Time-domain impedance boundary conditions for computational acoustics and aeroacoustics[END_REF] have proposed a causal formulation of TDIBC based on the wall softness S. As all Recursive Convolution methods (RC), Fung and Ju's method evaluates the convolution integral (Eq. ( 11)) in a low-storage manner. To avoid recomputing the whole convolution integral at all times (which would be increasingly costly in time) it is evaluated by [START_REF] Fung | Broadband time-domain impedance models[END_REF] computing the change of the convolution integral's value over a time step and (2) store the updated value. Such method bypasses the need to store the time history of S(t) and A out n (t) as only the value of the convolution integral calculated up to the previous time step needs to be stored. The method relies on partial fraction modeling in the frequency domain, as discussed in Sec. 3.1. This modeling relies in identifying the poles and residues of a complex function -the wall softness coefficient S in Fung and Ju's case. The update of the convolution integral is done by solving Eq. ( 11) over a time step so that a low-cost evaluation is possible using a trapezoidal quadrature rule. Fung and Ju's formulation has been applied Fig. 2. Typical case where the computational domain can be cut at x = x 1 . The propagation of acoustic waves in the truncated portion (between x 1 and x 2 ) is accounted for by the delayed reflection coefficient R τ (Eq. ( 12)).

to one-dimensional inviscid Euler calculations [START_REF] Fung | Broadband time-domain impedance models[END_REF][START_REF] Fung | Time-domain impedance boundary conditions for computational acoustics and aeroacoustics[END_REF] and, more recently, applied in fully compressible Navier-Stokes simulations [START_REF] Scalo | Compressible turbulent channel flow with impedance boundary conditions[END_REF][START_REF] Lin | High-fidelity simulation of a standing-wave thermoacoustic-piezoelectric engine[END_REF], with some modifications required to ensure numerical stability of the implementation, namely setting the phase parameter to zero (see section 6.3 in Lin et al. [START_REF] Lin | High-fidelity simulation of a standing-wave thermoacoustic-piezoelectric engine[END_REF]).

Another low-storage alternative to RC methods, is referred to in the literature as the Auxiliary Differential Equation (ADE) method. Such method also relies on partial fraction modeling but it involves the temporal integration of a set of first-order ordinary differential equations derived from the poles and residues equivalent to the convolution integrals in Eqs. [START_REF] Özyörük | A time-domain implementation of surface acoustic impedance condition with and without flow[END_REF] and [START_REF] Sjoerd | Impedance models in time domain including the extended Helmholtz resonator model[END_REF]. The major advantage of this method is that the same time-marching accuracy of the Navier-Stokes solver can be used for the TDIBC, thus ensuring consistency in the overall time accuracy order. A comparison between the ADE and the RC methods can be found in Dragna et al. [START_REF] Dragna | A generalized recursive convolution method for time-domain propagation in porous media[END_REF]. The ADE method was also used to implement a TDIBC in Troian et al. [START_REF] Troian | Broadband liner impedance eduction for multimodal acoustic propagation in the presence of a mean flow[END_REF].

In the combustion community, a state-space approach, called Characteristic Based state-space Boundary Condition (CB-SBC), has been developed by Schuermans et al. [START_REF] Schuermans | Thermoacoustic analysis of gas turbine combustion systems using unsteady CFD[END_REF], Kaess et al. [START_REF] Kaess | Time-domain impedance boundary condition for compressible turbulent flows[END_REF] and Jaensch et al. [START_REF] Jaensch | On the robust, flexible and consistent implementation of time domain impedance boundary conditions for compressible flow simulations[END_REF], and used in Navier-Stokes simulations [START_REF] Jaensch | Hybrid CFD/low-order modeling of nonlinear thermoacoustic oscillations[END_REF][START_REF] Kaess | Time-domain impedance boundary condition for compressible turbulent flows[END_REF][START_REF] Jaensch | On the robust, flexible and consistent implementation of time domain impedance boundary conditions for compressible flow simulations[END_REF][START_REF] Lourier | Scale adaptive simulation of a thermoacoustic instability in a partially premixed lean swirl combustor[END_REF]. CBSBC relies on control theory of linear time-invariant (LTI) systems to build a reflection model through first-order differential matrix equations [START_REF] Jaensch | On the robust, flexible and consistent implementation of time domain impedance boundary conditions for compressible flow simulations[END_REF]. Two formulations of CBSBC were proposed. The first formulation of CBSBC is a state-space approach under the so-called controllable canonical form, i.e. based on the modeling of a transfer function as a rational polynomial. The same modeling procedure was used by Zhong et al. [START_REF] Zhong | A controllable canonical form implementation of time domain impedance boundary conditions for broadband aeroacoustic computation[END_REF] to model the impedance. The transfer function represents the broadband reflection coefficient. Similarly to Fung and Ju's method [START_REF] Fung | Broadband time-domain impedance models[END_REF][START_REF] Fung | Time-domain impedance boundary conditions for computational acoustics and aeroacoustics[END_REF], this formulation of CSCBC aims at imposing wave reflection with near-zero time delay as, for practical reasons, it fails to handle pure time delay. Indeed, for the pure time delay imposition problem, CSCBC modeling requires a high number of Padé polynomial coefficients leading to ill-conditioned matrices, making this formulation difficult to use. A second formulation of CBSBC was proposed by Jaensch et al. [START_REF] Jaensch | Hybrid CFD/low-order modeling of nonlinear thermoacoustic oscillations[END_REF] specifically to impose pure time delays. This approach consists in implementing a Linearized Euler Equation (LEE) solver for each of the impedance boundary conditions in the domain and performing the temporal integration using a first-order upwind scheme. The 1D spatial discretization needed by the LEE solver is done by storing in memory a state matrix of dimension 2N × 2N, and three matrices of dimension 2 × N, where N is the number of nodes considered in the spatial discretization. In their work, Jaensch et al. [START_REF] Jaensch | Hybrid CFD/low-order modeling of nonlinear thermoacoustic oscillations[END_REF] have shown that to impose a pure time delay for a 1D wave propagation problem it was necessary to consider a number of nodes as high as N = 1000 in order to avoid acoustic energy dissipation.

In this work, a reformulation of Fung and Ju's method based on the reflection coefficient is presented and used to impose pure time delays in DNS/LES. Although mathematically equivalent to the TDIBC multi-oscillator strategy used by Lin et al. [START_REF] Lin | High-fidelity simulation of a standing-wave thermoacoustic-piezoelectric engine[END_REF], which targets the wall softness, the formulation adopted in the present work focuses on the reflection coefficient as it is a more natural and convenient choice for problems where delayed temporal responses of acoustic boundaries are of interest.

Section 2 will present the definition of delayed reflection, which allows truncating long portions of the computational domains. Section 3 will introduce the TDIBC formulation used in this work and the fitting technique necessary to model a delayed reflection coefficient. In Sec. 4 D-TDIBC will be used to impose a time delay in a 1D wave propagation test case. Finally, in Sec. 5 a 2D thermoacoustically unstable combustion chamber is studied via DNS, using the D-TDIBC technique to reduce the domain size.

Problem formulation: delayed impedance boundary conditions

In many simulations of industrial systems, longitudinal waves propagate over elongated portions of the computational domain (e.g. in the exhaust pipe of a car engine, in the chimney on an industrial furnace or, in general, for ducts of both constant section and speed of sound), merely resulting in a time lag of the acoustic response. In such cases the computational domain can be significantly reduced, by replacing such elongated portions with a boundary condition accounting for the temporal delay associated with the acoustic wave propagation time.

Fig. 2 shows an acoustic problem in which the portion from x 1 to x 2 is truncated (grey area). In this computational setup, the boundary condition located in x = x 1 must account for the acoustic properties of the truncated domain. This is possible if the boundary condition at x = x 1 imposes a time delay to the wave reflection located at x = x 2 : the resulting reflection coefficient is called delayed reflection coefficient [START_REF] Poinsot | Theoretical and Numerical Combustion[END_REF][START_REF] Munjal | Theory of a two source-location method for direct experimental evaluation of the four pole parameters of an aeroacoustic element[END_REF]. Assuming planar and linear wave propagation yields the following expression for the delayed reflection coefficient R τ : where τ is the total time delay (i.e. the time required by a wave to travel from x 1 to x 2 and back) and i = √ -1 is the imaginary unit. In this formulation of the delayed reflection coefficient, it is assumed that the base pressure in the truncated section is uniform. An different formulation should be derived for cases where, for example, significant pressure drops are present in the truncated section. In this work we build upon the multi-oscillator TDIBC methodology of Fung and Ju [START_REF] Fung | Time-domain impedance boundary conditions for computational acoustics and aeroacoustics[END_REF], with the phase-parameter correction by Lin et al. [START_REF] Lin | High-fidelity simulation of a standing-wave thermoacoustic-piezoelectric engine[END_REF], to create a novel TDIBC formulation, called Delayed-Time Domain Impedance Boundary Condition (D-TDIBC), relying on the determination of the partial fraction representation (via poles and residues) of the delayed complex reflection coefficient, as outlined in Sec. 3.1. It is important to note that by adopting Fung and Ju's formalism, which focuses on the causal reconstruction of the softness (or absorption) coefficient, it is indeed appropriate to talk about a multi-oscillator fitting strategy, where the fitted softness is the linear superposition of the softnesses of various causal oscillators, which can be each representative of real acoustically absorptive panels; on the other hand, the decomposition of the reflection coefficient in partial fractions as done in the present work (see Eq. ( 13)) is not tied to a particular physical interpretation. Due to the peculiar mathematical nature of the delayed complex reflection coefficient (Eq. ( 12), plotted in Fig. 3), exhibiting periodic variations in the frequency domain, the proposed D-TDIBC strategy was developed to handle a large number of poles and residues. While the focus herein is on the reflection coefficient R, the same fitting methodology outlined in Sec. 3 can be extended to the wall softness S.

R τ = R(ω)e -iωτ with τ = 2(x 2 -x 1 )/c 0 (12)

Methodology

Time domain imposition of complex reflection coefficient

The direct numerical integration of Eq. ( 10) does, in principle, allow to impose a complex reflection coefficient in the time domain; however, this results in rapidly unsustainable memory and CPU costs, especially in DNS/LES simulations. Building upon the methodologies of Fung and Ju [START_REF] Fung | Broadband time-domain impedance models[END_REF][START_REF] Fung | Time-domain impedance boundary conditions for computational acoustics and aeroacoustics[END_REF], Scalo et al. [START_REF] Scalo | Compressible turbulent channel flow with impedance boundary conditions[END_REF] and Lin et al. [START_REF] Lin | High-fidelity simulation of a standing-wave thermoacoustic-piezoelectric engine[END_REF], we approximate a generic target reflection coefficient R(ω) with the reflection coefficient of the boundary condition R BC expressed as a sum of rational functions: [START_REF] Zhong | A controllable canonical form implementation of time domain impedance boundary conditions for broadband aeroacoustic computation[END_REF] where R BC (ω, n 0 ) is the reflection coefficient model and the poles and residues (p k , µ k ) must come as n 0 conjugate pairs:

R(ω) ≃ R BC (ω, n 0 ) = 2n 0 k=1 µ k iω -p k
p 2k = p * 2k-1 ; µ 2k = µ * 2k-1 ; k ∈ [1; n 0 ] (14)
These conjugate pairs ensure that the reflection coefficient R is real-valued in the time domain. The time domain equivalent of Eq. ( 13) is:

R BC (t, n 0 ) = 2n 0 k=1
µ k e p k t H(t) [START_REF]Acoustics of Ducts and Mufflers[END_REF] where the Heaviside function H(t) indicates that the causality constrained is fulfilled. The time-domain incoming wave value A in n (t) corresponding to R BC (ω, n 0 ) is obtained by substituting Eq. ( 15) into Eq. ( 10) yields

A in n (t) = 2n 0 k=1 t 0 µ k e p k τ A out n (t -τ )dτ = 2n 0 k=1 I k (t) (16)
comprising 2n 0 convolution integrals called I k , which can be split into two contributions [START_REF] Fung | Broadband time-domain impedance models[END_REF][START_REF] Fung | Time-domain impedance boundary conditions for computational acoustics and aeroacoustics[END_REF]:

I k (t) = I k (t -t)e p k t + µ k t t-t e p k τ A out n (t -τ )dτ (17)
Equation ( 17) shows that each I k (t) can be evaluated recursively: the first term corresponds to the temporal integration over the interval τ ∈ [0, tt] and the second term to the temporal integration over the interval τ ∈ [t - and A out n (t). The resulting CPU and memory cost does not increase with the number of time steps, and the computational overhead associated with the integral in Eq. ( 17) over a time step t is minimal. We adopt a trapezoidal quadrature rule for the discretization of the integral in Eq. ( 17), yielding [START_REF] Fung | Broadband time-domain impedance models[END_REF]:

I k (t) = I k (t -t)e p k t + α k A out n (t) + β k A out n (t -t) (18) α k = µ k e p k t -1 p 2 k t - 1 p k ; β k = µ k e p k t -1 p 2 k t - e p k t p k ( 19
)
It should be noted that Eqs. ( 18) and ( 19) limit the overall temporal accuracy order to second order; this can be overcome by developing an alternative time integration strategy based on the Auxiliary Differential Equation (ADE) method (see Background), which is out of the scope of the present paper.

The initial values of I k (t) and A out n are set to zero in the simulations [START_REF] Fung | Broadband time-domain impedance models[END_REF][START_REF] Fung | Time-domain impedance boundary conditions for computational acoustics and aeroacoustics[END_REF][START_REF] Scalo | Compressible turbulent channel flow with impedance boundary conditions[END_REF][START_REF] Lin | High-fidelity simulation of a standing-wave thermoacoustic-piezoelectric engine[END_REF]. However, one can make a checkpoint restart by initializing the I k (t) using stored values from a previous simulation. Using Eqs. ( 18) and [START_REF] Korn | Mathematical Handbook for Scientists and Engineers: Definitions, Theorems, and Formulas for Reference and Review[END_REF], complex impedances can be imposed in DNS/LES solvers.

The methodology to determine the model inputs (µ k , p k ) is now presented.

Multi-Pole Modeling

As discussed in Sec. 3.1, in order to impose the reflection coefficient R(ω) we need to determine the set of poles and residues (p k , µ k ). There is no evidence of the uniqueness of the solution for a set of (µ k , p k ). We propose here a method to determine the poles and residues that is viable when modeling delayed reflection coefficients (Eq. ( 12)). It is a tedious task as the delayed reflection coefficient R τ (ω) has many zeros introduced by the complex exponential e -iωτ . The method consists in an iterative least-square fit. At each iteration, a rational fraction is added to the model until the model reflection coefficient R BC (ω) converges to the desired reflection coefficient R(ω).

In order to develop this fitting method, we will first describe the properties of the rational fractions in Eq. ( 13).

Pole Base Function properties

The reflection coefficient of the boundary condition R BC (ω, n 0 ) is a sum of n 0 rational fractions, called the Pole Base

Functions (PBF), denoted R k : R BC (ω, n 0 ) = n 0 k=1 R k (ω) = n 0 k=1 µ k iω -p k + µ * k iω -p * k Pole Base Function ( 20
)
Equation ( 20) can be recast by expressing the poles and residues in their algebraic forms, i.e. as µ k = a k + ib k and p k =

c k + id k with a k , b k , c k , d k ∈ R. Each R k is causal if c k < 0 [11,35].
Although Fung & Ju [START_REF] Fung | Time-domain impedance boundary conditions for computational acoustics and aeroacoustics[END_REF] suggested that the phase parameter B k = b k d k + a k c k could be not null, Lin et al. [START_REF] Lin | High-fidelity simulation of a standing-wave thermoacoustic-piezoelectric engine[END_REF] have experienced spurious oscillations, leading to unstable numerical simulations, when the parameter B k was unconstrained B k = 0. Although the TDIBC of interest is based on the reflection coefficient R rather than on the wall softness coefficient S, similar behavior was observed in preliminary numerical trials by the authors. Thus, the phase parameter was constrained to B k = 0, i.e. b k = -a k c k /d k , leading to: 21)). Three properties are visualized: ω 0,k , h 0,k and ω ǫ . At the resonant angular frequency ω 0,k , the real part of R k is maximum. The peak height h 0,k is the value of R k (ω) at ω = ω 0,k . For a given percentage ǫ (ǫ ∈ [0, 1]), the half-width ω ǫ is defined on the right-hand-side. ω -and ω + are the solutions of the equation ℜ (R k ) = ǫ and are used to define the width.

R k (ω) = 2a k iω -ω 2 -2c k iω + c 2 k + d 2 k ( 21
)
There are three degrees of freedom in Eq. ( 21), namely a k , c k and d k . The objective of this section is to link these degrees of freedom to the properties of the PBF.

A typical PBF is presented in Fig. 4 where three properties of the PBF can be identified: ω 0,k , h 0,k and ω ǫ . ω 0,k is the frequency where the real part of R k is maximum: it is referred as the resonant frequency. h 0,k is the peak height, that is, the value of R k (ω) at ω = ω 0,k . The last property shown in Fig. 4 is the half-width ω ǫ defined for a given percentage ǫ of h 0,k (ǫ ∈ [0, 1]). Because the PBF is not symmetric, at a given percentage ǫ, two half-widths ω ǫ can be defined: one on the left-hand side (using ω -) and one on the right-hand side (using ω + ). A choice is made to consider only the right half-width shown in Fig. 4, hence overestimating the width.

One can show that [START_REF] Scalo | Compressible turbulent channel flow with impedance boundary conditions[END_REF][START_REF] Lin | High-fidelity simulation of a standing-wave thermoacoustic-piezoelectric engine[END_REF]:

ω 2 0,k = c 2 k + d 2 k ( 22
)
Equation ( 22) gives a first constraint, say on d k . Injecting Eq. ( 22) into Eq. ( 21) and evaluating the resulting equation at

ω = ω 0,k yields: R k (ω 0,k ) = h 0,k = - a k c k ( 23
)
Using the constraints in Eqs. ( 22) and ( 23), R k becomes:

R k (ω) = 2h 0,k c k iω ω 2 + 2c k iω -ω 2 0,k (24)
In Eq. ( 24), c k is the only remaining degree of freedom. Fig. 5(a) illustrates its effect on the width of a PBF R k for a peak height h 0,k = 1 at the resonant frequency f 0,k = ω 0,k /2π = 100 Hz. In the case of a pure delay the width is known: R(ω) is a periodic function of period T = 2π /τ (Fig. 5(b)). As shown in Fig. 5(a), the real part of R k goes to zeros for ω << ω 0,k and ω >> ω 0,k . A consequence of this property is that a single PBF can fit a single peak between two consecutive zeros.

For the pure delay the frequency range between two consecutive zeros is T /2. To adjust the width to the pure delay case, we need to find the value of c k for a given ǫ such that:

ℜ R k ω 0,k + T 4 = ℜ R k ω 0,k + π 2τ = ǫ (25)
where ℜ(•) and ℑ(•) are the real and imaginary parts, respectively. Solving for c k in Eq. ( 25), we obtain the criterion for a given time delay τ :

c k (ǫ, ω 0,k ) = - τ ω 0,k + π 4 2π ω 0,k + τ ǫ 1 -ǫ (26)
Equations ( 22), ( 23) and [START_REF] Lourier | Scale adaptive simulation of a thermoacoustic instability in a partially premixed lean swirl combustor[END_REF] give constraints allowing to control the height and the width of a PBF for any resonant angular frequency ω 0,k . These constraints are used in the algorithm allowing the determination of the set of poles and residues (p k , µ k ) presented in the next section.

Iterative Multi-Pole Modeling Technique

Now that the properties of a single PBF are known, it is possible to go back to Eq. ( 20) where a sum of n 0 PBFs are used to match a reflection coefficient R(ω). In the case of pure delay, it is a tedious task as the complex exponential function, introduced by the time delay (see Eq. ( 12)), has a large number of zeros. A common practice to obtain the multi-pole approximation of a complex function is to use the Vector Fitting method [START_REF] Gustavsen | Rational approximation of frequency domain responses by vector fitting[END_REF][START_REF] Gustavsen | Improving the pole relocating properties of vector fitting[END_REF]. It was used to model TDIBC using partial fractions by several authors [START_REF] Cotté | Time-domain impedance boundary conditions for simulations of outdoor sound propagation[END_REF][START_REF] Dragna | A generalized recursive convolution method for time-domain propagation in porous media[END_REF][START_REF] Reymen | Time-domain impedance formulation based on recursive convolution[END_REF][START_REF] Reymen | Time-domain impedance formulation suited for broadband simulations[END_REF][START_REF] Li | Improved multipole broadband time-domain impedance boundary condition[END_REF]. However, preliminary tests have shown only limited results when applied to complex exponentials as convergence could not be reached when modeling pure time delays. To tackle this issue an iterative fitting technique has been devised.

The least-square fit algorithm is used to minimize the distance of the target function R(ω) to the fit function R BC (ω) (Eq. ( 21)) on both real and imaginary parts. The global least-square residual ξ is defined as the sum of the squared distance between R(ω) and R BC (ω):

ξ = ξ R + ξ I ( 27
)
where ξ R and ξ I are the least-square residuals on the real and imaginary parts, respectively. They are defined as the squared sum of the point-to-point oriented distance, E R and E I , between R(ω) and R BC (ω):

ξ R = m i=1 E 2 R (ω i ) = m i=1 ℜ R(ω i ) -R BC (ω i ) 2 (28) ξ I = m i=1 E 2 I (ω i ) = m i=1 ℑ R(ω i ) -R BC (ω i ) 2 ( 29
)
where m is the number of discrete values taken into account in the frequency array, i is the index of the points, R(ω i ) is the value of the target reflection coefficient at the ith point of the discrete angular frequency array ω i and R BC (ω i ) is the value of the fit function at the ith point.

The fitting procedure is presented in the Algorithm 1. It starts with only one term: n 0 = 1 in Eq. ( 20) (line 1 in the Algorithm 1). The conditional loop iterates until the number of PBF n reaches the final value n = n 0 . In line 3, we seek for the frequency where the error on the real part is maximum. The resonant frequency of R n is then set to this frequency (line 4) and its peak height h 0,n is chosen to cancel the error on the real part at that point (line 5). The c n parameter is initialized using ǫ = 1% in Eq. [START_REF] Lourier | Scale adaptive simulation of a thermoacoustic instability in a partially premixed lean swirl combustor[END_REF]. Equations ( 22) and ( 23) are used to determine the values of d n and a n (line 7 and 8). Finally, the pole and residue (p n , µ n ) of the nth PBF are initialized (lines 9 and 10). The optimization stage (line 11) minimizes the least-square residual ξ (see Eq. ( 27)) and all of the values of the parameters p k and µ k are optimized for k ∈ [1; n]. Finally, the order of the model n is incremented so that an additional PBF can be added (line 12).

An example of iterative fit is presented in Sec. 4.2 where the model R BC (ω) is shown at several iterations of the Algorithm 1.

Neither the uniqueness of the solution for a set of (p k , µ k ) nor the convergence of the Algorithm 1 have been mathematically proven. Such derivations are beyond the scope of this work and require further studies. In practice, the convergence is left at the mercy and robustness of the least square algorithm, which is very reliable, and convergence was always reached for the delayed reflection coefficients considered in this study.

The coefficients 1/c k and 1/d k are time constants, and, as such, they might constrain the maximum allowed time step size of the overall computation due to numerical stability if either |c k t| or |d k t| is too large [START_REF] Cotté | Time-domain impedance boundary conditions for simulations of outdoor sound propagation[END_REF]. In the proposed fitting method, there is no explicit constraint for such parameters. However, in practice these methods are used in explicit DNS and LES codes where t is very small so that the no numerical stability issues were encountered by the authors while using poles and residues provided by Algorithm 1.

Algorithm 1: Preprocessing for D-TDIBC: find (p k , µ k ) in Eq. ( 20) to model the target function R(ω).

Input : n 0 , ω i , R(ω i ), τ Output: (p k , µ k ) for k ∈ [1, n 0 ] 1 n ← 1 2 while n < n 0 do 3
Find ω max such that the error on the real part is maximum

≡ {ω max | ∀ω i : E R (ω i ) ≤ E R (ω max )} 4
Set ω 0,n = ω max using Eq. ( 22)

5

Set the resonant peak height to cancel the error on the real part h 0,n = E R (ω max ) using Eq. ( 23)

6
Set the parameter c n using Eq. ( 26)

7
Set the parameter d n using Eq. ( 22)

8
Set the parameter a n using Eq. ( 23)

9 µ n ← a n + i(-a n c n /d n ) 10 p n ← c n + id n 11
Least-square fit: modify the values of p k and µ k for k ∈ [1, n] to minimize ξ (see Eq. ( 27))

12 n ← n + 1 13 end

Validation for purely reflective time-delayed impedance

Section 3 has provided the modeling methodology necessary to account for acoustic delays. In this section, the objective is twofold: (1) to demonstrate the applicability of the modeling procedure on a limit case: the delayed pure reflection, (2) to validate the ability of TDIBC to impose a time delay to account for acoustic wave propagation in the truncated portion of the domain. In Sec. 4.1 we present the one-dimensional numerical setup used for validation. The methodology proposed in Sec. 3 will be used to model a reflection coefficient corresponding to a pure delay. In Sec. 4.3 a time domain simulation is used to demonstrate that D-TDIBC imposes the correct time delay τ , wave amplitude and phase.

Test case presentation

The aim of the simulation is to study the propagation of a Gaussian acoustic wave in the domain defined on

x ∈ [-1, 1.75] m (Fig. 6 to travel back at its initial position is T = 2(x rx 0 )/c 0 = 10.0 ms.

Delayed impedance modeling

To impose a pure time delay, the D-TDIBC method first requires to build a model for the delayed reflection coefficient using the Algorithm 1 (see Sec. 3). In theory, the delayed reflection coefficient R r | x=1 m = -e -iωτ must be modeled up to infinite frequencies. In practice, it is not feasible as it would require an infinite number of PBFs. In numerous applications, the frequencies of interest are limited to a given range (Fig. 6(b)) and it is sufficient to model R(ω) below a cutoff frequency f c . For f > f c , a non-reflecting boundary condition is used. For example, for the setup presented in Fig. 6(a), the frequency content of the initial Gaussian wave is below 1 kHz so that the chosen cutoff frequency is f c = 1 kHz. A low-pass filter is applied to the theoretical delayed reflection by multiplying the delayed reflection coefficient by an envelope function ψ allowing a smooth transition from 1 to 0 over a frequency range:

ψ( f ) = 1 2 1 -tanh f -f c δ ( 30
)
where f is the frequency and δ is a constant that specifies the frequency range of the transition from 1 to 0 (see Fig. 6(b)).

The envelop function ψ keeps the phase of the delayed reflection coefficient unchanged. The target function to be modeled is then: the method finds the biggest point-to-point oriented distance between the real parts (Eq. ( 28)) of R(ω) and R BC ω). Fig. 7(a) shows the model at the end of the first iteration. The initial solution of R BC (ω) targets the first peak in the real part where ℜ(R(ω)) = 1. In Fig. 7(b), the targeted peak is where ℜ(R(ω)) = -1. At an iteration i, initial guesses for p i and µ i are found to target a single peak at the time. Fig. 7 illustrates how the model propagates from low to high frequencies. The algorithm is considered converged after 20 iterations as the maximum point-to-point error on the modulus is below 1%, where we define the error as

R(ω) = ψ( f ) × R r | x=1 m (31)
E = max (|R(ω) -R BC (ω, n 0 = 20)|).
The D-TDIBC model used here is, thus, made of 20 poles and residues (p k , µ k ), that is to say only 40 complex constants (cf. Appendix A).

At ω = 0, the modulus of the reflection coefficient to be modeled is |R(ω = 0)| = 1. Such a condition cannot be fulfilled by the D-TDIBC as each of the PBF has a null modulus at the zero-frequency limit: |R k (ω = 0)| = 0. However, the frequency The reflection coefficient modeled by the poles and residues fulfills the causality and reality conditions by construction of the method. However, it is interesting to note that for a small number of poles and residues (hence for a low accuracy of the fit) the passivity condition is violated, as shown in Fig. 7(a) where |R(ω)| > 1 around 0.15 kHz. Nonetheless, for a higher number of poles and residues, the passivity condition is also verified.

One-dimensional wave propagation

In this section, the behavior of the D-TDIBC is investigated first on a basic one-dimensional wave propagation problem. Then the time domain response of each PBF is inspected to highlight the mechanism of D-TDIBC.

The one-dimensional wave propagation is computed using the AVBP solver. AVBP is a three-dimensional fully compressible Navier-Stokes equation solver. A two-step Taylor-Galerkin scheme, called TTGC is used. TTGC is third-order accurate in space and time [START_REF] Quartapelle | High-order Taylor-Galerkin methods for nonlinear multidimensional problems[END_REF][START_REF] Olivier Colin | Development of high-order Taylor-Galerkin schemes for LES[END_REF]. Characteristics boundary conditions (NSCBC) are used [START_REF] Poinsot | Boundary conditions for direct simulations of compressible viscous flows[END_REF][START_REF] Lodato | Three-dimensional boundary conditions for direct and large-eddy simulation of compressible viscous flow[END_REF][START_REF] Nicoud | Defining wave amplitude in characteristic boundary conditions[END_REF][START_REF] Selle | The actual impedance of non-reflecting boundary conditions: implications for the computation of resonators[END_REF]. In the NSCBC framework the characteristic waves are evaluated at the boundaries. TDIBC prescribes the ingoing characteristic wave from the outgoing characteristic wave and is, thus, consistent with the NSCBC formalism.

As the reflection coefficient modeled in x r is R r = -1, the wave is a priori expected to be fully reflected: the amplitude of the reflected wave must be equal to the one of the incident wave. As R r is real-valued, the reflected wave is expected to be centered in x = 0 m after a time t = T (see Sec. 4.1).

Fig. 8 illustrates the propagation of the pressure wave from time t = 0 up to t = T . At t = 0 the pressure field corresponds to the initial solution shown in Fig. 6(a). At t = 0.25T the wave has propagated in the positive x direction and is crossing the boundary condition at x = x BC = 1 m. At t = 0.5T the pressure wave has completely left the computational domain. This result is expected: at time t = 0.5T the wave is propagating inside the modeled domain m , i.e. the truncated portion (gray area in Fig. 6(a)). If one were to compute the complete domain, at t = 0.5T the acoustic pressure p ′ (x) would be zero for x ∈ [-1, 1] m. At t = 0.75T the pressure wave is re-injected in the computational domain by the boundary condition. The reflected pressure wave has the same amplitude than the incident wave (1 Pa) as expected. At t = T the reflected Gaussian wave is centered at x = 0 m which means that the time delay is precisely prescribed by D-TDIBC.

The acoustic pressure and velocity at the boundary are recorded at the boundary condition at x BC = 1 m during the simulation. The results are shown in Fig. 9(a). Additionally, the acoustic energy E a contained in the domain is shown. For the sake of clarity, each variable is normalized by its maximum. At first no acoustic activity is seen as the wave propagates inside the domain, the acoustic energy is maximum. The acoustic pressure p ′ at the boundary in Fig. 9(a) is consistent with the results in Fig. 8. At first no acoustic pressure is seen by the boundary until the Gaussian wave crosses the boundary at t ≃ 0.25T ≃ 2.5 ms. After the wave has crossed the boundary, the acoustic pressure returns to zero. After the time delay τ (Fig. 9(a)) at t ≃ 0.75T ≃ 7.5 ms the reflected wave is re-entering the computational domain with the same amplitude as the incident wave. Finally, the acoustic pressure p ′ returns to zero. The amplitude of the reflected acoustic velocity u ′ (t ≃ 0.75T ≃ 7.5 ms) is also the same as the incident wave (t ≃ 0.25T ≃ 2.5 ms) although the sign is changed: u ′ is negative. This indicates that the wave is traveling in the negative x direction after the reflection at x r = 1.75 m. This result is consistent with the imposition of a reflection coefficient R r = -1 in x r . When the wave crosses the boundary the acoustic energy initially contained in the domain is lost due to the acoustic flux at the boundary condition in x BC . The acoustic energy contained in the domain retrieves its initial value when the reflected wave re-enters the domain. This stresses the fact D-TDIBC conserves the acoustic energy while imposing a time delay.

In their work, Jaensch et al. [START_REF] Jaensch | Hybrid CFD/low-order modeling of nonlinear thermoacoustic oscillations[END_REF] solve a similar problem: they impose a pure time delay at a boundary condition with a one-dimensional incident Gaussian wave. As discussed in Sec. 1, their method uses a first order upwind numerical scheme to discretize the one-dimensional LEE over the truncated domain. Their method is, thus, inherently dissipative and the dissipation increases linearly with the cell size of the spatial discretization. In order to achieve a dissipation level comparable to the one obtained here, it was necessary to use 1000 points. As a consequence, in the state-space model used, the state matrix is of dimension 2000 × 2000, i.e. 4 million scalar values. D-TDIBC seems to yield comparable results for a lower memory storage requirement as only 40 complex constants were used to model the delayed reflection coefficient.

Fig. 9(b) illustrates the contribution of each PBF R k . The ingoing wave A in n (t) imposed by D-TDIBC is the sum of all the R k . Before t ≃ 7 ms the ingoing wave is null: the PBFs are canceling each other out. At t ≃ 7 ms the individual contributions of the PBFs stop canceling themselves out: they add up to create the re-entering characteristic wave. This stresses how the modeling procedure of R(ω) is critical for an accurate prediction of acoustic delays.

Validation for a combustion chamber

In this section, we investigate the thermoacoustic instability [START_REF] Rayleigh | On the theory of resonance[END_REF] in a laminar experimental setup called INTRIG Burner (IMFT, Toulouse). First, a DNS of the full setup (called "FULL") is conducted. A second DNS (called "TRUNCATED") in which the computational domain is truncated after the flame is conducted. In the TRUNCATED case, D-TDIBC is used to model the acoustic properties of the truncated domain as in Sec. 4. The results obtained using D-TDIBC will be compared with those of the FULL configuration.

Experimental setup

The INTRIG Burner is used to study a lean premixed laminar methane-air flame attached on a cylinder. The operating point corresponds to an equivalence ratio φ = 0.75 and a bulk velocity of the fresh gases of u b = 0.8 m s -1 . The associated laminar flame speed and adiabatic temperature are s l = 0.23 m s -1 and T ad = 1920 K.

The experimental rig is shown in Fig. 10. The gaseous methane-air premixture is injected upstream of the glass ball array located at x = -0.367 m. The flow is then laminarized by the glass balls and the honey comb panels. A lean premixed methane-air laminar flame then attaches to a cylindrical stainless-steel flame holder of with a diameter d = 8 mm. The combustion chamber has a constant cross section of h = 34 mm by l = 94 mm.

Numerical setup

As the flow is two-dimensional, 2D Direct Numerical Simulation (DNS) approach can be used [START_REF] Miguel-Brebion | Joint experimental and numerical study of the influence of flame holder temperature on the stabilization of a laminar methane flame on a cylinder[END_REF][START_REF] Xavier | Experimental and numerical investigation of flames stabilised behind rotating cylinders: interaction of flames with a moving wall[END_REF][START_REF] Brebion | Joint Numerical and Experimental Study of Thermoacoustic Instabilities[END_REF]. A 19-species mechanism (called LU19) is used [START_REF] Lu | A criterion based on computational singular perturbation for the identification of quasi steady state species: a reduced mechanism for methane oxidation with no chemistry[END_REF] for the modeling of chemical kinetics. Schmidt and Prandtl numbers are assumed constant: Pr = 0.6 and Sc = 0.6. The full mesh is composed of 625 000 cells. The reduced mesh is composed of 545 000 cells, i.e. a reduction of 13%. The flame thickness is δ 0 L = 680 µm and at least 11 points are used to resolve the flame front with a mesh resolution at the flame of x = 60 µm [START_REF] Miguel-Brebion | Joint experimental and numerical study of the influence of flame holder temperature on the stabilization of a laminar methane flame on a cylinder[END_REF][START_REF] Xavier | Experimental and numerical investigation of flames stabilised behind rotating cylinders: interaction of flames with a moving wall[END_REF][START_REF] Brebion | Joint Numerical and Experimental Study of Thermoacoustic Instabilities[END_REF]. A typical flame, here at a stable operating point, is shown in Fig. 10 (bottom).

The inlet acoustic boundary (glass balls at x = -0.367 m) is a hard wall (i.e. u(ω) = 0 m s -1 at all frequencies). It was checked in the experiment that this is a good approximation [START_REF] Brebion | Joint Numerical and Experimental Study of Thermoacoustic Instabilities[END_REF]. This assumption is close to realistic conditions as the inlet BC corresponds to the glass balls array. The outlet acoustic boundary (at x = 0.35 m) is an open-end (i.e. p(ω) = 0 Pa at all frequencies). This issue should be taken care of if one were willing to compute accurately the thermoacoustic stability of the configuration. However, the goal of this paper is more modest: it aims to validate D-TDIBC. The simple acoustic conditions used here are thus acceptable as the reference case is a numerical simulation. Two simulations are carried out and investigated in this section (Fig. 11): the "FULL" and the "TRUNCATED" cases. The outlet boundary condition of the FULL domain is located as x = 0.35 m and the outlet D-TDIBC boundary condition of the TRUNCATED domain is located at x = 0.1 m. Fig. 11 shows the computational domain and the boundary conditions used in the three simulations. The grey area in the TRUNCATED case corresponds to the domain where the acoustic wave propagation is modeled by a time delay imposed by D-TDIBC. Table 1 summarizes the parameters used in the simulation. The truncation of the gray domain in Fig. 11 is possible as the flow is strongly one-dimensional and homogeneous. For instance, the axial velocity is 4 orders of magnitude higher than the maximum of the vertical velocity in the truncated portion.

D-TDIBC: reflection coefficient model

The Iterative Multi-Pole Modeling Technique presented in Sec. 3 is used to model the acoustic wave propagation in the truncated part of the INTRIG Burner. The outlet boundary conditions of the Full and Truncated domains are located at x max = 0.35 m and x BC = 0.10 m, respectively. Consequently, the length of the domain to be modeled by D-TDIBC is L = 0.25 m. The value of x BC must be large enough to ensure that combustion is complete and only acoustics take place between x BC and x max . Moreover, the sound speed must be homogeneous in this zone. This condition has been verified a posteriori by analyzing the results of the FULL configuration simulation: the sound speed fluctuations are less than 0.5% of the mean value c 0 = 830 m s. According to duct acoustics theory, plane waves propagates in ducts at frequencies lower that a cutoff frequency f c [START_REF] Munjal | Theory of a two source-location method for direct experimental evaluation of the four pole parameters of an aeroacoustic element[END_REF]. This frequency is given by the speed of sound c 0 and the lowest height in the cross section.

In the INTRIG setup, the plane waves propagate at frequencies lower that f c = 3.5 kHz. The filtering procedure used in (Sec. 4.2) is applied here so that the frequencies lower than f c = 3.5 kHz are accurately modeled. The results are shown in Fig. 12. The reflection coefficient model used here consists in a set of 18 (p k , µ k ) and verifies causality, reality and passivity conditions (cf. Appendix B). 13)) for a time-delay of τ = 0.30 ms in the laminar combustor of Fig. 10.

Results and discussion

The flow setup in Fig. 10 is thermoacoustically unstable: the flame oscillates and couples to the acoustic modes of the setup [START_REF] Rayleigh | On the theory of resonance[END_REF]. The instability is sustained in DNS of the FULL and TRUNCATED configurations (Fig. 11). Fig. 13(a) illustrates the Sound Pressure Level (SPL) spectra expressed in dB and defined as SPL = 20 log(p/p ref ) based on a reference pressure level of p ref = 2 • 10 -5 Pa. The FULL simulation, i.e. computing the full domain, shows a strong acoustic activity with many amplified acoustic modes. Four principle modes are found, with frequencies 203 Hz, 534 Hz, 782 Hz and 1150 Hz, respectively. The mode shapes can be obtained by (1) performing a spectral analysis at several axial locations of the burner, (2) extracting the amplitude of the Fourier transform at the peak frequency in Fig. 13(a) and (3) normalizing by the maximum pressure amplitude. The three first modes observed in the FULL simulation are a quarter wave mode (Fig. 13(b)), a three-quarter wave mode (Fig. 13(c)) and a five-quarter wave mode (Fig. 13(d)).

The TRUNCATED simulation also exhibits thermoacoustic instabilities and is in good agreement with the FULL simulation. Fig. 13(a) illustrates the ability of D-TDIBC to accurately predict the broadband acoustic activity, despite the clear nonlinear excitation of higher harmonics, due to the high pressure amplitudes levels. In particular, the agreement between values of the excited frequencies (Fig. 13(a)) and the mode shapes (Figs. 13(b),(c) and (d)) is remarkable.

Some discrepancies are visible in the magnitude of the SPL spectrum, where the TRUNCATED simulations overall underpredict the intensity of the SPLs. This discrepancy can be attributed to two concurrent issues. First, the reflection coefficient used at the outlet of the FULL simulation is |R FULL | = 1 for all frequencies, while this is true only in a selected frequency range of interest (due to the filtering procedure used, cf. Fig. 12 and Fig. 6(b)) for the TRUNCATED simulation. Outside of that range, the magnitude of the reflection coefficient drops below 1, hence causing spurious acoustic energy loss at the D-TDIBC boundary. Second, the large amplitude of the pressure oscillations at the limit cycle may trigger losses due to nonlinear acoustic effects which are directly simulated in the FULL simulations in a portion of the domain that is not present in the TRUNCATED simulations.

Overall, Fig. 13 demonstrates that D-TDIBC can satisfactorily recover the acoustic eigenmodes of a domain while simulating only a part of it. It also gives the opportunity to investigate the thermoacoustic stability of different setups while simulating the same computational domain: the length of the truncated domain can be adjusted within the same DNS simply by introducing a different time delay at the boundary x = x BC (Fig. 10), without the need to create a new computational mesh.

Conclusion

A novel numerical strategy called Delayed-Time Domain Impedance Boundary Condition (D-TDIBC) has been derived from an existing TDIBC approach [START_REF] Fung | Broadband time-domain impedance models[END_REF][START_REF] Fung | Time-domain impedance boundary conditions for computational acoustics and aeroacoustics[END_REF][START_REF] Scalo | Compressible turbulent channel flow with impedance boundary conditions[END_REF][START_REF] Lin | High-fidelity simulation of a standing-wave thermoacoustic-piezoelectric engine[END_REF] in order to account for time delays in Navier-Stokes simulation. The method allows to mimic acoustic wave propagation in domains larger than the computational domain used in the numerical study. The benefits are twofold: (1) computational cost can be saved and spatial resolution can focus on areas of interest, (2) any additional length can be added to the computational domain allowing straightforward investigations of geometry changes downstream or upstream of the zone of interest. While the former mainly concerns structured meshes, the latter applies to all types of grids. D-TDIBC is an extension of Time Domain Impedance Boundary Condition (TDIBC) methods [START_REF] Fung | Broadband time-domain impedance models[END_REF][START_REF] Fung | Time-domain impedance boundary conditions for computational acoustics and aeroacoustics[END_REF][START_REF] Scalo | Compressible turbulent channel flow with impedance boundary conditions[END_REF][START_REF] Lin | High-fidelity simulation of a standing-wave thermoacoustic-piezoelectric engine[END_REF]. The proposed methodology works with the reflection coefficient R(ω) and gives accurate models for time delayed fully reflecting acoustic boundary conditions such as closed and open-end conditions. These models were used in the time domain simulations and were found to give accurate acoustic wave propagation properties for both one-dimensional and two-dimensional reactive Navier-Stokes simulations.

The one-dimensional wave propagation results illustrate the capacity of D-TDIBC to accurately impose a time delay while conserving the acoustic energy in case of a full reflection. The investigation on individual Pole Base Functions (PBFs) stresses the need for an accurate modeling strategy. Indeed, the time delayed response of the boundary results in a cancellation of individual responses: the PBFs. This is accomplished through the diversity of their amplitudes, pseudo-periods and phases.

The two-dimensional simulation focuses on the thermoacoustic stability of a laminar combustion chamber (the INTRIG Burner [START_REF] Miguel-Brebion | Joint experimental and numerical study of the influence of flame holder temperature on the stabilization of a laminar methane flame on a cylinder[END_REF][START_REF] Xavier | Experimental and numerical investigation of flames stabilised behind rotating cylinders: interaction of flames with a moving wall[END_REF][START_REF] Brebion | Joint Numerical and Experimental Study of Thermoacoustic Instabilities[END_REF]) with a lean premixed laminar methane-air flame attached on a cylinder. A portion of the outlet domain was truncated and the acoustic wave propagation was modeled with D-TDIBC. The D-TDIBC simulation was compared to a reference simulation where the full domain is taken into account. D-TDIBC was found to give accurate predictions of the thermoacoustic stability while saving computational resources. The frequencies, amplitudes and mode shapes are in excellent agreement with the reference simulation.

Fig. 1 .

 1 Fig. 1. Orientation of ingoing and outgoing acoustic waves A in n and A out n with respect to an Impedance Boundary (IB) with outward pointing normal n.

Fig. 3 .

 3 Fig. 3. Acoustic variables for two canonical limit cases: (a) the flanged open-end, and (b) the closed-end (hard wall) cases. Both configurations are fully reflective and have a purely real reflection coefficient R at all frequencies in x = x 2 . The truncation of the grey areas can be modeled by the imposition of a time delay τ . The resulting delayed reflection coefficient R τ at x = x 1 is complex. At the bottom of the figure both real (solid lines) and imaginary parts (dashed lines) of R τ and R are shown for both cases.

Fig. 3

 3 Fig.3illustrates the impact of an arbitrary time delay τ on two classical limit cases: the flanged open-end, i.e. p(ω) = 0 (Fig.3(a)), and the closed-end, i.e. u(ω) = 0 (Fig.3(b)). These cases have purely real reflection coefficients R(ω) at the acoustic boundary located in x = x 2 . For all frequencies ω, the flanged open end corresponds to a reflection coefficient of R(ω) = 1 and the closed end to R(ω) = -1. The delayed reflection coefficient R τ (ω) becomes complex-valued and frequency-dependent even though R(ω) is purely real and frequency-independent.In this work we build upon the multi-oscillator TDIBC methodology of Fung and Ju[START_REF] Fung | Time-domain impedance boundary conditions for computational acoustics and aeroacoustics[END_REF], with the phase-parameter correction by Lin et al.[START_REF] Lin | High-fidelity simulation of a standing-wave thermoacoustic-piezoelectric engine[END_REF], to create a novel TDIBC formulation, called Delayed-Time Domain Impedance Boundary Condition (D-TDIBC), relying on the determination of the partial fraction representation (via poles and residues) of the delayed complex reflection coefficient, as outlined in Sec. 3.1. It is important to note that by adopting Fung and Ju's formalism, which focuses on the causal reconstruction of the softness (or absorption) coefficient, it is indeed appropriate to talk about a multi-oscillator fitting strategy, where the fitted softness is the linear superposition of the softnesses of various causal oscillators, which can be each representative of real acoustically absorptive panels; on the other hand, the decomposition of the reflection coefficient in partial fractions as done in the present work (see Eq. (13)) is not tied to a particular physical interpretation. Due to the peculiar mathematical nature of the delayed complex reflection coefficient (Eq. (12), plotted in Fig.3), exhibiting periodic variations in the frequency domain, the proposed D-TDIBC strategy was developed to handle a large number of poles and residues. While the focus herein is on the reflection coefficient R, the same fitting methodology outlined in Sec. 3 can be extended to the wall softness S.

Fig. 4 .

 4 Fig. 4. Real part of a typical PBF R k (Eq. (21)). Three properties are visualized: ω 0,k , h 0,k and ω ǫ . At the resonant angular frequency ω 0,k , the real part of R k is maximum. The peak height h 0,k is the value of R k (ω) at ω = ω 0,k . For a given percentage ǫ (ǫ ∈ [0, 1]), the half-width ω ǫ is defined on the right-hand-side. ω -and ω + are the solutions of the equation ℜ (R k ) = ǫ and are used to define the width.

Fig. 5 .

 5 Fig. 5. (a) Influence of the c k parameters for f 0,k = 1 2π ω 0,k = 100 Hz and h 0,k = 1. (b) Real part of a purely delayed reflection coefficient R τ for a time delay τ . The real part is a periodic function (cosine) of period T = 2π τ .

  (a)). The mean speed of sound is c 0 = 350 m s -1 . This perturbation will propagate on the positive x direction. The physical domain is split into two sub-domains: = c + m . The computational domain c is x ∈ [-1, 1] m and the modeled domain m is x ∈ [1, 1.75] m (gray area in Fig.6(a)). The D-TDIBC is used to model the wave propagation in m . At the left boundary condition (x = x l = -1 m) the reflection coefficient is R l = 1 and at the right boundary condition (x = x r = 1.75 m) the reflection coefficient is R r = -1. In the computational domain, the reflection coefficient at x = 1 m accounting for R r has to be modeled by a time delay (see Sec. 1) leading to R r | x=1 m = -e -iωτ . The time delay corresponding to this modeled domain is τ = 2(x rx BC )/c 0 = 4.28 ms. The acoustic time for the initial wave located in x 0

Fig. 7

 7 Fig.7shows the model at iterations 1, 2, 5 and 20 of the Iterative Multi-Pole Modeling Technique (Algorithm 1). R(ω) is the filtered delayed reflection coefficient shown in Fig.6(b) and R BC ω) is the reflection coefficient model at a given iteration defined in Eq.[START_REF] Zhong | A controllable canonical form implementation of time domain impedance boundary conditions for broadband aeroacoustic computation[END_REF]. Their real and imaginary parts are shown in the top and bottom plots, respectively. At each iteration,

Fig. 6 .

 6 Fig. 6. (a) Initial solution and visualization of the computational ( c ) and modeled ( m ) domains. (b) Theoretical (gray solid line) and filtered (black solid line) delayed reflection coefficients. The filter shape (Eq. (30)) is shown (black dashed line).

Fig. 7 .

 7 Fig. 7. Reflection coefficient R BC (ω) modeled by the Iterative Multi-Pole Modeling Technique algorithm (Algorithm 1 in Sec. 3.3) a several iterations. An accurate model is obtained for n 0 = 20, where n 0 is the number of PBF R k as in Eq. (20).

Fig. 8 .

 8 Fig. 8. Propagation and reflection of a Gaussian pressure wave in the domain c + m by simulating only c . The acoustic properties of m are imposed using D-TDIBC.

Fig. 9 .

 9 Fig. 9. (a) Acoustic velocity (dashed line) and pressure (solid line with circles) at the boundary condition. (b) Temporal response of 20 PBF (dashed and dotted lines) and of the ingoing wave imposed by the D-TDIBC (black solid line).

Fig. 10 .

 10 Fig. 10. Sketch of a transverse cut of the INTRIG Burner (top). The computational domain of the FULL configuration goes from the glass balls array (x = -0.367 m) to the exhaust (x = 0.35 m) while it is cut (at x = x BC = 0.1 m) in the TRUNCATED configuration. A zoom (bottom) of the flame region of the initial solution used in DNS is shown. The light grey and dark grey regions represent the zero and high CO 2 mass fraction levels. The streamlines are plotted in white.

Fig. 11 .

 11 Fig. 11. Schematic of the computational domains and BC used for the two simulations: FULL and TRUNCATED cases.

1 R

 1 TRUNCATED (ω) = e -iωτD-TDIBC

Fig. 12 .

 12 Fig. 12. Model of the delayed reflection coefficient R BC (ω) with n 0 = 18 (see Eq. (13)) for a time-delay of τ = 0.30 ms in the laminar combustor of Fig.10.

Fig. 13 .

 13 Fig. 13. (a) Spectra obtained for the FULL and the TRUNCATED simulations. D-TDIBC allows an accurate prediction of the thermoacoustic stability observed in the FULL configuration. The spectra are based on a probe located at x = -8.5 mm. The shapes of the first (b), second (c) and third (d) modes correspond to the acoustic eigenmodes of a left-hand side closed and right-hand side open duct modified by the jump of characteristic impedance between the fresh and burnt gases. The markers correspond to the probes in the simulation (lines are added for the sake of clarity).

  Poles and residues used in the two-dimensional reactive DNS of the INTRIG Burner.

					3.48e+03	4.75e+03	7.56e+02
	k = 2 k = 3 k = 4 k = 5 k = 6 k = 7 k = 8 k = 9 k = 10 k = 11 k = 12 k = 13 k = 14 k = 15 k = 16 k = 17 k = 18 k = 19 k = 20	-1.10e+04 -8.21e+04 -1.33e+05 -6.17e+04 -3.22e+04 -1.73e+04 -1.19e+02 -2.11e+00 -9.15e+03 -1.32e+05 6.55e+02 -5.04e+03 -2.25e+03 -3.80e+01 -7.67e+03 6.89e+00 -1.50e+01 3.82e+03 -6.06e+00	-7.38e+03 -3.92e+04 -4.66e+04 -1.39e+04 -5.12e+03 -2.04e+03 2.65e+07 6.77e+04 -8.40e+02 -1.15e+05 5.49e+02 -4.09e+02 -1.70e+02 -8.77e+07 -9.34e+02 1.17e-01 1.75e+06 3.80e+02 -7.54e+05	-9.52e+02 -1.38e+03 -1.50e+03 -1.30e+03 -1.15e+03 -1.03e+03 -2.17e+02 -1.89e+00 -9.34e+02 -7.95e+03 -5.20e+02 -9.40e+02 -9.78e+02 -7.84e+01 -1.80e+03 -1.61e+02 -2.93e+01 -1.50e+03 -9.40e+00	1.42e+03 2.89e+03 4.29e+03 5.79e+03 7.26e+03 8.73e+03 -9.77e-04 -5.88e-05 1.02e+04 9.11e+03 6.20e+02 1.16e+04 1.30e+04 3.40e-05 1.48e+04 9.46e+03 -2.51e-04 1.50e+04 7.56e-05	1.71e+03 3.20e+03 4.54e+03 5.93e+03 7.35e+03 8.79e+03 2.17e+02 1.89e+00 1.02e+04 1.21e+04 8.09e+02 1.16e+04 1.30e+04 7.84e+01 1.49e+04 9.46e+03 2.93e+01 1.51e+04 9.40e+00	2.72e+02 5.09e+02 7.23e+02 9.44e+02 1.17e+03 1.40e+03 3.46e+01 3.01e-01 1.63e+03 1.92e+03 1.29e+02 1.85e+03 2.07e+03 1.25e+01 2.37e+03 1.51e+03 4.67e+00 2.40e+03 1.50e+00

Appendix B. D-TDIBC model: 2D reactive DNS of the INTRIG Burner (Table B.3)

Table B.3

This is true when the e +iωt convention is adopted.

Acknowledgements

The research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP/2007-2013)/ERC Grant Agreement ERC-AdG 319067-INTECOCIS. This work was granted access to the high-performance computing resources of CINES under the allocation A0012B07036 made by Grand Equipement National de Calcul Intensif. The support of Calmip for access to the computational resources of EOS is acknowledged.

Carlo Scalo acknowledges the support of Purdue's Rosen Center for Advanced Computing (RCAC), the Air Force Office of Scientific Research (AFOSR) grants FA9550-16-1-0209 (Core Program) and FA9550-18-1-0292 (2018 Young Investigator Program), and NSF's Fluids Dynamics Program (Award # 1706474).

Appendix A. D-TDIBC model: 1D wave propagation problem (Table A.2) Table A.2

Poles and residues used in the one-dimensional pure delayed Gaussian wave propagation problem.


F Deschaux-Beaume

F Schmidt

N Fréty

T Cutard

J C Boyer

C Levaillant

NUMERICAL SIMULATION OF FAILURE PREDICTION FOR CERAMIC TOOLS: COMPARISON WITH FORGING TEST RESULTS

Keywords: finite element method, fracture prediction, ceramics, forging

A fracture prediction criterion for brittle materials has been introduced in the POLLUX finite element code, in order to predict the risk-of-rupture of ceramic tools during a forging operation. POLLUX is a software dedicated to the simulation of forging operations, initially developed by INSA (Lyon). The chosen probabilistic fracture model is based on the weakest link theory and on the statistical theory of Weibull. A surface approach or a volume approach can be retained, from the type of critical flaws in the ceramic. Two different criteria are available in order to characterise the stress state considering the tension normal stresses and neglecting the compression stresses. An identification procedure of critical flaw type is presented for a particular ceramic material. Statistical parameters of ceramic fracture have been experimentally determined using bending tests performed in environmental conditions close to the forging ones. A constitutive equation of the workpiece material has been proposed, issued from torsion tests. In order to validate the model in the case of ceramic tools subjected to multiaxial stress states, a particular configuration has been defined to compare the simulation predictions with the experimental results. A forging test has been then developed, in which a billet of superalloy is formed in a ceramic tool up to its fracture at the temperature of 1423K. The experimental distribution of tool fracture, according to the strain of the billet, are in good agreement with fracture predictions computed by the simulation.

INTRODUCTION

Isothermal forging is a very effective forming process for some refractory alloys like nickelbased materials. This technique should allow the production of near net shape parts, directly from simple shape preforms. However, the process is based on the superplastic properties of the forming materials, and it must be carried out at high temperature and at low strain rate. This means that the material for the forging dies, which will be used at the same temperature as the workpiece, must conserve high mechanical properties at high temperature. The forging tools are actually made in Mo-based alloys, but the high reactivity of this element with oxygen imposes to carry out the process under an inert gas. The use of ceramic dies should be a less expensive way, as the good oxidation resistance of these materials does not require a controlled atmosphere. However, ceramics are inherently brittle. This characteristic leads in particular to a wide variation in the material strength. Usual design techniques are therefore not suited to ceramic components, and statistic and reliability analysis are needed. Statistical fracture models, based on the weakest link theory, have been previously developed. The first probabilistic approach used to account for the scatter in fracture strength of brittle materials was introduced by Weibull [1]. This model has been initially proposed for constant uniaxial stress fields. New criteria were later added to extend the Weibull approach to multiaxial stress fields [START_REF]Fracture[END_REF][START_REF] Evans | [END_REF]. An other model, developed by Batdorf, combines the weakest link theory and the linear elastic fracture [4]. Except in some simple cases, it is necessary to use numerical simulation in order to predict the risk of failure of ceramic parts with these models. Some studies have been previously conducted, in which brittle fracture models were coupled to finite element simulations in order to calculate the failure probability of particular ceramic structures. In some cases, a post processor has been developed for a specific software, to analyse the probabilistic behaviour of simple ceramic structures [START_REF] Berdin | Proc. of Fracture Processes in Brittle Materials[END_REF] or of uniaxially loaded ceramic composites [START_REF] Guillaumat | Proc. of Journée nationales sur les composites[END_REF]. However, most of the studies in this field have been made by the National Aeronautics and Space Administration [START_REF] Gyekenyesi | [END_REF][8][9]. A computerised design program for reliability evaluation of ceramic components has been developed from this research work. This software can be coupled to the ANSYS, ABAQUS, NASTRAN and COSMOS commercial finite element programs. Nevertheless, any program previously quoted is not convenient for the simulation of the forging process, as large strains take place in the workpiece.

The main objective of this study is to develop a finite element model for the failure prediction of ceramic forging dies as a function of the service conditions in order to improve the design of these components.

DESCRIPTION OF THE FAILURE PREDICTION SOFTWARE 2.1. Brittle fracture modelling

The wide variation of the strength measured for brittle materials is due to the nature and distribution of intrinsic microscopic flaws. These flaws are unavoidably present as a result of the material manufacturing and machining operations, or of a damage process. The distinction must be made between volume and surface flaws, which lead to different strength distributions. In the case of a coarse material microstructure, the volume flaws are generally more critical than the surface ones. However, the surface flaws may be more critical for example in the case of rough machining or of oxidised surface states. The identification of the nature of critical flaws allows one to choose the analysis approach. If a surface one is chosen, the volume flaws are assumed not to be responsible for the fracture. As the flaw distribution cannot be determined, some statistical functions are used to describe the strength variability of a component. The Weibull model has been chosen in the present study because the statistical parameters are experimentally easy to determine. The Weibull theory associates a strength  c to each elementary volume, Ve (or elementary surface Se for the surface approach). The strength of the material is characterised by the distribution function of  c , denoted F(), which is defined as the number of all elementary volumes (surfaces) having  c ≤, divided by the total volume (external surface). This function also gives the probability P of randomly choosing an elementary volume (surface) having a strength equal to or less than , or also the probability that this volume fails under the stress :

P( c ≤)=F() (1)
Any distribution function may be written in the form:

F()=1-exp(-()) (2)
where () must be a positive and non decreasing function, vanishing at a value  u , called the threshold stress.

The advantage of this formula appears when considering the failure probability of the structure constituted by all the elementary volumes (surfaces), using the weakest link theory. When supposing that the total volume V (surface S) fails if any one of its parts fails, then the probability of non failure P NF , under an uniform uniaxial load , is equal to:

P NF =(1-P) n =exp(-n()) (3)
where n is the number of elementary volumes (surfaces) of the structure (V=nVe or S=nSe)

The failure probability P n of the structure takes then the simple form:

P n =1-exp(-n()) (4)
Weibull has chosen a power form for the function (), which is the most simple function satisfying the general conditions previously described:

( )     u  0       m ( 5
)
where  u ,  0 and m are parameters depending on the material characteristics.

The threshold stress  u is usually taken equal to zero. The two other parameters can then be easily identified by comparison with results of fracture tests.

In the case of a non uniform uniaxial load , the failure probability takes the integral form:

P v  1 exp  1 Ve   0       m dV V        (6.1)
for the volume approach, or

P s  1 exp  1 Se   0       m dS S        (6.2)
for the surface approach.

The previous relationship can be generalised to multiaxial stress states. Two strategies are available: the first consists in considering an equivalent stress which characterises the loading at each point. This criteria leads sometimes to unsafe estimates of the failure probability; the second consists in assuming the independence of the actions of each principal tensile stress [START_REF]Fracture[END_REF].

Taking into account the first approach, Weibull proposed to calculate an equivalent stress by averaging the normal tensile stresses in all directions. This model leads to a complex integral formulation, which may increase the computation time.

In order to obtain shorter computation times, only the maximal normal stress has been considered in our model, as other investigators did [START_REF] Berdin | Proc. of Fracture Processes in Brittle Materials[END_REF]. This variable is also the maximal principal stress at each point of the structure.

In the second approach, the fracture probability depends on all the tensile principal stresses. The independence of actions of the principal stresses means that the reliability of the structure is the product of the reliability of the individual principal stresses acting independently. The fracture probability is then given by the following relationship: for the surface approach, (

P v  1 exp  1 Ve  1  0       m   2  0       m  
where  1 ,  2 and  3 are the positive principal stresses or take the zero value for the compressive principal stresses.

It is interesting to note that in the both models, the compressive stresses are not taken into account. The fracture is supposed to be always the result of tensile stresses. This choice has been made because the compressive fracture strength of ceramics is 3 or 4 times higher than the tensile strength. It is so very unlikely that structure fails under the effect of compressive stresses.

Numerical introduction of the model

The model previously described has been introduced as a subroutine in the POLLUX finite element code [10]. POLLUX is a specialised software developed by INSA-Lyon [11] in order to predict the stress state in the workpiece and the surrounding tools during a forging process. First, the thermomechanical analysis is achieved, leading to the calculation of the stress tensor at each point of the tools, which are supposed to show a perfectly elastic behaviour. Then, the failure probabilities of the tools are determined, using the Weibull parameters and the criteria chosen by the user.

If Ve or Se are arbitrary taken equal to 1, the discrete decomposition of the models gives: 2) for the surface approach, if the principle of the independent actions (PIA) is chosen.

P v  1 exp   1  0       m V ig   ig  iel  
 1 ,  2 ,  3 are the principal stresses in the Gauss point ig of element iel or in the contour node inc( 1 ≥ 2 ≥ 3 ) or are equal to zero for the compressive principal stress. V(ig) is the volume associated to the Gauss point ig, and S(inc) is the surface associated to the contour node inc.

The elementary failure probability DPr, or local risk-of -rupture, is also calculated at each Gauss point in the case of the volume approach, and interpolated to the nodes. This variable corresponds to the fracture probability of an elementary volume of constant size equal to 1, subjected to the stress state of the considered point. The values can be stored in a data file for graphical display of the critical regions. An other data file, containing the failure probability of the overall structure as a function of the time or tool displacement, is also generated. At the end of each increment, the program compares failure probability of the tools to the maximal permissible probability, chosen by the user. The computation is interrupted when this limit or the maximal tool course is achieved.

RELIABILITY ANALYSIS PROCEDURE FOR CERAMIC FORGING TOOLS

Identification of critical flaw type

Reliability analysis of ceramic structures with the software previously described, first requires to choose a failure criteria corresponding to the material and the use conditions. If the choice between MPS and PIA criteria can not be made by considering physical observations, the choice of the approach type (surface or volume) must be based on these observations. In the case of the forging process with ceramic tools, the material strength is modified, because of damage inherent to the use conditions. Two types of flaws can be created: surface flaws, essentially produced by the oxidation of the material at the forging temperature, and internal flaws (microcracks), which can be formed because of the mechanical loading of the tool or because of the thermal cycles. In order to identify the critical flaw type formed during the forging process, we have compared the fracture strength of the initial material and of the material subjected to a damaging factor (oxidation, mechanical or thermal cycles). The failure probability of the material, as a function of the stress, has been determined using four point bending tests (see section 3.2), at room temperature. Fig. 1 shows the experimental points and the corresponding interpolated Weibull functions for the different states of the material. First, it is to be noticed that the oxidised material has a weaker strength than initial material. Then it appears that mechanical and thermal cycles have no significant effects on the material characteristics. The strength distribution of the material after mechanical cycles is indeed very close to the initial material one, and the strength distribution of the material after oxidation and thermal cycles is approximately the same as for the uncycled oxidised material. This means that the surface flaws, created by the oxidation of the material, are more critical than the initial flaws, or the internal flaws resulting of the cyclic loading. According to this result, the surface approach must be chosen for the calculation of the ceramic tools reliability during the forging process.

Determination of statistic parameters

The second step of the reliability analysis of the ceramic tool consists in the experimental determination of the statistical parameters of the Weibull model corresponding to the material and to the processing conditions. Four points bending tests (Fig. 2) have then been carried out at 1423K on a set of N specimens previously oxidised at 1423K during 50 hours. The failure probability of each sample is given by the formula:

P ex  i N  1 ( 10
)
where i is the rank of the maximal stress  max in the corresponding specimen at the fracture point.  max is given by: where S eff is the effective surface, determined from the surface in tension, S, of a specimen:

 max  3 2 F L  l   bh 2 (
S eff    max       m dS S  (13)
Eq. ( 12) can also be written: The Weibull plot, which expresses the evolution of ln(ln(1/1-P)) as a function of ln( max ), is usually used to represent the experimental results.

The Weibull parameters are estimated by fitting the equation ( 14) in the experimental points of a Weibull plot, using the least squares method (Fig. 3). The Weibull parameters m and  u obtained are respectively equal to 12.4 and 580 MPa.

C e g ra

Determination of the workpiece rheology

The further set of data required for the reliability analysis is the rheological behaviour of the workpiece material. In the forging conditions, the superalloys have a superplastic behaviour. In consequence, the flow stress depends on the strain rate.

In order to determine the constitutive equation of this material, torsion tests have been carried out at the forging temperature, with 3 different strain rates ranging from 0.03 to 0.3 s -1 . During the test, the experimental torque  is measured as a function of the turn number T. The equivalent strain rate ú and strain , calculated at the periphery of the specimen, are given by:

ú   2 ú T R 3b with ú T  dT dt (15)  = 2 TR 3b (16)
where R and b are respectively the radius and the length of the specimen. The Von Mises equivalent stress  is given by:

   3 2R 3  ˜ n  ˜ m   (17)
where ˜ n and ˜ m represent the sensitivity of the torque to the strain and to the strain rate.

The evolution of  versus  are plotted in Fig. 4 for the different strain rates. The curves show a stress peak ( p), followed by a softening and a stabilisation of the stress values ( s). According to Berdin et al. [12], the stress peak observed in this type of material at the low strain levels is due to a dynamic recrystallisation phenomenon. The material behaviour can be approximated by the model:

  1    p  s (18) with p  Ap Ý  mp , s  As Ý  ms ( 19
)
where the different constants have been determined by fitting the model to the experimental results. mp and ms are the strain rate sensitivity coefficients, respectively equal to 0.22 and 0.15. Ap and As depend only on the temperature and are respectively equal to 150.1 and 109.6.  represents the recrystallized fraction, given by an Avrami type equation [12]. Figure 4 shows a good agreement between the experimental results and the curves calculated from this model.

CONFRONTATION WITH FORGING-FRACTURE TEST RESULTS

Description of the forging test

In order to validate the model and the approach previously described, fracture tests have been carried out on ceramic tools in near forging conditions. The results were compared to the numerical results of the tool reliability. Fig. 5 shows the samples chosen for this confrontation. The ceramic and the workpiece were heated until 1423K in an induction furnace. During the test under isothermal conditions, a superalloy workpiece is pressed by a punch in the ceramic specimen until its fracture. The forging load and the punch displacement are measured all along the tests. After the tests, realised in the same conditions on a set of ceramic specimens, the failure probability calculated from Eq. ( 10) is expressed as a function of the punch displacement. A program of 18 fracture tests has been planned to obtain an experimental distribution of ceramic fracture.

Failure probability calculation

At the moment, only the two first tests have been carried out, but encouraging results have been obtained. A good agreement between the experimental results and the numerical simulation predictions is indeed observed. First, the calculation of the stress state in the specimen has shown that the orthoradial stress is the maximal component of the stress tensor in the ceramic (Fig. 6). This predicts that at the fracture, the cracks in the ceramic specimen must be longitudinal. This has been verified experimentally and can clearly be observed on Figure 5. Both MPS and PIA criteria available for the calculation have been tested, in order to determine the more efficient for the present case. Nevertheless, the choice seems to be difficult, because of the weak difference obtained from the both criteria. This is due to the fact that for the ceramic geometry chosen, one of the principal stresses is always predominant, and the two others have no significant effects on the risk-of-rupture. The more appropriated criteria will probably be easier to determine when considering the totality of the 18 experimental tests.

CONCLUSION

A fracture prediction program for brittle materials has been introduced in the finite element code POLLUX, in order to predict the risk-of-rupture of ceramic tools during high temperature forging operations. For an efficient use, the following characteristics must be known: -the type of the critical flaws in the ceramic; -the statistical parameters of the fracture strength distribution of the ceramic; -the metal workpiece rheology. These characteristics can be measured by classical mechanical tests (bending, torsion or tension tests). The software is then a very useful prediction tool for the fracture of ceramic parts: -it permits to predict the failure probability of a ceramic tool with a given geometry, as a function of its use conditions.

-it can also give elements to predict the localisation and the orientation of the fracture cracks. Such a software can so be very helpful for the design of the geometry and of the dimension of large ceramic tools. The comparison between the numerical simulation results and the forging fracture tests has shown a good agreement. Complementary experimental tests will be performed in a mean future in order to have a more complete comparison , particularly in terms of the statistical distribution of the fracture loads and displacements. These results are expected to be very helpful to choose the more appropriated fracture criteria (PIA or MPS) for the ceramic material used in the present study. Other failure criteria could also be considered, for example taking into account the shear stresses in the ceramic. Future developments of the software forecast the introduction of this type of criteria.

  for the surface approach, if the maximal principal stress criteria (MPS) is chosen, and:P v  1 exp   i

1 Fig. 1 .

 11 Fig.1. Fracture strength distribution of the initial and of the damaged material.

Fig. 2 .

 2 Fig.2. Four point bending testsThe Weibull distribution corresponding to the bending tests, in the case of surface critical flaws, can be expressed as follow:

4 N o m b r e d 1 Fig. 3 .

 413 Fig.3. Weibull's plot of the four point bending strength at 1423K of ceramic specimens previously oxidised at 1423K during 50 hours.

1 Fig. 4 .

 14 Fig.4. Stress-strain curves obtained from torsion test at 1423K and calculated from the model.

Fig. 5 .

 5 Fig.5. Ceramic specimen and workpiece used for the fracture tests

1 Fig. 6 . 1 Fig. 7 .

 1617 Fig.6. Orthoradial stress distribution obtained from the forging test simulation after a punch displacement of 4 mm.