Mickaël Rivier

Jean Lachaud

Pietro Marco Congedo

Ablative thermal protection system under uncertainties including pyrolysis gas composition

Keywords: Sensitivity analysis, Uncertainty quantification, Porous media, Heat-shield ablation Analyse de sensibilité, Quantification des incertitudes, Milieu poreux, Abalation de bouclier thermique

Spacecrafts such as Stardust (NASA, 2006) are protected by an ablative Thermal Protection System (TPS) for their hypersonic atmospheric entry. A new generation of TPS material, called Phenolic Impregnated Carbon Ablator (PICA), has been introduced with the Stardust mission. This new generation of low density carbon-phenolic composites is now widely used in the aerospace industry. Complex heat and mass transfer phenomena coupled to phenolic pyrolysis and pyrolysis gas chemistry occur in the material during atmospheric entry. Computer programs, as the Porous material Analysis Toolbox based on OpenFoam (PATO) released open source by NASA, allow to study the material response. In this study, a non-intrusive Anchored Analysis of Variance (Anchored-ANOVA) method has been interfaced with PATO to perform low-cost sensitivity analysis on this problem featuring a large number of uncertain parameters. Then, a Polynomial-Chaos method has been employed in order to compute the statistics of some quantities of interest for the atmospheric entry of the Stardust capsule, by taking into account uncertainties on effective material properties and pyrolysis gas composition. This first study including pyrolysis gas composition uncertainties shows their key contribution to the variability of the quantities of interest.

Introduction

Space exploration missions often include entering a planetary atmosphere at hypersonic speed. A high enthalpy hypersonic shock forms around the spacecraft and kinetic energy is progressively dissipated into heat [START_REF] Anderson | Hypersonic and high temperature gas dynamics[END_REF]. Heat is transferred to the surface of the spacecraft by radiation and convection. A suitable heat shield is needed to protect the payload. The level of heat flux increases with entry speed and atmospheric density. For fast hypersonic entries, typically faster than 8 km/s from earth orbit, ablative materials are used as Thermal Protection Systems (TPS). These materials mitigate the incoming heat through phase changes, chemical reactions, and material removal [START_REF] Duffa | Ablative Thermal Protection Systems Modeling[END_REF]. A low-density porous carbon/phenolic composite called PICA was used for the Stardust comet-dust sample-return capsule, which reentered the Earth's atmosphere at 12.7 km/s [START_REF] Stackpoole | Post-flight evaluation of stardust sample return capsule forebody heatshield material[END_REF]. PICA is made of a carbon fiber preform partially impregnated with phenolic resin. Chemistry mechanisms (simplified illustra;on) Scanning Electron Microscopy (SEM) of a carbon preform before phenolic impregna;on SEM of PICA (virgin) [START_REF] Stackpoole | Post-flight evaluation of stardust sample return capsule forebody heatshield material[END_REF]

Boundary layer

Non--viscous flow

Porous microstructure

Flow/surface interac;ons --Abla;on (oxida;on, sublima;on, spalla;on), --Heat and mass surface balance (1) Subsurface phenomena (not always) --Boundary layer gases entering in the material (3) --In--depth oxida;on [START_REF] Stackpoole | Post-flight evaluation of stardust sample return capsule forebody heatshield material[END_REF] Heat and mass transport --Solid density conserva;on (1) --Steady--state gas conserva;on (1) --Gas density conserva;on (2) --Gas momentum conserva;on (2) --Gas element/species conserva;on (3) --Energy conserva;on [START_REF] Anderson | Hypersonic and high temperature gas dynamics[END_REF] Produc;on phenomena --Pyrolysis produc;on rate (1) --Pyrolysis species/elements produced (3) --Chemistry of the gases [START_REF] Stackpoole | Post-flight evaluation of stardust sample return capsule forebody heatshield material[END_REF] Long distance effects --Radia;on coking rad. Heat ≈ 6000 K 1 : in all material--response models (type 1) 2 : in some material--response models (type 2) 3 : in detailed material--response models (type 3) During atmospheric entry, carbon/phenolic materials undergo thermal degradation and ultimately recession captured by the following physico-chemical phenomena (Figure 1). The phenolic polymer thermally decomposes and progressively carbonizes into a low density carbon form, losing mass while releasing pyrolysis gases. The pyrolysis gases percolate and diffuse to the surface through the network of pores. Reactions within the pyrolysis-gas mixture (homogeneous reactions) and between pyrolysis gases and the char take place with possible coking effects (heterogeneous reactions). Mixing and reaction of the pyrolysis gases with boundary layer gases into the pores of the material occur when boundary layer gases penetrate in the material by forced convection or due to fast diffusion at low pressures [START_REF] Lachaud | Detailed chemical equilibrium model for porous ablative materials[END_REF]. At the surface, the material is removed by ablation and the outer surface recedes. Depending on entry conditions, ablation may be caused by heterogeneous chemical reactions (oxidation, nitridation), phase change (sublimation), and possibly mechanical erosion (often called spallation).

Material response (PATO)

A detailed heat and mass transfer model is required to estimate the performance of the porous material and design the thermal protection system. Two important design criteria are the expected level of recession and the maximum back wall temperature. The key parameter uncertainties are propagated to obtain the design uncertainties to be used in the margin policy [START_REF] Wright | Sizing and margins assessment of the mars science laboratory aeroshell thermal protection system[END_REF]. A Monte Carlo approach has been developed and used to propagate uncertainties on material properties and areoheating conditions for the design of the NASA Mars Science Laboratory [START_REF] Sepka | Monte carlo approach to fiat uncertainties with applications for mars science laboratory[END_REF] and Orion [START_REF] Sepka | A reliability comparison of classical and stochastic thickness margin approaches for the orion heat shield resulting from material property uncertainties[END_REF] spacecrafts. In these study however, no uncertainty is attributed to the pyrolysis gas composition. Recent publications have shown that the pyrolysis gas composition strongly varies depending on temperature and heating rate [START_REF] Wong | Quantitative determination of species production from the pyrolysis of a phenol-formaldehyde resin[END_REF][START_REF] Bessire | Pyrolysis of phenolic impregnated carbon ablator (pica)[END_REF][START_REF] Wong | Detailed analysis of species production from the pyrolysis of phenolic impregnated carbon ablator[END_REF]. The pyrolysis gases are composed of carbon, oxygen and hydrogen elements. The pyrolysis gas composition influences the pyrolysis gas enthalpy -which impacts heat transfer in the porous material -and the boundary layer chemistry -which controls the ablation rates and the surface temperature. For this first analysis we will allow an uncertainty of 10% on these elements. The composition in term of species is then computed in each cell of the mesh and at each time step using an equilibrium chemistry solver [START_REF] Lachaud | Detailed chemical equilibrium model for porous ablative materials[END_REF]. This makes the computation very costly and requires the use of the low-cost uncertainty quantification methods. In the literature, low-cost uncertainty propagation has been already performed alongside Global Sensitivity Analysis for problems of natural convection in [START_REF] Fajraoui | Analyzing natural convection in porous enclosure with polynomial chaos expansions: Effect of thermal dispersion, anisotropic permeability and heterogeneity[END_REF]. Uncertainty analyses have also been performed on surface ablation rates and their effect on aeroheating predictions for Mars entry in [START_REF] Wang | Uncertainty analysis of laminar and turbulent aeroheating predictions for mars entry[END_REF], and on ablation problems in plasma wind tunnel [START_REF] Turchi | Thermochemical ablation modeling forward uncertainty analysis-part i: Numerical methods and effect of model parameters[END_REF][START_REF] Turchi | Thermochemical ablation modeling forward uncertainty analysis-part ii: Application to plasma wind-tunnel testing[END_REF].

In section 2, we present the problem studied and the physical hypotheses. In section 3, we present the inverse analysis method implemented in the study. The results of the uncertainty quantification analysis are presented in section 4. Finally, section 5 draws some conclusions and perspectives.

Definition of the uncertainty analysis problem and hypotheses

For this first analysis, we chose to study the entry of Stardust, that was the first mission using a low-density carbon-phenolic ablator in 2006. The thermal response of the TPS has been be studied at the stagnation point during the whole reentry, from entry interface to cool down. As in the state-of-the-art design approach we assumed that the problem is locally mono-dimensional. The actual thickness of the ablative material was two inches [START_REF] Stackpoole | Post-flight evaluation of stardust sample return capsule forebody heatshield material[END_REF], therefore we used this value. Adiabatic conditions are used at the bondline. A convective boundary condition is used at the surface of the ablative material. Surface total pressure and heat flux were taken from reference [START_REF] Trumble | Post-flight aerothermal analysis of the stardust sample return capsule[END_REF]. Figure 2 illustrates the temperature evolution during the atmospheric entry at the stagnation point and in-depth under the stagnation point with nominal TACOT properties.

The analysis is performed using the properties of the Theoretical Ablative Composite for Open Testing (TACOT). Its composition and properties are comparable to PICA. Nominal TACOT properties are available in the open literature [START_REF] Lachaud | Ablation test-case series #2[END_REF]. Volume-wise, TACOT is made of 10% of carbon fibers, 10% of phenolic resin, and is 80% porous.

Inria

Model

A generic heat and mass transfer model for porous media has been recently developed and documented [START_REF] Lachaud | A generic local thermal equilibrium model for porous reactive materials submitted to high temperatures[END_REF]. It is suitable to model ablative heat shields. For the sake of conciseness, we only present a short summary in this section. The model was developed for porous materials containing several solid phases and a single gas phase. The detailed chemical interactions occurring between the solid phases and the gas phase are modeled at the pore scale assuming local thermal equilibrium. Homogenized models were obtained for solid pyrolysis, pyrolysis species injection in the gas phase, heterogeneous reactions between the solid phases and the gas phase, and homogeneous reactions in the gas phase. The chemistry models were integrated in a macroscopic model making use of volume-averaged governing equations for the conservation of solid mass, gas mass, species (finite-rate chemistry) or elements (equilibrium chemistry), momentum, and energy. The model is implemented in the Porous material Analysis Toolbox (PATO), distributed Open Source by NASA. First-order implicit finite-volume schemes in time and space [START_REF] Jasak | Error analysis and estimation for the finite volume method with applications to fluid flows[END_REF][START_REF] Weller | A tensorial approach to computational continuum mechanics using object-oriented techniques[END_REF], which have been shown to provide excellent convergence and accuracy [START_REF] Lachaud | Porous material analysis toolbox based on openfoam and applications[END_REF][START_REF] Lachaud | Detailed chemical equilibrium model for porous ablative materials[END_REF], were used for the simulations presented in Section 4. In this study, we used an equilibrium chemistry model that is equivalent to the reference NASA TPS design model [START_REF] Omidy | Code-to-code comparison, and material response modeling of Stardust and MSL using PATO and FIAT[END_REF]. The current approach is to assume that the elemental pyrolysis gas composition is fixed. To save on computational time, precomputed tables are used to obtain the gas composition (species) and properties (enthalpy, viscosity, molar mass). In the current study, we wish to vary the elemental pyrolysis gas composition. The pyrolysis gas composition in term of species is therefore computed from pyrolysis gas elemental composition in each cell of the mesh and at each time step using an equilibrium chemistry solver, as described in [START_REF] Lachaud | Detailed chemical equilibrium model for porous ablative materials[END_REF]. The surface ablation rate is computed using the thin film coefficient approach, also known as B' approach, accounting for the change in pyrolysis gas composition injected in the boundary layer [START_REF] Lachaud | Detailed chemical equilibrium model for porous ablative materials[END_REF].

Uncertain parameters and associated uncertainties

In previous studies, uncertain material property parameters have been identified for PICA [START_REF] Sepka | Monte carlo approach to fiat uncertainties with applications for mars science laboratory[END_REF][START_REF] Wright | Sizing and margins assessment of the mars science laboratory aeroshell thermal protection system[END_REF][START_REF] Mahzari | An inverse parameter estimation methodology for the analysis of aeroheating and thermal protection system experimental data[END_REF]. We decided to include the same uncertain parameters in our study. We also added a set of new parameters to assess the effect of the pyrolysis gas composition on the material response as described in the introduction. The nominal elemental composition of the pyrolysis gases for TACOT are, in mole fractions, C (0.206), H (0.679), O (0.115). In total, we have used twenty-seven uncertain parameters in the TACOT material model. We have attributed 5 to 10 % uncertainty to each of them as follows (the number in brackets is the label used to identify each uncertainty in the following of this paper):

• Density (1) and volume fraction (2) of the fibrous preform (5% uncertainty),

• Density (3) and volume fraction (4) of the phenolic matrix (5% uncertainty),

• Virgin's ( 5) and char's ( 6) permeability (5% uncertainty),

• Pyrolysis model (10% uncertainty):

-Elementary composition of the pyrolysis gases in Carbon [START_REF] Sepka | A reliability comparison of classical and stochastic thickness margin approaches for the orion heat shield resulting from material property uncertainties[END_REF] , Hydrogen [START_REF] Wong | Quantitative determination of species production from the pyrolysis of a phenol-formaldehyde resin[END_REF] and Oxygen ( 9), -Pyrolysis reaction 1: pre-exponential factor [START_REF] Wong | Detailed analysis of species production from the pyrolysis of phenolic impregnated carbon ablator[END_REF], activation energy [START_REF] Fajraoui | Analyzing natural convection in porous enclosure with polynomial chaos expansions: Effect of thermal dispersion, anisotropic permeability and heterogeneity[END_REF], pyrolysis enthalpy [START_REF] Wang | Uncertainty analysis of laminar and turbulent aeroheating predictions for mars entry[END_REF], -Pyrolysis reaction 2: pre-exponential factor (13), activation energy [START_REF] Turchi | Thermochemical ablation modeling forward uncertainty analysis-part ii: Application to plasma wind-tunnel testing[END_REF], pyrolysis enthalpy (15),

• Thermal properties of virgin material (5% uncertainty): heat capacity ( 16), orthogonal conductivity [START_REF] Lachaud | A generic local thermal equilibrium model for porous reactive materials submitted to high temperatures[END_REF], radial conductivities [START_REF] Jasak | Error analysis and estimation for the finite volume method with applications to fluid flows[END_REF][START_REF] Weller | A tensorial approach to computational continuum mechanics using object-oriented techniques[END_REF], emissivity [START_REF] Lachaud | Porous material analysis toolbox based on openfoam and applications[END_REF], reflectivity (21),

• Thermal properties of charred material (5% uncertainty): heat capacity [START_REF] Mahzari | An inverse parameter estimation methodology for the analysis of aeroheating and thermal protection system experimental data[END_REF], orthogonal conductivity [START_REF] Archer | Sensitivity measures, anova-like techniques and the use of bootstrap[END_REF], radial conductivities [START_REF] Congedo | Shape optimization of an airfoil in a bzt flow with multiple-source uncertainties[END_REF][START_REF] Tang | Sensitivity analysis using anchored anova expansion and high-order moments computation[END_REF], emissivity [START_REF] Yang | Adaptive {ANOVA} decomposition of stochastic incompressible and compressible flows[END_REF], reflectivity [START_REF] Xiu | The Wiener-Askey Polynomial Chaos for Stochastic Differential Equations[END_REF].

We chose to simplify the constraint of elementary fractions summing to one through maintaining the relative ratios and normalizing the elementary composition. Practically, given a random draw of the mass fractions y i within the ±10% interval, normalized mass fractions

y i = y i k y k
will actually be given to PATO. Moreover, with 1 and 2 the uncertain volume fraction of the fibrous preform and the virgin matrix respectively, the virgin and charred porosities ε v and ε c are computed as follows:

ε v = 1 -1 -2 ε c = 1 -1 - 2 2
Finally, one may note the presence of 3D conductivities in the list of uncertain parameters, which clashes with the mono-dimensional assumption made earlier. Such parameters are left in Inria the study to artificially increase the input dimension and verify the capability of the proposed approach to detect their null impact and discard them.

Sensitivity and uncertainty analysis theory and tools

Let us consider a stochastic differential equation of the form:

L (x, ξ, φ) = f (x, ξ) (1)
where L is a non-linear spatial differential operator (for instance, the steady Navier-Stokes operator) depending on a set of uncertainties, designated with the random vector ξ (of dimension the number of uncertain parameters in the problem) and where f (x, ξ) is a source term depending on x and ξ. In the following, we drop the dependence on x in order to simplify the notation.

The solution of the stochastic equation ( 1) is φ(ξ), which is a function of the space variable

x ∈ R d and of ξ ∈ Ξ = Ξ 1 × • • • × Ξ N (Ξ ⊂ R N ) and ξ ∈ Ξ -→ φ(ξ) ∈ L 2 (Ξ, p(ξ)),
where

p(ξ) = N i=1 p(ξ i )
is the probability density function of ξ. One of the objective of Uncertainty Quantification is to compute the statistics of the quantity of interest, i.e. φ(ξ).

We can define the central statistical moment of φ of order n as

µ n (f ) = Ξ (φ(ξ) -E(φ)) n p(ξ)dξ, (2)
where E(φ) indicates the expected value of φ

E(φ) = Ξ φ(ξ)p(ξ)dξ. ( 3
)
In the following, we indicate with σ 2 , the variance (second-order moment). We illustrate the main concepts of the ANOVA-decomposition in Section 3.1. Then, to clearly present the context of uncertainty analysis theory and provide a comprehensive understanding of the approach followed in this work, we will use as illustration a mono-dimensional heat transfer problem presented in Section 3.2. The UQ methods are then described in Section 3.3.

ANOVA-based decomposition

Let us suppose that the response of a given system of interest can be represented by a N -dimensional function:

y = φ(ξ) = φ(ξ 1 , ξ 2 , • • • , ξ N ) (4)
We consider Eq. ( 4) in its functional expansion form as follows

y = φ 0 + N 1 i N φ i (ξ i ) + N 1 i<j N φ ij (ξ i , ξ j ) + • • • + φ 1,2,••• ,N (ξ 1 , ξ 2 , • • • , ξ N )
or in compact form using a multi index system:

y = φ s0 + 2 N -1 j=1 φ sj (ξ sj ) (5)
The multi indices s j are defined such as

s 0 = (0, 0, 0, • • • , 0) s 1 = (1, 0, 0, • • • , 0) s 2 = (0, 1, 0, • • • , 0) . . . s N = (0, 0, 0, • • • , 1) s N +1 = (1, 1, 0, • • • , 0) s N +2 = (1, 0, 1, • • • , 0) . . . s N = (1, 1, 1, • • • , 1) (6)
where N = 2 N -1. The representation of Eq. ( 5) is called ANOVA (Analysis Of Variance) decomposition [START_REF] Archer | Sensitivity measures, anova-like techniques and the use of bootstrap[END_REF] of

φ(ξ), if for any j ∈ {1, • • • , N }, R φ sj (ξ sj )p(ξ i ) dξ i = 0 for ξ i ∈ {ξ sj } (7)
It follows from Eq. ( 7) the orthogonality of ANOVA component terms, namely

E(φ sj φ s k ) = 0 for j = k (8)
ANOVA allows identifying the contribution of a given stochastic parameter to the total variance of an output quantity. Meanwhile, we obviously have

E(φ sj ) = 0 for j = 1, • • • , N
Note that the terms in the ANOVA decomposition can be expressed as integrals of φ(ξ). Indeed, we have

E(Y ) = φ 0 E(Y |ξ i ) = φ 0 + φ i (ξ i ) E(Y |ξ i , ξ j ) = φ 0 + φ i (ξ i ) + φ j (ξ j ) + φ ij (ξ i , ξ j ) (9)
and so on, where E(Y |•) denotes the conditional expectation.

Analytical solution for transient heat transfer

Deterministic problem

Effective heat transfer is the main mode of energy transport in most porous materials. Let us consider a homogeneous semi-infinite unidimensional medium. Under the assumption of constant material properties, the transient heat transfer equation is given by

∂ t T = α ∂ 2 x T ( 10
)
where α = k/(ρ × c p ) is the diffusivity. We will consider a medium initially at the room temperature T (x, t = 0) = T 0 = 300 K. Its surface temperature is held at T (x = 0, t) = T w = 1646 K during the experiment. Laplace transform is used for its resolution as presented in Appendix A.

The temperature profile as a function of time and space is given by

T (x, t) = T 0 + (T w -T 0 ) erfc x 2 √ αt ( 11
)
Inria
where erfc is the complementary error function. The temperature profiles computed for times of 1, 10 and 60 seconds are plotted in Fig. 3 for a representative medium of diffusivity of 10 -7 m 2 /s. We will only consider in the illustrations that follow the first centimeter of the medium.

We see here that the hypothesis that the medium is semi-infinite does not play a role on the result, as the heat wave hasn't reached the one centimeter mark after one minute of heating.

Formulation under uncertainty

Let us now formulate the problem presented in Eq. 10 under an uncertainty quantification perspective. In particular, let assume that two parameters are affected by some variability and/or are not well-known: a 5% of variation in terms of min/max is then imposed on T w and α considering a uniform distribution (with respect to the deterministic values previously used, denoted in the following as T wm and α m , respectively). The problem is now formulated as follows:

∂ t T = α ∂ 2 x T = α m × (0.95 + (1.05 -0.95) * ξ 2 ) ∂ 2 x T ( 12
)
where T w = T wm × (0.95 + (1.05 -0.95) * ξ 1 ), and

ξ 1 , ξ 2 vary in [0, 1].
Several methods can be used in order to solve the problem defined in Eq. 12. In this work, we use systematically the so-called non-intrusive methods: this means that a single deterministic computation (used to solve for example the differential operator defined in Eq. 10) is replaced with a whole set of such computations, each one of those being run for specific values of the uncertain conditions. The final solution can be then written as follows:

T (x, t, ξ 1 , ξ 2 ) = T 0 + (T wm (0.95 + 0.1ξ 1 ) -T 0 ) erfc x 2 α m (0.95 + 0.1ξ 2 )t . ( 13
)
Now, let us show how the computation of the variance and the computation of the contribution of each source of uncertainty can be reduced only to the computation of some integrals on the analytical solution shown in Eq. 13 for some fixed values of x and t:

The ANOVA functional expansion (more details are provided in the next subsection) is a unique tool for assessing the contribution of each uncertainty (and of the interactions) to the global variance. This is computed as follows (variables x and t are dropped since this does not change the following developments)

T (ξ) = T 0 + T ξ1 + T ξ2 + T ξ1ξ2 , ( 14
)
where

T 0 = Ξ 2
T (ξ)p(ξ)dξ;

T ξ1 = Ξ T (ξ)p(ξ 2 )dξ 2 -T 0 ; T ξ2 = Ξ T (ξ)p(ξ 1 )dξ 1 -T 0 ; T ξ1ξ2 = T (ξ) -T ξ1 -T ξ2 -T 0 . ( 15
)
The overall variance σ 2 can be computed by means of the ANOVA expansion as

σ 2 = σ 2 ξ1 + σ 2 ξ2 + σ 2 ξ1ξ2 , (16)
where

σ 2 ξ1 = Ξ T 2 ξ1 p(ξ 1 )dξ 1 ; σ 2 ξ2 = Ξ T 2 ξ2 p(ξ 2 )dξ 2 ; σ 2 ξ1ξ2 = Ξ 2 T 2 ξ1ξ2 p(ξ)dξ. ( 17
)
Note that σ 2 ξ1 , σ 2 ξ2 represent the unique contribution of ξ 1 and ξ 2 to the global variance σ 2 , respectively. Moreover, σ 2 ξ1ξ2 represents the contribution given by the interaction between ξ 1 and ξ 2 .

Note that only integrals of the expression defined in Eq. 13 are required, in order to compute the contributions to the variance for fixed values of (x, t). In Figure 3 (on the right), the solution is then represented in terms of mean and the associated error bars (square root of the variance, i.e. standard deviation). Figure 4 illustrates the variance of the temperature t = 60 s, induced by each uncertainty. Note that the contribution T w is predominant and explains most of the global variability of the temperature. This simple example illustrates the interest in propagating some physical input uncertainties through numerical models.

Non-intrusive formulations for expensive computer codes

Unfortunately, generally, it is not possible to compute an analytical solution of the problem defined in Eq. 1. This could require the resolution of a complex system of equation, relying on a numerical approximation of the solution on some discretized grid of the numerical domain. Note then that computing the integrals of Eq. 17 can be very costly. Moreover, some additional issues Inria could arise in the presence of a large number of uncertainties or if the quantity of interest features some discontinuities. The real challenge is then to formulate an efficient numerical algorithm permitting to build an accurate representation of the quantity of interest as a function of input uncertainties.

As previously mentioned, only non-intrusive strategy are targeted in this work. In particular, here, we tackle a problem featuring a large number of uncertainties, that can be very challenging to solve, due to the so-called Curse of Dimensionality. It refers to the loss of convergence and the infeasible number of calculations needed when the number of parameters increases, for any chosen method. We have partially cured this problem with a two-steps approach. First, we applied an anchored-ANOVA approach on the complete problem. This analysis permits to compute the hierarchy and detect the most important uncertainties. Note that this approach only needs a very reduced number of deterministic simulations to perform uncertainty propagation and sensitivity analysis. In a second step, we applied a Polynomial-Chaos approach for treating the subspace including only the predominant parameters, in order to provide a good representation of the quantity of interest in the reduced stochastic space.

As mentioned before, due to the non-intrusivity of the stochastic methods considered here, the coupling with PATO, or any other heat and mass transfer computational model, is very straightforward: it reduces to the creation of a small interface for building automatically PATO input parameters files for each set of uncertain conditions.

Both methods are described briefly in the following. For more details, refer to [START_REF] Congedo | Shape optimization of an airfoil in a bzt flow with multiple-source uncertainties[END_REF] for the Polynomial-based method and to [START_REF] Tang | Sensitivity analysis using anchored anova expansion and high-order moments computation[END_REF] for the anchored-ANOVA approach.

Anchored-ANOVA approach: Definitions and basic notions

In order to introduce the less expensive anchored ANOVA, the Dirac measure is used for the integrals of Eq. ( 9):

p(ξ i ) dξ i = δ(ξ i -c i ) dξ i for i = 1, • • • , N (18)
Thus, p(ξ) dξ = δ(ξ -c) dξ. The point c = (c 1 , • • • , c N ) is called "anchor point". Hence, the ANOVA component terms in Eq. ( 9) can be expressed as follows:

φ(c) = φ 0 φ(c|ξ i ) = φ 0 + φ i (ξ i ) φ(c|ξ i , ξ j ) = φ 0 + φ i (ξ i ) + φ j (ξ j ) + φ ij (ξ i , ξ j ) . . . ( 19
)
The formulae in Eq. ( 19) are used to quantify the expectation and variance of the component functions, by simply evaluating the model outputs at chosen sampling points. For more details, see [START_REF] Tang | Sensitivity analysis using anchored anova expansion and high-order moments computation[END_REF]. This permits a strong reduction of the computational cost, since this avoids the computation of several integrals. Moreover, a variance-based adaptive criterion (see for more details [START_REF] Yang | Adaptive {ANOVA} decomposition of stochastic incompressible and compressible flows[END_REF]) is used in order to compute the so-called effective dimension and to evaluate highorder interactions with a reduced computational cost. The order at which the ANOVA model is truncated, is called effective dimension, beyond which the difference between the ANOVA model and the truncated expansion in a certain measure is very small. This implies that we will ignore terms in the ANOVA model corresponding to interactions exceeding the fixed threshold.

In this work, a covariance decomposition of the output variance has been considered, as proposed in [START_REF] Tang | Sensitivity analysis using anchored anova expansion and high-order moments computation[END_REF], in order to accurately compute the statistics using the anchored-ANOVA expansion. The covariance decomposition makes the result less sensitive to the choice of the anchor point if a full expansion of the anchored ANOVA is employed.

Polynomial-chaos based approach

Under specific conditions, a stochastic process can be expressed as a spectral expansion based on suitable orthogonal polynomials, with weights associated to a particular probability density function. The first study in this field is the Wiener (1938) process. The basic idea is to project the variables of the problem onto a stochastic space spanned by a complete set of orthogonal polynomials Ψ that are functions of random variables ξ. For example, the unknown variable φ has the following spectral representation:

φ (ξ) = ∞ i=0 φ i Ψ i (ξ)) . ( 20
)
In practice, the series in Eq. ( 20) has to be truncated in terms of the polynomial degree p 0 , where the total number of terms of the series M is determined by:

M + 1 = (N + p 0 )! N ! p 0 ! , ( 21
)
where N is the dimensionality of the uncertainty vector ξ. Each polynomial Ψ i (ξ) is a multivariate polynomial form which involves tensorization of 1D polynomial forms. The polynomial basis is chosen accordingly to the Wiener-Askey scheme [START_REF] Xiu | The Wiener-Askey Polynomial Chaos for Stochastic Differential Equations[END_REF] in order to select orthogonal polynomials Inria with respect to the probability density function p(ξ) of the input. The orthogonality can be advantageously used to compute the coefficients of the expansion in a non-intrusive PC framework

φ i = φ(ξ), Ψ i (ξ) Ψ i (ξ), Ψ i (ξ) , ∀i. (22)
Several approaches can be used to estimate PC coefficients. The approach used in this study is based on quadrature formula. As a consequence, the solution of a deterministic problem for each quadrature point is required.

For further details, see Congedo et al. [START_REF] Congedo | Shape optimization of an airfoil in a bzt flow with multiple-source uncertainties[END_REF]. In both cases, once the chaos polynomials and the associated φ i coefficients are computed, the expected value and the variance of the stochastic solution φ i (ξ) are obtained from:

E P C = φ 0 (23)
V ar P C = N i=1 φ 2 i Ψ 2 i ( 24
)
Another interesting property of PC expansion is to make easier sensitivity analysis based on the analysis of variance decomposition (ANOVA). It can be easily computed by using some interesting properties of the previous development [START_REF] Crestaux | Polynomial chaos expansion for sensitivity analysis[END_REF]. Let us recall here that the contribution to the variance of a given random variables with index k, i.e. the first order Sobol's index, can be obtained by:

S k = i∈α φ 2 i Ψ 2 i (ξ) V ar P C ( 25
)
where α represent the set of indexes associated to a given uncertainty k. For more details, Ref.

[28] is strongly recommended.

Results

The Anchored-ANOVA method is applied to the problem presented in section 2. We present first the convergence analysis with respect to the number of samples considered to get well-converged statistics of all the quantities of interest. In Anchored-ANOVA, the first-order analysis is based on a chosen number of points per direction. We have reported in Table 1 the outcome of this analysis for the temperature computed at a depth of 1.5 cm at a time equal to 80s, in terms of decreasing contributions to the variance, for eight and sixteen points along each direction, respectively (which makes a total number of runs of the solver of 216 and 432, respectively). As it can be observed, errors are quite small, indicating that quantities are converged with only eight points per direction. Same conclusions can be drawn for all the other quantities of interest considered in this work, i.e. temperatures at different depths, the virgin front and the char front.

Results shown in the following rely then on this analysis.

Let us now analyse the results from a quantitative and qualitative point of view. Figure 5 presents the transformation of the ablative material with error bars, which represent the spreading of the quantities of interest in terms of standard deviation. The plotted results are the surface recession -due to ablation -and the propagation of the pyrolysis front. The char 2% and virgin 98% are used in the ablation community to identify almost completely charred material (2% left of virgin matrix) and almost pristine virgin material (98% left of virgin matrix). The charring zone is considered to be between these two accepted limits [START_REF] Lachaud | Ablation test-case series #2[END_REF]. Note that, since the heating is quite soft at the beginning of the simulation, the material is fully pyrolysed before being ablated. This explains why the char 2% and the wall curves are clearly separated in this case. The observed variability is very small for the three curves. In fact, for the virgin 98%, which is the worst case, the Coefficient of Variation (standard deviation to mean ratio) is of the order of 3%.

A physical analysis of the results is now presented in terms of the contribution of each uncertainty to several quantities, namely the recession, the virgin front, the char front and the surface temperature, as a function of time. The results are presented in Figure 6. The uncertainties propagated on the material properties are in agreements with previous studies [START_REF] Sepka | Monte carlo approach to fiat uncertainties with applications for mars science laboratory[END_REF][START_REF] Wright | Sizing and margins assessment of the mars science laboratory aeroshell thermal protection system[END_REF][START_REF] Mahzari | An inverse parameter estimation methodology for the analysis of aeroheating and thermal protection system experimental data[END_REF]. Concerning the recession and the location of the Char 2%, the maximal standard deviation is of the order of 0.03-0.04 cm. As it can be observed in Figure 6a and b, the location of the recession and char 2% mostly depend on both the parameters influencing the heat transfer in the charred material and the pyrolysis parameters. The virgin 98% is predominantly influenced by the uncertainty on the Activation energy 1, which controls the initiation of pyrolysis. Finally, the surface temperature is mostly driven by the char's emissivity except at the start of the pyrolysis reaction, where the activation energies and gas composition are predominant.

We would like to point out a new result. The uncertainty on the elemental pyrolysis gas composition (Carbon, Hydrogen and Oxygen) clearly induces variability on the studied quantities. The effect of the these uncertainties is observed on the in-depth temperature evolution as well, as shown in Figure 7, where the contributions of each uncertainty to the variance of the temperature is computed over the time at different depths in the material. The standard deviation of the temperature takes the highest value at a depth of 0.7 cm (100 K), while it remains quite small for the other depths considered here.

With the advancement of the ablation front, contributions of the virgin's parameters are quickly overtaken by char's ones. Parameters of the pyrolysis reaction and material composition also show decreasing contributions as the pyrolysis reaction comes to an end. In order to make more evident the different contributions, Figure 8 illustrates the different contributions gathered in terms of different groups of uncertainty, i.e. fibers and matrices properties, TACOT's composition, pyrolysis parameters, etc.

Inria

By using the low-cost sensitivity analysis technique, the hierarchy of the most important uncertain parameters contributing to the variance of the temperature can be computed as a function of the depths at a fixed time. For example, in Figure 9 we show the hierarchy at a time of 80 s. As it can be observed, the trend is highly non linear, and this could be particularly useful to build reliable design margin policies and to guide material model development efforts.

Construction of the Polynomial-Chaos based surrogate

As explained in Section 3.3.2, the interest of low-cost sensitivity analysis technique is twofold. More than only identifying a ranking of main uncertainties, it can be used in order to build a surrogate model (a Polynomial-Chaos based one in this case) on a reduced set of uncertainties, i.e. the predominant ones. This is applied here to the temperature computed at a depth of 1.6 cm for different times. Note that this can be easily applied to a whatever quantity of interest, but the surrogate will be not the same since the most important uncertainties can be different with respect to the time, the depth and the quantity of interest.

In the case under consideration, first the uncertainties contributing the most to the variance of the temperature at different times are computed (Figure 10 illustrates the ranking for a time t = 80 s where the variance is maximal). Secondly, the PC-based surrogate is constructed on a set of input parameters defined by the predominant uncertainties. Here, uncertainties labeled as 24, 3, 4, 26, 1, 2, 9, 11 and 22 are chosen as input parameters (See Section 2.2 for identifying each uncertainty).

Obviously, the reduction of the problem yields a loss of accuracy with respect to the statistics computation, which could be estimated. Variance reduction is around 11%, which is mainly due to the high number of significant parameters, which are neglected in the reduced problem.

Using the surrogate model, a whatever post-processing statistical analysis can be done for free. As an example, the Probability Density Function (PDF) of the temperature at different times are computed and represented in Figure 11. As it can be observed, gaussian-like PDF can be observed.

One can see here an efficient method for studying TPS response under uncertainties.

Conclusion

The objective of the study was to propagate both material and pyrolysis gas composition uncertainties on the thermal response of the TPS of a spacecraft during atmospheric entry. We chose to study the entry of Stardust that was the first mission using the new generation of

Inria

The solution of this second order ordinary differential equation is θ * (x) = A exp( s/α x) + B exp(s/α x) (30)

where A and B are determined with the boundary conditions. On the semi-infinite domain, we have

• A = 0 as the temperature has to have a limit when x tends towards infinity.

• θ * (x = 0) = B. From θ(0, t) = T w -T 0 , one obtains B = (T w -T 0 )/s after applying Laplace transform.

Hence, we have

θ * (x) = T w -T 0 s exp(-s/α x) (31)
Returning to the temporal space is done using the space change relation: exp(-a √ s)/s ⇔ 1 -erf (a/(2 √ t), where erf is the error function. After some algebra and returning to variable T , we obtain the physical temperature profile as a function of time

Figure 1 :

 1 Figure 1: Phenomenology of porous carbon/phenolic ablative materials

Figure 2 :

 2 Figure 2: Surface and in-depth temperatures obtained with nominal material and pyrolysis gas composition parameters.

Figure 3 :

 3 Figure 3: Analytical solution for unidimensional transient heat transfer with fixed surface temperature, for a semi-infinite medium of thermal diffusivity 10 -7 m 2 /s. Initial temperature of the body: 300 K; surface temperature: 1646 K. Left: temperature profiles for 1, 10, and 60 seconds. Right: error bars in terms of standard deviation when considering two uncertainties.

Figure 4 :

 4 Figure 4: Variance of the temperature (including the contribution of each uncertainty) at a time of 60 seconds.

Figure 5 :

 5 Figure 5: Recession, char at 2% and virgin at 98% with the ± standard deviation envelope.

Figure 6 :

 6 Figure 6: Total variance divided by contributor, computed for several quantities, as a function of time.

Figure 7 :

 7 Figure 7: Total variance divided by contributor, computed for several quantities, as a function of time.

Figure 8 :

 8 Figure 8: Total variance divided by contributor, computed for several quantities, as a function of time. Here, uncertainties are classified into different groups.

Figure 9 :

 9 Figure 9: Hierarchy of the most relevant uncertain parameters to the temperature variance (sorted: the higher, the bigger Sobol index) as a function of depth at t=80 s.

Figure 10 :

 10 Figure 10: Sorted contributions at a depth of 1.6 cm and at a time of t = 80 s

Table 1 :

 1 Contribution to the variance of each uncertainty and error analysis with respect to the number of points per direction. Absolute and relative differences are also given.

	Unc.	Contrib. q16 (%)	Contrib. q8 (%)	∆ abs	∆ rel
	24	16.9	16.9	0.0	0.0
	3	14.8	14.8	0.0	0.0
	4	14.7	14.8	0.1	6.8•10 -3
	26	9.03	9.01	0.02	2.2•10 -3
	1	8.31	8.32	0.01	1.2•10 -3
	2	8.29	8.3	0.01	1.2•10 -3
	9	7.63	7.64	0.01	1.3•10 -3
	11	5.13	5.13	0.0	0.0
	22	4.05	4.05	0.0	0.0
	8	2.88	2.88	0.0	0.0
	14	2.88	2.88	0.0	0.0
	12	2.48	2.48	0.0	0.0
	15	1.48	1.48	0.0	0.0
	7	1.15	1.15	0.0	0.0
	16	0.166	0.166	0.0	0.0
	18	5.05•10 -2	5.03•10 -2	0.2•10 -3	4.0•10 -3
	10	3.78•10 -2	3.78•10 -2	0.0	0.0
	13	5.81•10 -3	5.81•10 -3	0.0	0.0
	20	6.68•10 -4	6.7•10 -4	0.2•10 -5	3.0•10 -3
	6	7.22•10 -5	7.21•10 -5	0.1•10 -6	1.4•10 -3

  T (x, t) = T 0 + (T w -T 0 ) 1 -erf

	x/ 2 √ √ t α	(32)

Inria

RESEARCH CENTRE BORDEAUX -SUD-OUEST

351, Cours de la Libération

Bâtiment A 29 33405 Talence Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr

RR n°9175

Figure 11: PDF of the temperature at a depth of 1.6 cm at different times.

Appendix A

Detailed demonstrations of the analytical resolution of the transient heat transfer equation given by

where α is a constant diffusivity, appears to be hard to find in text books. We will present its resolution using the Laplace transform for a semi-infinite unidimensional medium, with the following boundary conditions

corresponding to a homogeneous medium for which the surface temperature is suddenly raised. Let us apply a variable change that will simplify the integration: θ(x, t) = T (x, t) -T 0 . Applying the change of variable and the Laplace transform, Eq. 26 rewrites

In the Laplace space, the variable is defined as θ * = ∞ 0 e -st θdt. After permutation (between

) and integration of the first integral, we obtain

Thanks to the change of variable, we have θ(x, 0) = 0 and Eq. 28 simplifies into


Joao Guilherme

Caldas Steinstraesser
email: joao-guilherme.caldas-steinstraesser@inria.fr

Gaspard Kemlin
email: gaspard.kemlin@inria.fr

Antoine Rousseau
email: antoine.rousseau@inria.fr

A domain decomposition method for linearized Boussinesqtype equations

Keywords: Boussinesq-type equations, finite differences scheme, transparent boundary conditions, domain decomposition, interface conditions, Schwarz alternating method AMS Subject Classifications: 65M55

In this paper, we derive discrete transparent boundary conditions for a class of linearized Boussinesq equations. These conditions happen to be non-local in time and we test numerically their accuracy with a Crank-Nicolson time-discretization on a staggered grid. We use the derived transparent boundary conditions as interface conditions in a domain decomposition method, where they become local in time. We analyze numerically their efficiency thanks to comparisons made with other interface conditions.

Introduction

Among the main challenges faced in the mathematical framework of coastal engineering is the study of wave propagation in the nearshore area. One field of research in this topic makes use of the Boussinesq equations for water of varying depth that describe the nonlinear propagation of waves in shallow water. The work of Peregrine [START_REF] Peregrine | Long waves on a beach[END_REF], Green and Naghdi [START_REF] Green | A derivation of equations for wave propagation in water of variable depth[END_REF] laid the basis for many Boussinesq-type equations used nowadays. The dispersion properties of these equations have then been improved by Nwogu [START_REF] Nwogu | Alternative Form of Boussinesq Equations for Nearshore Wave Propagation[END_REF] for practical numerical simulation of ocean wave processes from deep to shallow water. In this paper, we work on the equations derived by Nwogu [START_REF] Nwogu | Alternative Form of Boussinesq Equations for Nearshore Wave Propagation[END_REF] and that can be recalled as follows. Consider a three-dimensional wave field with surface elevation η (x,y,t) over a non-constant water depth h(x,y) and with speed u(x,y,z,t) = (u,v), respectively the speeds along the x and y axis, defined at a reference depth z = z(x,y). Figure 1 provides a sketch where we represented these quantities. With τ being the bottom shear stress, Nwogu [START_REF] Nwogu | Alternative Form of Boussinesq Equations for Nearshore Wave Propagation[END_REF] obtained (1.1), which consists of a continuity equation and a momentum equation. These equations have been used in a C/Matlab program known as the Boussinesq Ocean and Surf Zone (BOSZ) model, developed by Roeber and Cheung [START_REF] Roeber | Boussinesq-type model for energetic breaking waves in fringing reef environment[END_REF].

         η t +∇•[(h+η)u]+∇• z 2 2 - h 2 6 h∇(∇•u)+ z+ h 2 h∇(∇•(hu)) = 0, u t +(u•∇)u+g∇η +z z 2 ∇(∇•u t )+∇(∇•(hu t )) +τ = 0.
(1.1) To simplify the framework, we consider in this paper the 1D equations and ignore the bottom shear stress τ. We also consider a constant flat bottom h = h 0 . Thus, the total height H can be recovered by the relation H = h 0 +η. Then, we perform a linearization around the equilibrium point ( η, ū) = (0, ū), with ū ∈ R, to get η t + ūη x +h 0 u x + hu xxx = 0, u t +gη x + ūu x + hu xxt = 0, (

x y z u v η h H
to which we now refer as the linearized Boussinesq equations. Here, h and h are constants defined by

h = z 2 2 +h 0 z+ h 2 0 3 h 0 , h = z z 2 +h 0 . (1.3)
Using the value z = -0.53753×h 0 (see [START_REF] Nwogu | Alternative Form of Boussinesq Equations for Nearshore Wave Propagation[END_REF]) and h 0 = 1, we have that h and h are both negative. Setting ū = 0, h = 0, h = -ε a small parameter, g = 1 and h 0 = 1 leads to the formulation of the linearized Green-Naghdi equations, for which discrete transparent boundary conditions have been derived by Kazakova and Noble [START_REF] Kazakova | Discrete transparent boundary conditions for the linearized green-naghdi system of equations[END_REF]. In this paper, we focus on the case where ū = 0 for the sake of simplicity, but the conclusions of section 3 remain the same for ū = 0. The equations we are going to focus on in this paper are η t +h 0 u x + hu xxx = 0, (1.4a)

u t +gη x + hu xxt = 0. (1.4b)
The first objective of this paper is to derive transparent boundary conditions for equations (1.4a) - (1.4b). Indeed, this system is set on the whole space R and thus we need to restrict the area of computation to a bounded domain for practical applications. This requires to find suitable boundary conditions. We focus on transparent boundary conditions in order to let waves escaping the domain without any reflection (a phenomenon that we observe for instance with Dirichlet conditions). From a mathematical point of view, we set the problem as follows: given a compactly supported initial data, one derives suitable conditions at the boundaries so that the solution on the bounded domain coincides with the restriction to this domain of the solution computed on the whole domain. In practice, this study can be done in either a continuous or discretized framework. A review of these techniques can be found in [START_REF] Antoine | A Review of Transparent and Artificial Boundary Conditions Techniques for Linear and Nonlinear Schr ödinger Equations[END_REF] where the authors build such conditions for the Schr ödinger equation. In the linear case, the study in the continuous framework is carried out by applying the Laplace transform in time and adapt boundary conditions to keep solutions bounded. The Laplace-inverse transform of these conditions results into non-local in time operators. The adaptation of this technique to the discretized framework uses the Z-transform, which is the discrete equivalent of the Laplace transform. Again, the numerical inversion of the obtained discrete conditions yields non-local operators.

The second and main objective of this paper is to test the efficiency of the discrete transparent conditions as interface conditions in a domain decomposition method with the alternating Schwarz method. The interest of this method lies in the possibilities, with few modifications in the original code, to split the original computational domain and/or couple different models (for instance farshore and nearshore). Our work in this paper can be seen as a step towards more efficient coupling between different numerical models in coastal engineering. The main difficulty is then to find suitable interface conditions to exchange information between the different subdomains so that the convergence is achieved as fast as possible.

The paper is organized as follows. In section 2, we apply the study from [START_REF] Arnold | Discrete Transparent Boundary Conditions for the Schr ödinger Equation[END_REF][START_REF] Besse | Discrete Artificial Boundary Conditions for the Korteweg-de Vries Equation[END_REF][START_REF] Kazakova | Discrete transparent boundary conditions for the linearized green-naghdi system of equations[END_REF] to the linearized Boussinesq equations and obtain discrete transparent boundary conditions. Then, we provide some numerical tests to evaluate the accuracy of such conditions before implementing them in a domain decomposition method. In section 3, we briefly recall the idea of the additive Schwarz method that we are going to use before going on with the adaptation of transparent boundary conditions as interface conditions. We conclude with some numerical experiments to evaluate their efficiency.

Derivation of transparent boundary conditions

First, we introduce the initial value problem that we seek to solve using transparent boundary conditions (TBC):

           η t +h 0 u x + hu xxx = 0, ∀ x ∈ R, t > 0, u t +gη x + hu xxt = 0, ∀ x ∈ R, t > 0, u(x,0) = u 0 (x), η(x,0) = η 0 (x), ∀ x ∈ R, u(x,t) ----→ x→±∞ 0, η(x,t) ----→ x→±∞ 0, ∀ t > 0. (2.1)
The goal of TBC is to find boundary conditions for a finite domain (for instance [0,L]) such that the solution on this domain coincides with the restriction to this domain of the solution of problem (2.1). Hence, the problem we will work on is

           η t +h 0 u x + hu xxx = 0, ∀ x ∈ [0,L], t > 0, u t +gη x + hu xxt = 0, ∀ x ∈ [0,L], t > 0, u(x,0) = u 0 (x), η(x,0) = η 0 (x), ∀ x ∈ [0,L],
+ transparent boundary conditions at x = 0 and x = L. Note that it is possible to decouple equations (1.4a) -(1.4b) to obtain an equation on u only. Taking the cross derivatives of (1.4a) and (1.4b), we get

u tt + hu xxtt -gh 0 u xx -g hu xxxx = 0.
(2.3)

Thus, system (1.4a) -(2.
3) is equivalent to (1.4a) -(1.4b). We can also note a fourth order space derivative. Therefore, we will need to derive four conditions. In this section, we follow the steps proposed in [START_REF] Arnold | Discrete Transparent Boundary Conditions for the Schr ödinger Equation[END_REF][START_REF] Besse | Discrete Artificial Boundary Conditions for the Korteweg-de Vries Equation[END_REF][START_REF] Kazakova | Discrete transparent boundary conditions for the linearized green-naghdi system of equations[END_REF] to derive discrete TBC that are adapted to the discretized problem. We end this section with some numerical tests to analyze the efficiency of the obtained conditions.

Deriving discrete transparent boundary conditions

We introduce here the derivation of discrete TBC (DTBC) for the discretized problem. These are boundary conditions directly derived from the discretized problem. We adapt here the method used in [START_REF] Arnold | Discrete Transparent Boundary Conditions for the Schr ödinger Equation[END_REF] for the Schr ödinger equation, in [START_REF] Besse | Discrete Artificial Boundary Conditions for the Korteweg-de Vries Equation[END_REF] for the KdV equation and in [START_REF] Kazakova | Discrete transparent boundary conditions for the linearized green-naghdi system of equations[END_REF] for the linearized Green-Naghdi equations. It follows four main stpdf:

1. Discretize the equations (1.4a) and (1.4b).

2. Use the Z-transform on the problem on the complementary set.

3. Find the conditions at the two boundaries in the Z-space.

4. Use the inverse Z-transform to find the transparent boundary conditions.

Discretization

The first step is to find a discretization in time on a staggered grid. Let δx be the spatialstep and δt be the time-step. We build this grid such that J = L/δx and x j = jδx:

0 = x 0 < x 1 < x 2 < ••• < x J-1 < x J = L.
j < 0 and j > J denote nodes that are out of the domain we want to work on. The time t n then stands for nδt. Using a Crank-Nicolson scheme, (1.4a) yields

η n+1 j+ 1 2 -η n j+ 1 2 δt + h 0 2 u n+1 j+1 -u n+1 j δx + u n j+1 -u n j δx + h 2 u n+1 j+2 -3u n+1 j+1 +3u n+1 j -u n+1 j-1 δx 3 + u n j+2 -3u n j+1 +3u n j -u n j-1 δx 3 = 0.
(2.4)

Note that the finite differences operator used for the third spatial derivative is centered around j+ 1 2 . Discretizing (1.4b) is straightforward:

u n+1 j -u n j δt + g 2   η n+1 j+ 1 2 -η n+1 j-1 2 δx + η n j+ 1 2 -η n j-1 2 δx   + h δt u n+1 j+1 -2u n+1 j +u n+1 j-1 δx 2 - u n j+1 -2u n j +u n j-1 δx 2 = 0.
(2.5)

Problem on the complementary set

To solve the initial value problem (2.2), we assume that the initial conditions u 0 and η 0 are compactly supported in [0,L]. The derivation of the DTBC associated to our problem can be done by studying the problem on the complementary set of [0,L]:

           η t +h 0 u x + hu xxx = 0, ∀ x ∈ R\[0,L], t > 0, u t +gη x + hu xxt = 0, ∀ x ∈ R\[0,L], t > 0, u(x,0) = 0, η(x,0) = 0, ∀ x ∈ R\[0,L], u(x,t) ----→ x→±∞ 0, η(x,t) ----→ x→±∞ 0, ∀ t > 0.
(2.6)

Z-transform

The second step is to compute the Z-transform of (2.6). The Z-transform of the sequence (u n ) n 0 is defined as a function of the complex variable z:

∀ |z| > R > 0, û(z) = Z {(u n )} := ∑ n 0 u n z -n , (2.7)
where R is the convergence radius of the series. One important property of the Z-

transform is that Z {(u n+1 )} = z û(z)-zu 0 . (2.8)
In our case, as ∀ j 0, j J, u 0 j = 0 and η 0 j = 0 (by compactness of the initial conditions in [0,L]), we have

∀ j 0, j J, Z u n+1 j = z ûj (z) and Z η n+1 j+ 1 2 = z ηj+ 1 2 (z).
(2.9)

From now on, ûj and ηj+ 1 2 will stand for ûj (z) and ηj+ 1 2 (z). Using (2.9) on (2.4) and (2.5), we get

         ηj+ 1 2 = - 1 s(z) h 0 δx ûj+1 -ûj + h δx 3 ûj+2 -3 ûj+1 +3 ûj -ûj-1 , (2.10a)
g δx ηj+ 1 2 -ηj-1 2 + h δx 2 s(z) ûj+1 +s(z) 1-2 h δx 2 ûj + h δx 2 s(z) ûj-1 = 0, (2.10b)
where

s(z) = 2 δt z-1 z+1 . (2.11)
Note that s(z) has a singularity at z = -1. For this reason, we will work only on C 1 := {z ∈ C; |z| > 1}. On this set, (s(z)) >0. Injecting (2.10a) in (2.10b) gives, in a similar way than the derivation of (2.3), the following recurrence (in space) relation for û:

- g h δx 4 ûj+2 + h δx 2 s(z) 2 - gh 0 δx 2 +4 g h δx 4 ûj+1 + s(z) 2 1-2 h δx 2 +2 gh 0 δx 2 -6 g h δx 4 ûj + h δx 2 s(z) 2 - gh 0 δx 2 +4 g h δx 4 ûj-1 - g h δx 4 ûj-2 = 0.
(2.12)

We look for a solution of the form

ûj = 4 ∑ i=1 λ i r i (z) j , (2.13)
where the r i are the roots of the corresponding characteristic polynomial, which reads

P(r) := r 4 + α-βs(z) 2 -4 r 3 + s(z) 2 (2β-γ)-2α+6 r 2 + α-βs(z) 2 -4 r+1, (2.14)
where we used

α = h 0 δx 2 h < 0, β = hδx 2 g h > 0, γ = δx 4 g h < 0. (2.15)
Note that α is adimensional whereas β and γ are homogeneous to s 2 , resulting into the coefficients of P being adimensional.

Solving the recurrence relation

Now, we are interested in the study of the polynomial P, and especially its roots as the explicit formulation of ûj will rely on them.

Proposition 1. For all z ∈ C 1 , the roots of the polynomial P given in (2.14) verify (by ordering them with

|r i | |r i+1 |, i = 1,2,3) |r 4 (z)| |r 3 (z)| > 1 > |r 2 (z)| |r 1 (z)|.
Proof. The proof of this proposition is inspired from [START_REF] Besse | Discrete Artificial Boundary Conditions for the Korteweg-de Vries Equation[END_REF]. First, let us specify that, as the constant coefficient is equal to 1, we know that |r 1 (z)r 2 (z)r 3 (z)r 4 (z)| = 1. Moreover, we can see in (2.14) that the coefficients of P are symmetric, in the sense that, as 0 is not a root, P(r) = 0 ⇔ P(1/r) = 0. From these two facts, we can say that

|r 4 (z)| |r 3 (z)| 1 |r 2 (z)| |r 1 (z)|,
where r 2 = 1/r 3 and r 1 = 1/r 4 . Next, we show that there is no root of modulus 1. To this end, assume there exists a θ such that e iθ is a root of P. Then, e 4iθ + α-βs(z) 2 -4 e 3iθ + s(z) 2 (2β-γ)-2α+6 e 2iθ + α-βs(z) 2 -4 e iθ +1 = 0. (2.16)

Multiplying by e -2iθ gives, after using Euler formulas, 2cos2θ +2cosθ α-βs(z) 2 -4 +s(z) 2 (2β-γ)-2α+6 = 0.

(2.17)

Finally, we find

s(z) 2 = 2α(1-cosθ)+8cosθ -4cos 2 θ -4 2β(1-cosθ)-γ = 2α(1-cosθ)-4(cosθ -1) 2 2β(1-cosθ)-γ 0, (2.18)
so that s(z) ∈ iR. Hence, as z ∈ C 1 ⇒ (s(z)) > 0, we reached a contradiction and P has no root of modulus 1, which concludes the proof. Now that we know the roots of the polynomial P, we can explicit û in both the right and the left regions:

         ûj = 4 ∑ i=1 λ l i r i (z) j , ∀ j 2, ûj = 4 ∑ i=1 λ r i r i (z) j , ∀ j J -2.
(2.19)

As we are aiming at a bounded solution, we can eliminate two components in each region (and drop the r and l):

ûj = λ 3 r 3 (z) j +λ 4 r 4 (z) j , ∀ j 2, ûj = λ 1 r 1 (z) j +λ 2 r 2 (z) j , ∀ j J -2.
(2.20)

Finding the discrete transparent boundary conditions

Let us recall that we need two relations at each interface, as mentioned at the beginning of this section. It can be easily verified by a straightforward computation using (2.20) that, on the left, ûj satisfies the following two conditions (we dropped the z for the sake of clarity):

∀ j 2,        1 r 3 r 4 ûj - r 3 +r 4 r 3 r 4 ûj-1 + ûj-2 = 0, 1 (r 3 r 4 ) 2 ûj+2 -2 r 3 +r 4 (r 3 r 4 ) 2 ûj+1 + (r 3 +r 4 ) 2 (r 3 r 4 ) 2 ûj -ûj-2 = 0.
(2.21)

For the right interface, the conditions are

∀ j J -2, ûj+2 -(r 1 +r 2 ) ûj+1 +r 1 r 2 ûj = 0, ûj+2 -2(r 1 +r 2 ) ûj+1 +(r 1 +r 2 ) 2 ûj -(r 1 r 2 ) 2 ûj-2 = 0. (2.22)
Before applying the inverse Z-transform , we define the nine following kernels:

k 1 := r 1 +r 2 , k 2 := (r 1 +r 2 ) 2 , k 3 := (r 1 r 2 ), k 4 := (r 1 r 2 ) 2 , k 5 := r 3 +r 4 r 3 r 4 , k 6 := (r 3 +r 4 ) 2 (r 3 r 4 ) 2 , k 7 := 1 r 3 r 4 , k 8 := 1 (r 3 r 4 ) 2 , k 9 := r 3 +r 4 (r 3 r 4 ) 2 .
(2.23) Hence, equations (2.21) and (2.22) become

∀ j 2, ûj-2 -k 5 ûj-1 +k 7 ûj = 0, ûj-2 -k 6 ûj +2k 9 ûj+1 -k 8 ûj+2 = 0, (2.24) ∀ j J -2, ûj+2 -k 1 ûj+1 +k 3 ûj = 0, ûj+2 -2k 1 ûj+1 +k 2 ûj -k 4 ûj-2 = 0. (2.25)
Finally, by applying the inverse Z-transform at j = 2 and j = J -2, we find two conditions at each interface, with {u n } = u 0 ,...,u n :

Left    Γ l 1 ({u n }) := u n 0 -(Y 5 * u 1 ) n +(Y 7 * u 2 ) n = 0, Γ l 2 ({u n }) := u n 0 -(Y 6 * u 2 ) n +2(Y 9 * u 3 ) n -(Y 8 * u 4 ) n = 0, (2.26) Right    Γ r 1 ({u n }) := u n J -(Y 1 * u J-1 ) n +(Y 3 * u J-2 ) n = 0, Γ r 2 ({u n }) := u n J -2(Y 1 * u J-1 ) n +(Y 2 * u J-2 ) n -(Y 4 * u J-4 ) n = 0, (2.27)
where Y n i := Z -1 {k i (z)} for i = 1,...,9 and * denotes the discrete convolution:

∀ i ∈ 1,9 , (Y i * u j ) n := n ∑ m=0 Y m i u n-m j .
(2.28)

It is possible to compute the convolution coefficients Y i with an inverse Fast Fourier transform, as described in [START_REF] Arnold | Discrete Transparent Boundary Conditions for the Schr ödinger Equation[END_REF][START_REF] Besse | Discrete Artificial Boundary Conditions for the Korteweg-de Vries Equation[END_REF]. Indeed, the explicit formula for the inversion of û(z) is

u n = 1 2iπ C r û(z)z n-1 dz, r > R, (2.29)
where C r denotes any circle of radius r > R. Applying this formula to the kernels k i and discretizing gives

∀ i ∈ 1,9 , n ∈ 0, N -1 , Y n i = r n N N-1 ∑ k=0 k i rω k N ω nk N = r n F -1 k i rω k N , (2.30)
where F -1 is the inverse discrete Fourier transform, N is the number of nodes to discretize the circle and ω N := e i 2π N . This results into non-local in time operators Γ l,r 1,2 but only at some points near the boundaries and [START_REF] Arnold | Discrete Transparent Boundary Conditions for the Schr ödinger Equation[END_REF][START_REF] Besse | Discrete Artificial Boundary Conditions for the Korteweg-de Vries Equation[END_REF] give techniques to transform the convolution into recurrence relations that can be updated at each time-step. The goal of this paper being to study the application of transparent boundary conditions to domain decomposition methods, we will use the convolution coefficients as they are.

Remark 2.1. Some Y i are represented in figure 3. Instabilities can be observed on the Y i for large n (they are not represented here). However, we will see when applying the DTBC to domain decomposition techniques in section 3 that only the first coefficients Y 0 i are needed. The interested reader can refer to [START_REF] Besse | Discrete transparent boundary conditions for the mixed KDV-BBM equation[END_REF], where a better sampling is introduced to avoid these instabilities for large n.

Final numerical scheme

We now have all the tools in our hands to derive a global numerical scheme. To emphasize the fact that the operators at the boundaries Γ l,r 1,2 act on u, we decide to solve numerically the decoupled equation (2.3), that we recall here:

u tt + hu xxtt -gh 0 u xx -g hu xxxx = 0.
Using a Crank-Nicolson scheme, we get the scheme (2.39). Note that we get a 4 th order differential equation, which resulted into a 5-points stencil. As it is also a two time-stpdf 10 -6 10 -5 scheme, we need to specify how to compute u 1 j knowing only the initial conditions u 0 and η 0 . This can be done with a Taylor extension in the vicinity of t = 0:

|Y n 1 | n -3/2 |Y n 1 | n -3/2
10 -6 |Y n 9 | n -3/2 |Y n 9 | n -3/2
u+ hu xx t=δt = u+ hu xx t=0 +δt u+ hu xx t t=0 + δt 2 2 u+ hu xx tt t=0 +o(δt 2 ). (2.31)
As u(•,0) ≡ 0, the first term is 0. Then, using equation (1.4b), we find that Finally, we get the approximation

u+ hu xx t t=0 = -gη x | t=0 . ( 2
u+ hu xx t=δt = -gη x | t=0 +o(δt 2 ), (2.34)
that we can discretize to compute u 1 . One may thus wonder what conditions to use so that we also have TBC at the first time-step? The idea is to apply again the same method and solve the problem on the complementary set (as η x (•,0) is compactly supported in [0,L]):

u+ hu xx = 0, ∀ x ∈ R\[0,L], u(x) ----→ x→±∞ 0. (2.35)
Given the conditions at infinity, we have

u(x) = c 1 e 1 √ - h x , ∀ x < 0, u(x) = c 2 e -1 √ - h x , ∀ x > L, (2.36)
which gives the following boundary conditions:

       u (0) = 1 √ - h u(0), u (L) = - 1 √ - h u(L).
(2.37) Knowing u 0 and η 0 , the final numerical scheme is:

• n = 0                            u 1 0 --h -u 1 0 +u 1 1 δx = 0, u 1 j + h u 1 j+1 -2u 1 j +u 1 j-1 δx 2 = -gδt η 0 j+ 1 2 -η 0 j-1 2 δx , ∀ 1 j J -1, u 1 J + -h -u 1 J-1 +u 1 J δx = 0.
(2.38) Remark 2.3. The numerical scheme (2.39) can be used to find u only, as all our study relied only on equations where we had eliminated the unknown η (except to perform the first time-step). Note that at each iteration in time, we can recover η using (2.4). It is also possible to solve the system (2.4) -(2.5) directly with the same transparent boundary conditions (2.26) -(2.27) and find u and η simultaneously. The results are the same and the main advantage is that, as it is first order system in time, we don't need to perform the first step independently. Using X n =(η n ,u n )∈R 2N-1 , where N is the number of nodes in the discretization, the numerical system reads

• n 1                                                          Γ 1 u n+1 = 0, Γ l 2 u n+1 = 0, u n+1 j -2u n j +u n-1 j δt 2 + h δt 2 u n+1 j+1 -2u n+1 j +u n+1 j-1 δx 2 -2 u n j+1 -2u n j +u n j-1 δx 2 + u n-1 j+1 -2u n-1 j +u n-1 j-1 δx 2 - gh 0 4 u n+1 j+1 -2u n+1 j +u n+1 j-1 δx 2 +2 u n j+1 -2u n j +u n j-1 δx 2 + u n-1 j+1 -2u n-1 j +u n-1 j-1 δx 2 - g h 4 u n+1 j+2 -4u n+1 j+1 +6u n+1 j -4u n+1 j-1 +u n+1 j-2 δx 4 +2 u n j+2 -4u n j+1 +6u n j -4u n j-1 +u n j-2 δx 4 + u n-1 j+2 -4u n-1 j+1 +6u n-1 j -4u n-1 j-1 +u n-1 j-2 δx 4 = 0, ∀ 2 j J -2, Γ r 1 u n+1 = 0, Γ r 2 u n+1 = 0. ( 2
1 δt Id A B C X n+1 = 1 δt Id -A -B Id 0 C X n + 0 R N-1 b n , ( 2

.40)

with

A:=         -h 2δx 3 -h 0 2δx 3 h 2δx 3 + h 0 2δx -3 h 2δx 3 h 2δx 3 -h 2δx 3 3 h 2δx 3 -h 0 2δx -3 h 2δx 3 + h 0 2δx h 2δx 3 (0) . . . . . . . . . . . . (0) -h 2δx 3 3 h 2δx 3 -h 0 2δx -3 h 2δx 3 + h 0 2δx h 2δx 3 -h 2δx 3 3 h 2δx 3 -3 h 2δx 3 -h 0 2δx h 2δx 3 + h 0 2δx         ∈R (N-1)×N , B:=            0 0 0 (0) -g 2δx g 2δx . . . . . . -g 2δx g 2δx (0) 0 0 0            ∈R N×(N-1) , Id 0 :=            0 0 (0) 1 . . . 1 (0) 0 0            ∈R N×N , C :=                     1 -Y 0 5 Y 0 7 1 0 -Y 0 6 2Y 0 9 -Y 0 8 h δtδx 3 1 δt -2 h δtδx 3 h δtδx 3 (0) . . . . . . . . . ( 0
) h δtδx 3 1 δt -2 h δtδx 3 h δtδx 3 -Y 0 4 0 Y 0 2 -2Y 0 1 1 Y 0 3 -Y 0 1 1                     ∈ R N×N , b n :=                    (Y 5 * u 1 ) n -(Y 7 * u 2 ) n (Y 6 * u 2 ) n -2(Y 9 * u 3 ) n +(Y 8 * u 4 ) n 0 . . . 0 2(Y 1 * u J-1 ) n -(Y 2 * u J-2 ) n +(Y 4 * u J-4 ) n (Y 1 * u J-1 ) n -(Y 3 * u ) n                    ∈ R N ,
where

∀ i ∈ 1,9 , (Y i * u j ) n := n ∑ m=1 Y m i u n-m j .

Numerical results

We present here some numerical results to analyze the efficiency of the conditions we just derived. We used the following values for the different parameters we can vary ( h and h resulted from h 0 and z):

L = 1m, g = 9.81m•s -2 , δt = 0.001s, δx = 0.01m, h 0 = 1m, z = -0.53753×h 0 , h = -0.39306m 2 , h = -0.05973m 3 .
(2.41)

For the initial conditions, we used

u 0 ≡ 0, η 0 (x) = exp -400× x- L 2 2 .
(2.42)

Finally, we used r=1.001 for the computation of the Y i in (2.30). The results are presented in figure 4 and table 1. They are quite satisfying as the solution computed with DTBC fits well the reference solution computed on a bigger domain. We used the error e n as relative error at time t n and e T as global l 2 error in time:

e n = u ref (•,t n )-u num (•,t n ) l 2 u ref (•,t n ) l 2 , e T = δt× n max ∑ n=1 (e n ) 2 .
(2.43)

We used a trapezoidal rule to compute the discrete l 2 -norm in space and the reference solution is computed on a bigger domain, with the same transparent conditions. We also encountered some conditioning issues when solving the linear system, but using a Jacobi pre-conditioner is enough to get a correct solution. The error is expected to be small by definition of the transparency of the conditions. Table 1: e n for different times and e T .

• • • • • • • • • • • • • • • • • x 0 0.2 0.4 0.6 0.8 1 u(•, t) -0.15 -0.1 -0.05 0 0.05 0.1 0.15 t = 0.25 uref u • uref u • • • • • • • • • • • • • • • • • • x 0 0.2 0.4 0.6 0.8 1 u(•, t) -0.15 -0.1 -0.05 0 0.05 0.1 0.15 t = 0.5 uref u • uref u • • • • • • • • • • • • • • • • • • x 0 0.2 0.4 0.6 0.8 1 u(•, t) -0.15 -0.1 -0.05 0 0.05 0.1 0.15 t = 0.75 uref u • uref u • • • • • • • • • • • • • • • • • • x 0 0.2 0.4 0.6 0.8 1 u(•, t) -0.15 -0.1

Partial conclusion

In this section, we generalized the work of Kazakova and Noble [START_REF] Kazakova | Discrete transparent boundary conditions for the linearized green-naghdi system of equations[END_REF] to the case of the linearized Boussinesq equations. The study is similar but what is new here is that the additional dispersive term hu xxx yields a polynomial of higher order and more conditions at the boundaries. The derived TBC revealed themselves, as expected, non-local in time.

Then, we performed numerical experiments where we noticed that the conditions are really satisfying. In the next section, we will focus on the application of these conditions to domain decomposition methods.

Application to domain decomposition methods

The discrete boundary conditions (2.26) and (2.27) will be used in this section as interface boundary conditions (IBC) in a domain decomposition method (DDM). We briefly describe the DDM that we are going to use before presenting analysis and numerical results.

The Schwarz method

Domain decomposition methods are used to split a domain Ω, on which we want to solve a given problem, in multiple domains Ω i , that can possibly overlap. Then, we can solve the problem in each domain. Hence, one must find functions that satisfy the PDE in each domain and that match with its neighbours on the interfaces, in a sense that has to be defined. The main difficulty of domain decomposition methods lies in the definition of efficient conditions at the interface between subdomains. The original DDM was developed by Schwarz in 1870 [START_REF] Schwarz | Ueber einen Grenz übergang durch alternirendes Verfahren[END_REF] for a stationary problem and consists in an iterative method: the solution in the i-th subdomain Ω i is computed as the limit of a sequence u k i , k 0. At each iteration k, we solve the problem in each subdomain with boundary conditions at the interfaces imposed using functions from the other subdomains. We will consider here the additive Schwarz method (ASM) in which the interface conditions is always constructed using solutions u k-1 j (j = i) from the previous iteration in the neighbour subdomains. Therefore, at each interface between two subdomains Ω i and Ω j , the IBC in Ω i is

Ω 1 ∂Ω ext 1 Mono-domain BC Ω 2 ∂Ω ext 2 Mono-domain BC Γ Interface BC
B i u k+1 i = B i u k j . (3.1)
Note that it is possible to impose several interface conditions when Ω i has several neighbours. Initially, the operators B i were Dirichlet conditions: B i (u)= u. For more details on the Schwarz method, the reader can refer to [START_REF] Gander | Schwarz methods over the course of time[END_REF][START_REF] Lions | On the Schwarz Alternating Method I[END_REF][START_REF] Lions | On the Schwarz Alternating Method III: A variant for Nonoverlapping Subdomains[END_REF]. We now look for an operator B i inspired from the DTBC we derived in the previous section. Without loss of generality, we consider the domain Ω = [0,L] divided into two subdomains Ω 1 (left subdomain) and Ω 2 (right subdomain) that can possibly overlap (cf. figure 6). Ω is discretized into N nodes, while Ω 1 and Ω 2 are discretized into N 1 and N 2 nodes. The nodes in the subdomains coincide with the nodes in Ω. The interface conditions B 1 and B 2 will have to be such that:

• At each Schwarz iteration k, there is a unique discrete solution u n,k i , at each timestep t n , in each subdomain Ω i .

• The solution in each subdomain will have to converge (in the sense of Schwarz) to the solution on Ω restricted to this subdomain.

Note that, as problem (2.2) is a time-dependent problem, we will perform the Schwarz method at each time-step. There are other approaches for solving time-dependent PDEs with Schwarz methods, e.g. the Schwarz waveform relaxation method [START_REF] Gander | Overlapping schwarz waveform relaxation for the heat equation in n dimensions[END_REF][START_REF] Halpern | Optimized schwarz waveform relaxation and discontinuous Galerkin time stepping for heterogeneous problems[END_REF][START_REF] Halpern | Optimized and quasi-optimal Schwarz waveform relaxation for the one-dimensional Schr ödinger equation[END_REF], in which the communications between the subdomains are made after the full simulations in time. Nevertheless, as shown here in the sequel, our local-in-time DDM has a very fast convergence, making it unnecessary to use global-in-time approaches, which would be more costly.

To avoid any confusion between the iteration in the Schwarz method and our timedependent problem, the word iteration and the integer k will refer to the Schwarz algorithm whereas time-step and n, m will refer to the evolution in time of the problem we are solving.

A Schwarz method with transparent conditions at the interface

From [START_REF] Japhet | The Best Interface Conditions for Domain Decomposition Methods: Absorbing Boundary Conditions[END_REF], we know that transparent boundary conditions are very good candidates for interface conditions in DDM. Thus, the TBC we derived in the previous section might inspire us to set up interface conditions in an additive Schwarz method. We recall that we look for discrete interface conditions, so that all the future reasonings will be done with the discrete equation and conditions (2.39). If we focus on the right interface of the left domain, we recall that the transparent boundary conditions is

   Γ r 1 u n+1 1 = 0, Γ r 2 u n+1 1 = 0. ( 3

.2)

A first heuristic is to use the left-hand-side with u n+1,k 2 to provide a right-hand-side for the computation of u n+1,k+1

1 . With u n+1 1 k = u 0 1 ,...,u n 1 ,u n+1,k 1
where all the u m 1 for 0 m n

1 = B 1 u n+1,k 2 . ( 3.5)
It is worth noting that the operator (3.5) used in (3.4) is local in time, contrarily to the operators Γ r 1,2 . Similarly, the IBC for the left boundary of the right domain is

     u n+1,k+1 2,0 -Y 0 5 •u n+1,k+1 2,1 +Y 0 7 •u n+1,k+1 2,2 = u n+1,k 1,O 1,2 -Y 0 5 •u n+1,k 1,O 1,2 +1 +Y 0 7 •u n+1,k 1,O 1,2 +2 , u n+1,k+1 2,0 -Y 0 6 •u n+1,k+1 2,2 +2Y 0 9 •u n+1,k+1 2,3 -Y 0 8 •u n+1,k+1 2,4 = u n+1,k 1,O 1,2 -Y 0 6 •u n+1,k 1,O 1,2 +2 +2Y 0 9 •u n+1,k 1,O 1,2 +3 -Y 0 8 •u n+1,k 1,O 1,2 +4 , (3.6) where J 2 = N 2 -1 and O 1,2 = N -N 2 is
the index such that the node O 1,2 of the left domain coincides with the node 0 of the right domain. Again, we rewrite this system as

B 2 u n+1,k+1 2 = B 2 u n+1,k 1 .
(3.7)

Remark 3.1. Let us notice in (3.4) and (3.6) that the interface conditions need at least 5 nodes in the overlap zone and that it does not allow a non-overlapping ASM. It is not really an issue as the number of nodes required in the overlap zone is small and does not depend on the mesh size.

Finally, the ASM along with these interface conditions reads

                 Γ l 1 u n+1 1 = 0, Γ l 2 u n+1 1 = 0, left, u n+1,k+1 1 = f 1 u n 1 ,u n-1 1 , interior, B 1 u n+1,k+1 1 = B 1 u n+1,k 2 , right,                  B 2 u n+1,k+1 2 = B 2 u n+1,k 1 , left, u n+1,k+1 2 = f 2 u n 2 ,u n-1 2 , interior, Γ r 1 u n+1 2 = 0, Γ r 2 u n+1 2 = 0, right, (3.8
) where f i stands for the discrete equation from the numerical scheme (2.39) applied to the interior of domain i.

Numerical results

We end this paper with numerical results obtained with the implementation of an additive Schwarz method to solve the linearized Boussinesq equations. We used the same parameters than in the previous tests (2.41). What we are interested in here is the number of Schwarz iterations required by the DDM (3.8) to converge to the reference solution (given by the solution computed on the domain Ω with discrete transparent conditions at the boundaries). We will use the following stopping criterion:

e n,k DDM ε, ( 3.9)
where ε = 10 -12 and e n,k DDM = δx To study the efficiency of this method, we analyse its convergence at fixed points in time and compare our conditions at the interface to the classical Dirichlet conditions.

Results are presented in figure 7 and table 2 for the minimum overlap size. They are very satifying as we can see that our interface conditions make the Schwarz algorithm converge in 2 iterations only, where the Dirichlet interface conditions require several hundreds of iterations. As these conditions correspond to the fastest possible convergence (2 iterations), it is not necessary to study the influence of the size of the overlap zone nor the efficiency of a global Schwarz algorithm as it will not converge in less than 2 iterations.

Conclusion

In this paper, we have been using discrete transparent boundary conditions for a class of Boussinesq equations. As expected these conditions (non-local in time) provide very satisfying results with respect to the wave reflection at boundaries. When implemented in a domain decomposition framework, the new conditions happen to be both very efficient and local in time, which provides the best possible framework for the simulation in large domains using decomposition techniques.

In future works, we shall both adapt these results to nonlinear and 2D Boussinesq-type equations and study the extension of domain decomposition algorithms to coupling techniques (e.g. to simulate a dispersive/non-dispersive transition).

Figure 1 :

 1 Figure 1: Definition of the quantities η, h, H, u.

Figure 2 :

 2 Figure 2: Transparent boundary conditions means that the solution we want to compute on [0,L] is a picture of the solution on R restricted to [0,L].

Figure 3 :

 3 Figure 3: Coefficients Y 1 and Y 9 used in the convolution of the DTBC (2.26) and (2.27), with parameters from (2.41). We can observe a decreasing rate of O n -3/2 . The behavior of the other coefficients is similar.

  .32) Combining (1.4a) and (1.4b), we also have u+ hu xx tt t=0 = -(gη xt )| t=0 = gh 0 u xx +g hu xxxx t=0 = 0.(2.33)

Figure 4 :

 4 Figure 4: Snapshot at different times of the reference solution u ref (computed on a bigger domain with the transparent conditions at the boundaries) and the solution u computed with the same transparent conditions.

Figure 5 :

 5 Figure 5: Decomposition of a domain Ω = Ω 1 ∪Ω 2 into two subdomains without overlapping. On ∂Ω = ∂Ω ext 1 ∪∂Ω ext 2 , we use the boundary conditions of the mono-domain problem. On Γ = ∂Ω 1 ∩∂Ω 2 , interface boundary conditions have to be defined.

2 Figure 6 :

 26 Figure 6: Decomposition of Ω into two domains with overlapping. Blue (green) nodes correspond to the left (right) domain. A filled circle represents the TBC when they are used on the external boundary (non-local in time) whereas an empty circle represents a TBC used at the interface between the two domains (local in time).

  1,2 is the index on the mono-domain of the node corresponding to the first node of the right domain.

Figure 7 :

 7 Figure 7: Evolution of e n,kDDM at different times. We used, for the subdomains, N 1 = 88 and N 2 = 18 so that the overlap zone is of size 5 (we have N = 101 from (2.41)).

Remark 3 . 2 .

 32 We also tested conditions (2.26) -(2.27) (designed for ū = 0) to solve numerically (1.2) when ū = 0. If they are not adapted when used as transparent boundary conditions, they are still very efficient when used as interface conditions in a DDM: the convergence of the Schwarz method with (3.4) -(3.6) and ū = 0 is reached in 3 iterations.

Table 2 :

 2 Number of iterations required for convergence in the ASM. We used, for the subdomains, N 1 =88 and N 2 = 18 so that the overlap zone is of size 5 (we have N = 101 from (2.41)).

Acknowledgements

The authors thank E. Blayo for fruitful discussions and helpful comments. This work was partially supported by the Marine Energy Research & Innovation Center (MERIC) project CORFO 14CEI2-28228 and the Fundaci ón Inria Chile.

have been computed by the Schwarz algorithm at the previous time-stpdf, we have

is the index such that the node J 2,1 of the right domain coincides with the node J 1 = N 1 -1 of the left domain. The operator defined on each side of this system will be denoted as B 1 , so that the IBC becomes