MULTIGRID PRECONDITIONING OF THE STEAM GENERATOR TWO-PHASE MIXTURE BALANCE EQUATIONS IN THE GENEPI SOFTWARE

Keywords: Two-phase flows, Steam generator, Multigrid, GENEPI

In the framework of the two-phase fluid simulations of the steam generators of pressurized water nuclear reactors, we present in this paper a geometric version of a pseudo-Full MultiGrid (pseudo-FMG) Full Approximation Storage (FAS) preconditioning of balance equations in the GENEPI code. In our application, the 3D steady state flow is reached by a transient computation using a semi-implicit fractional step algorithm for the averaged two-phase mixture balance equations (mass, momentum and energy for the secondary flow). Our application, running on workstation clusters, is based on a CEA code-linker and the PVM package. The difficulties to apply the geometric FAS multigrid method to the momentum and mass balance equations are addressed. The use of a sequential pseudo-FMG FAS twogrid method for both energy and mass/momentum balance equations, using dynamic multigrid cycles, leads to perceptibly improvements in the computation convergences. An original parallel red-black pseudo-FMG FAS three-grid algorithm is presented too. The numerical tests (steam generator mockup simulations) underline the sizable increase in speed of convergence of the computations, essentially for the ones involving a large number of freedom degrees (about 100 thousand cells). The two-phase mixture balance equation residuals are quickly reduced: the reached speed-up stands between 2 and 3 following the number of grids. The effects on the convergence behavior of the numerical parameters are investigated.

INTRODUCTION

This paper is devoted to the presentation of the multigrid preconditioning of the mixture momentum balance equations in the CEA GENEPI software [START_REF] Obry | An advanced steam generators design 3D code[END_REF] dedicated to the Steam Generator (SG) simulation. It is based on a geometric pseudo-full multigrid (FMG) version of the Full Approximation Storage method (FAS) to accelerate the solve of the mixture balance equations. The reader can get more details in the reference [START_REF] Belliard | Multigrid preconditionning in the genepi software[END_REF] and [START_REF] Belliard | Multigrid preconditioning of the steam generator two-phase mixture balance equations[END_REF].

The GENEPI code solves the balance equations of an equivalent mixture in a porous media. The following strong formulations of the mixture balance equations are issued from a homogenization process, [START_REF] Grandotto | Calculs des écoulements diphasiques dans les échangeurs par une méthode aux éléments finis[END_REF]:

1. mass balance (1) ✁ ∇ ✂ ✄ β ✁ G ☎ ✝✆ 0 2. momentum balance βρ∂ t ✁ v ✞ βρ ✄ ✁ v ✂ ✁ ∇ ☎ ✁ v ✞ div ✄ βx ✄ 1 x☎ ρ ✁ v R ✟ ✁ v R ☎ ✆ βρ ✁ g β Λρ ✁ v β ✁ ∇ P ✞ div ✄ βµ T ✄ ✁ ∇ ✁ v ✞ ✁ ∇ t ✁ v ☎ ✠☎ (2)
3. energy balance (enthalpy)

βρ∂ t H ✞ β ✄ ✁ G ✂ ✁ ∇ ☎ H ✞ div ✄ βx ✄ 1 x☎ ρL ✁ v R ☎ ✆ βQ ✞ div ✄ βχ T ✁ ∇ H ☎ (3)
Concerning the mixture energy and the mixture momentum balance equations ( 2) -( 3) , the steadystate flow regime is reached by mean of a pseudo-time marching. No time marching is applied for the mixture mass balance equation (1). Hence this approach leads to an algorithm very close to those used in the incompressible fluid dynamic framework: at each pseudo-time step, a Chorin like scheme [START_REF] Gresho | On the theory of semi implicit projection methods for viscous incompressible flow and its implementation via finite element method that also introduces a nearly consistent matrix. i, theory[END_REF] allows the simultaneous computation of the mixture mass flux and the mixture pressure. During a pseudo-time step, the balance equations are successively solved: energy, then momentum/mass balance equations.

Space discretization is done by mean of a Galerkine finite element method (FEM) leading to a weighted integral version of the above equations (weak formulation) in which the mechanical stress term and the energy diffusion one are integrated by part. The unknowns are : H (Q1), P (Q0) and ✁ v (Q1). The porosity β and the mass flux

✁ G ✆ ρ ✁
v are taken in Q1. Generally, the other physical quantities (i.e. ρ, x, µ T , Λ, ...) are taken in Q0 (i.e. by element). According to the hyperbolic type of the flow equations, Dirichlet boundary conditions are used at the inlets of the domain (mass flux and enthalpy) and Neumann ones at the outlets (pressure). The other boundaries of the domain are adiabatic and impermeable walls. At each pseudo-time step, the arising linear systems are partially solved (5 to 20 iterations) by an iterative method. In order to compute ρ, x and L in function of H and P, we need water thermodynamic tables and the µ T , χ T , Λ, ✁ v R terms are obtained by the use of a large set of semi-empirical closure relations [START_REF] Obry | An advanced steam generators design 3D code[END_REF]. The heat source Q in the enthalpy equation is linked to the resolution of an energy balance equation for the primary flow. To evaluate this term, other correlations on the heat exchange coefficient and the wall temperature are included.

The paper is organized as follow. The Section two is devoted to the presentation of the multigrid correction scheme and of some specific points concerning the FAS correction implementation for the mixture balance equations (energy, momentum and mass). The implementation itself is the object of the Section three. Finally, I present some numerical test cases in the Section four, followed by some concluding remarks.

In the framework of the weighted residual method, a typical problem is to find s (s V , V is a Hilbert space) solution of :

(4) r I ✄ s☎ :✆ ✁ T ✄ s☎ φ I ✄ x☎ dv ✆ 0 ✂ φ I V
where {φ I } is a basis of V , T ✄ s☎ a residual function of the solution s and r I ✄ s☎ the weighted integral of the residual. If u is an approximation of the solution s, the error is defined by e :✆ s u. Roughly speaking, a multigrid technique is an iterative method to solve this equation on a given grid Ω 0 (here, called the finest grid), based on successive estimations of the error on a hierarchy of nested coarser grids Ω l , 0 ✄ l ☎ l max where l max is the index of the coarsest grid. It uses a nested sequence of two-grid methods.

In a two-grid method, we perform only some passes of an iterative method on the fine grid Ω 0 . Then, on the coarse grid Ω 1 , an error equation is formed, involving the restricted fine grid residual ¯r 1 :✆ P 0 1 r 0 (and the restricted fine grid approximation ¯u 1 :✆ P 0 1 u 0 in the case of non-linear problems). Then, a coarse grid error is obtained. It is defined by e 1 :✆ u 1 ¯u 1 . This one is either exact or only estimated. This coarse grid error is then prolonged to the fine grid in order to correct the former estimation of the solution (eventually with a relaxation coefficient α): u 0 ✆ u 0 ✞ αP 1 0 e 1 . This relaxed scheme can be also written as:

(5) u 0 ✆ ✄ 1 αP 1 0 P 0 1 ☎ u 0 ✞ αP 1 0 u 1 ✂
In a FAS multigrid method, on each coarse grid Ω l , 0 ✄ l ☎ l max , the error equation is :

(6) ✁ T l ✄ u l ☎ φ I l ✄ x☎ dv ✆ S I l ✂ φ I l V l with : (7) S I l :✆ ✁ T l ✄ ¯u l ☎ φ I l ✄ x☎ dv ✞ ¯r l
There are two ways to construct the coarse grid operators T l ✄ ☎

. The first one is to build it only from the fine grid (Ω l✝ 1 ) equation discretization (arithmetic multigrid version). The second one is based on the discretizations of the balance equations on the coarse grid Ω l itself (geometric multigrid version). In this paper, we only consider the geometric multigrid version. Now, we perform some analyses concerning some specific aspects of the introduction of a geometric FAS method for the solving of the mixture balance equations (energy, momentum and mass). They are related to the strong variations of the forcing terms following the space discretization step, to the specific application of the mass flow Dirichlet condition near the walls and the use of the FEM and the Chorin's projection algorithm.

For the full mixture balance equation FAS preconditioning, we restrict the mixture mass flux, pressure and enthalpy, the primary fluid temperature and the mixture momentum and energy balance residuals. As explained here after, the restriction of the mixture mass balance residual is not necessary. Let notice that we choose to restrict the mixture mass flux instead of the mixture velocity for a compatibility reason with the divergence free mass flux constraint. A nodal variable u l ✝ 1 is restricted by a canonical restriction (in this case, the coarse node I and the fine one i are the same nodes) and a mixture weighted residual by a nodal weighted average restriction. The element variables u e l ✞ 1 l ✝ 1 (e.g. mixture pressure) are transferred by mean of an arithmetic average or by a volume weighted restriction. The errors concerning the mixture specific enthalpy, the mass flux, the pressure and the primary fluid temperature must be computed and prolonged. The prolongation operator used for the nodal quantities is the trilinear interpolation. For the element quantities, we use a direct affectation.

In our geometric version of the FAS multigrid method, the porosity field is independently computed on each grid. We use an unique U-tube bundle description and several computation domain meshes. Consequently, the porosity field β l and β l✝ 1 are different in discontinuity regions (boundaries of the U-tube bundle).

Divergence Free Mass Flux

The pressure evolution in the Chorin-Gresho algorithm [START_REF] Gresho | On the theory of semi implicit projection methods for viscous incompressible flow and its implementation via finite element method that also introduces a nearly consistent matrix. i, theory[END_REF]) is based on successive updates of pressure distributions associated at the divergence free space projections. If we denote

✁ G := ρ ✄ H n ✁ P n ☎ ✁
v based on the solution ✁ v of the mixture momentum balance equation (2) using the pressure value taken at the previous pseudo-time step P n , we have:

(8) ✁ G i✂ n ✄ 1 ✆ ✁ G i✂ ∑ e ☎ λ e ✁ ∇ ✄ βφ i ☎ dv ☎ βφ i dv and (9) P e✂ n ✄ 1 ✆ P e✂ n 2λ e δt
with λ e solution of a weak formulation of the following equation:

(10)

✁ ∇ ✂ β ✁ ∇ λ ✆ ✁ ∇ ✂ β ✁ G ✂
This scheme is also used to satisfied the free divergence constraint for a given mass flux field. For each pseudo-time step, the divergence free constraint equation ( 1) is solved to the computer precision. Hence, the term ¯r l ✆ P l✝ 1 l r l ✝ 1 in Equation ( 7) is zero. But, the mass flow restricted on the coarse grid ✁ Ḡl is not generally divergence free. Hence, the coarse grid FAS corrected mass balance equation is no longer compatible with the the coarse grid mass flow divergence free condition. So, the pressure Poisson equation (10) must be adapted. For the strong formulation of the balance equations, it becomes:

(11) ✁ ∇ ✂ β l ✁ ∇ λ l ✆ ✁ ∇ ✂ β l ✁ G l ✁ ∇ ✂ β l ✁ Ḡl ✂

Boundary Conditions

Usualy, the user specifies the mixture mass flow rate Q in and the mixture specific enthalpy H in at the inlets and the pressure is specified P out at the outlet. Concerning the outlet pressure and the inlet enthalpy, the locations and the values of the boundary conditions are similar for all the grids. It is not longer true for the inlet mass fluxes, since similar mass flow rates for all the grids lead to different inlet average mixture velocities (the inlet porosities may change). For each pseudo-time step, using the mixture density distribution ρ ✄ H in ✁ P☎ , new inlet mixture velocity values are deduced:

(12) ✁ v in ✆ Q in ρ ☎ S inlet βdS ✁ n ✂
and assigned to the inlet nodes. To avoid discrepancies in the boundary conditions between the several grids, the values of the fine grid inlet mass fluxes are restricted on the coarse grid inlet nodes.

Pressure Computation

As previously mentioned, the pressure evolution in the Chorin-Gresho algorithm is based on successive updates of the initial pressure distribution. This initial pressure must be in coherence with the initial mass flux (Consistent Pressure Poison Equation -CPPE-solver). Then, a pressure update is performed at each divergence free space projection. In the case of the FAS algorithm, two error terms are potentially prolonged on the fine grid: the mass flux error and the pressure error. Simultaneously applying these two error corrections breaks the coherence between the velocity and the pressure. Indeed, we apply the mass flux error correction and then compute a new fine grid pressure field solving the CPPE, before running the Picard iterations.

Forcing Terms

As previously mentioned, the porosity fields may be different on each grid: the β l field is build by intersecting several meshes of the inner technological devices (U-tube bundle, baffles, ...) with the computation domain mesh Ω l . Consequently, following the geometrical scale addressed, the forcing terms may be strongly different. For example, friction forces may be apply to a coarse grid node, but not to the equivalent fine grid node, see Fig. 1. Friction forces are induced by the U-tube bundle or the baffles included in the cells. Hence nodal friction forces are spread on the whole nodes belonging to cells "with friction forces" (see Fig. 1).

The forcing term formulation on a given grid has no direct relation with the formulations on the other grids (geometrical multigrid version). Hence, the coherence of the several formulations between the grids is not generally assumed. Clearly, this point may limit the error reduction in the multigrid algorithm as seen in the FAS correction of the mixture energy balance equation [START_REF] Belliard | Multigrid preconditionning in the genepi software[END_REF].

Primary Fluid Energy Balance Equation

We don't implement the FAS method for the the primary fluid energy balance equation. This equation is very easily and quickly solved on the primary curvilinear grid. However, because it is implied in the source term in the RHS of the mixture energy balance equation, we compute a coarse grid restriction and a coarse grid error of the primary fluid temperature. A fine grid correction is then performed.

✁ ✁ ✂ ✂ ✄ ✁✄ ✄ ✁✄ ☎ ☎ ✆ ✆ ✝ ✝ ✞ ✁✞ ✞ ✁✞ ✟ ✁✟ ✟ ✁✟ ✠ ✁✠ ✠ ✁✠ ✡ ✁✡ ✡ ✁✡ ☛ ✁☛ ☛ ✁☛ ☞ ☞ ✌ ✁✌ ✌ ✁✌ ✍ ✍ ✎ ✁✎ ✎ ✁✎ ✏ ✏ ✑ ✑ ✒ ✒ ✓ ✁✓ ✓ ✁✓ ✔ ✁✔ ✔ ✁✔ ✕ ✕ ✖ ✖ ✗ ✁✗ ✗ ✁✗ ✘ ✁✘ ✘ ✁✘ ✙ ✁✙ ✙ ✁✙ ✚ ✁✚ ✚ ✁✚ ✛ ✁✛ ✛ ✁✛ ✜ ✁✜ ✜ ✁✜ ✢ ✁✢ ✢ ✁✢ ✣ ✣ ✤ ✤ ✥ ✥

AN IMPLEMENTATION

In this Section, we review the expressions of the coarse grid correction terms (term S I l of Equation ( 6)) for the energy balance equation, the momentum balance equation and the mass balance equation, the definition of the error criterion driving the dynamic multigrid cycle and we present our FAS algorithm.

We can do a remark. Multigrid method are first of all useful in case of a big number of cells. To day, a computation involving 250,000 cells needs about 1 Gigabytes memory and three grids is enough even for this reference computation. In fact, with hexahedral elements in 3-D geometry, the coarsening ratio between two consecutive grids is 8 and, for a four level multigrid algorithm, the cell number ratio between the finest and the coarsest grid is 8 4 ✆ 4096. Hence, the coarsest grid for our reference computation contains only about 60 cells. It is too small for an efficient simulation.

Coarse Grid Correction Terms

We have to compute the S I l term of Equ. ( 6) . It is formed by two parts. The first one is the

☎ T l ✄ Xl ☎ φ I l ✄
x☎ dv term built using the restricted variables X (see below). The second one is completely defined by the chosen restriction (nodal weighted average restriction for the nodal residuals). It is worth to notice that the first part of the correction terms are the coarse grid non-linear residuals of the fully non-linear steady-state flow balance equations, built with the restrictions of the primary fluid and mixture variables. The values of the coefficients implied in the coarse grid correction terms (denoted restricted coefficients) are evaluated using these restricted variables. The restricted coefficients are: ¯x l , Ll , ρl , μ T , χT , Λl , ✁ v Rl .

Energy Balance Equation

The term

☎ T l ✄ Xl ☎ φ I l ✄
x☎ dv is computed by (subscript l stands the coarse grid) :

✁ T l ✄ Xl ☎ φ I l ✄ x☎ dv ✁ dv φ I l ✄ x☎ ✁ β l ✄ ¯ ✁ G l ✂ ✁ ∇ ☎ Hl ✞ div ✄ β l ¯x l ✄ 1 ¯x l ☎ ρl Ll ✁ v Rl ☎ ✄✂ ✁ dv φ I l ✄ x☎ β l Ql ✞ ✁ dv ✁ ∇ φ I l ✄ x☎ ✄ β l χT ✁ ∇ Hl ☎ (13)
with the restricted coefficients defined as previously mentioned and Ql the restricted source term, including the restricted volume thermal source computed with the restricted primary fluid temperature, the coarse grid boundary thermal flux and, eventually, the restricted pressure gradient terms (if it exists, this last term is a function of the restricted variables).

Momentum Balance Equation

For the momentum balance equation, the term

☎ T l ✄ Xl ☎ φ I l ✄ x☎ dv is computed by : ✁ T l ✄ Xl ☎ φ I l ✄ x☎ dv ✁ dv φ I l ✄ x☎ ✁ β l ρl ✄ ✁ ¯v l ✂ ✁ ∇ ☎ ✁ ¯v l ✞ div ✄ β l ¯x l ✄ 1 ¯x l ☎ ρl ✁ v Rl ✟ ✁ v Rl ☎ ✄✂ ✁ dv φ I l ✄ x☎ β l ρl ✄ ✁ g Λl ✁ ¯v l ☎ ✁ dv ✁ ∇ φ I l ✄ x☎ ✁ β l Pl β l μ T ✄ ✁ ∇ ✁ ¯v l ✞ ✁ ∇ t ✁ ¯v l ☎ ✄✂ ✞ BC (14)

Mass Balance Equation

The term :

☎ T l ✄ Xl ☎ ψ E l ✄ x☎ dv (a correction field by element) is computed by (E :✆ e l ) : ✁ T l ✄ Xl ☎ ψ E l ✄ x☎ dv ✁ dv ψ E l ✄ x☎ ✁ ∇ ✂ ✄ β l ✁ Ḡl ☎ ✂ (15)
In fact, we don't solve directly this equation, but we enforce the divergence free projection of the gap mass flux:

✄ ✁ G l ✁ Ḡl ☎ .

A Criterion To Stop The Dynamic Multigrid Cycle

An important feature is related to the going out of the multigrid cycles. In fact, the high efficiency of the multigrid solver is drastically reduced after some cycles, leading to a stalled regime in case of static cycle. To overcome this drawback, we need an ad hoc test on the relative variations of the errors (testing if the error no longer decreases). Doing this, we face to dynamic multigrid cycles. The goal is to stop the coarse grid runs (or reduce the associated computational cost) and to disconnect the corrections of the variables for the fine grid. The first point is motivated by the high cost of the computation in regard of the communication one.

For the mixture specific enthalpy and the mixture mass flux, we build the following indicator set (Ω l , l = 1, ... l max ):

(16) ind l ✄ e l ☎ ✆ abs ✄ ✁ e l m L2 e l m✝ 1 L2 ☎ e l re f L2 ☎
where e l is the coarse grid error, m is the multigrid cycle counter,

✂ L2 denotes the discrete L 2 -norm, e l re f L2 a reference L 2 -norm (here, e l 1 L2
) and abs(...) the absolute value function. Doing this, the corrections on the finer grid Ω l ✝ 1 is monitored by the stalled regime detections of the computed errors itself.

For each coarse grid task, we can give a general tolerance level ε MG l . The default value for the industrial SG simulations is 5 ✄ 10 ✝ 3 but it can be user managed. If during the computation, we reach the following condition for one of the coarse grid error (specific enthalpy error or mass flux error):

(17) ind l ✄ e l ☎ ✄ ε MG l ✁
the value of the computed error is bring to zero. When, this condition is reached for the specific enthalpy and the mass flux, the value of the coupling period for the task related to the grid M l is gradually reduced to two iterations (we successively multiply by 0.8 the coupling period).

The Two-grid Algorithm

A GENEPI two-grid cycle is a pseudo-FMG FAS two-grid slash-cycle (without post-smoothing) since no correction is applied during the first cycle. Hence, the effective computation begins by a coarse grid solving. The number of coarse grid pseudo-time steps involved in the first cycle cp 0 l is user managed. If this number is big enough (of the same magnitude that the pseudo-time step number required to solve the coarse grid problem alone), we get a true FMG method. If it is only equal to the current number of pseudo-time steps cp 1 , then we get a pseudo-FMG method.

A Red-black Parallel Version Of A Three-grid Algorithm

A two-grid FAS Method is essentially a sequential method. We have tested a parallel red-black three-grid method. As in a red-black Jacobi method, two groups of grids are set-up (here : {M 0 , M 2 } and {M 1 }). The two groups work sequentially and all the tasks of the same group work in parallel. The groups are set-up in the way that two consecutive grids do not belong to the same group. Following this rule, we apply a sequential non-ideal two-grid method to each couple of consecutive grids. See [START_REF] Belliard | Multigrid preconditionning in the genepi software[END_REF] for more details.

To test our implementation of the FAS method on steam generator two-phase fluid simulations, we present some sequential two-grid and parallel red-black three-grid computations of the CEA Clotaire mock-up [START_REF] Campan | Steam generator experiment for advanced computer code qualification : CLOTAIRE international program[END_REF]. The riser part forms a half cylinder of 0.62 m in diameter and 9.16 m in height. The inside is filled with U-shaped tube bundle, 7.2 m in height, into which the hot primary flow enters. One flow distribution baffle, nine tube support plates and one anti-vibration bar are fixed in respectively the bottom, upright and curved part of the bundle. The simulation fluid is Freon (r114). Except when it is mentioned, the boundary conditions, the physical and numerical parameters are identical for all the computations. We stop the computation when each variable u (specific enthalpy, mass flux, pressure, primary fluid temperature) has verified a steady-state flow criterion crit:

(18) u n ✄ 1 u n L2 u n L2 δt ☎ crit ✂
with crit = 10 ✝ 3 s ✝ 1 or 10 ✝ 4 s ✝ 1 . These computations differ by the number of grids (2 or 3), the number of fine grid cells (22,400 or 88,704), the smoother type (CG or CGS), the pseudo-time step number during the first multigrid cycle (cp 0 l ✁ 0 ✄ l ☎ l max ), the coupling periods cp l and the dynamic multigrid cycle cut-off criteria ε MG l (see Equation 16).

Our reference computation is a pseudo-FMG FAS simulation involving the following numerical parameters:

cp 0 ✆ 15, cp 1 ✆ 60 and cp 2 ✆ 120, preconditioned CG smoother, α ✆ 0✂ 7, ε MG 1 ✆ 10 ✝ 3 and ε MG 2 ✆ 10 ✝ 3 (

dynamic multigrid cycle cut-off criteria on the coarse grids).

Each fine grid is built by subdivision of the coarser one: each coarse grid cell is cut in eight parts. The 22,400 cells fine grid is shown on Figures 2 and3 with its associated coarse grid (2,800 cells). Computations involving 22,400 grid cells are run on a 1,700 MHz PIV Personal Computer (roughly one hour CPU time and 100 Mb). Computations involving 88,704 grid cells are run on the CEA supercomputer IXIA, a 64 Dec-Alpha ES40 stations with four 833 MHz EV68 processors (roughly height hours CPU time and 700 Mb). Results are compared in terms of pseudo-time step number and CPU time (eventually, elapsed time).

Clotaire Mock-up Sequential Two-grid FAS Simulations

The 22,400-cell Grid

Table 1 shows the reference computation results in comparison with the standard computation and the FMG FAS method ones. Concerning the fine grid, Figure 4 shows the convergence histories of the mixture flow variables and Figure 5 shows the evolutions of the corrections of these variables u

( ✁ u new ✝ u old ✁ ✁ u old ✁ L2
), and of the non-linear balance equation residuals in relative discrete L 2 norm. The CPU time speed-up is 1.5 for the reference computation. As a whole, the coarse grid solution does bring an acceleration of the fine grid computation: about 57% of the standard computation fine grid pseudo-time steps are saved (the iteration speed-up is about 1.7). The dynamic multigrid cycles are stopped before the obtention of the 10 ✝ 3 s ✝ 1 steady-state criteria and after this point the standard method is run. The number of multigrid cycles is 25 ( 15*25 = 375 fine grid pseudo-time steps or 1500 coarse grid pseudo-time Figure 4: Clotaire mock-up simulations: comparison of the convergence histories of the fine grid variables for the FAS method and the standard solver in L 2 norm. Pseudo-FMG FAS method: 22,400 cells,

cp 0 ✆ 15, cp 1 ✆ 60, ε MG 1 ✆ 10 ✝ 3 , α ✆ 0✂ 7; CG smoother.
The true application of the FMG algorithm does not increase the CPU time speed-up: the number of fine grid pseudo-time steps roughly is the same as for the pseudo-FMG FAS computation and the first 900 coarse grid pseudo-time steps lead to a CPU time overhead without compensation, see Table 1.

An important feature is the fine grid correction of the primary flow temperature. It provides a substantial increase in the mixture enthalpy convergence. Without correction, the specific enthalpy conver-Table 1: Clotaire mock-up simulations (steady-state: 10 ✝ 3 s ✝ 1 ); sequential two-grid FAS method: 22,400 cells,

cp 0 ✆ 15, cp 1 ✆ 60 (first cp 1 ✆ 900; FMG FAS) α ✆ 0✂ 7, ε MG 1 ✆ 10 ✝ 3 , CG smoother. Std Pseudo-FMG FAS FMG FAS Ω 0 Ω 1 Ω 0 Ω 1 Ω 0 Steady-
cp 0 ✆ 15, cp 1 ✆ 60, ε MG 1 ✆ 10 ✝ 3 , α ✆ 0✂ 7; CG smoother.
gence is not so fast, due to the coupling of the mixture energy balance equation with the primary fluid one and the physical time needed to propagate the primary flow information.

As a whole, the speed-up bring by the pseudo-FMG FAS method is similar to the nested iteration method one for this degres of freedom number. In this very much simpler method, we only use a coarse gird computation result to initialize the fine grid computation. In fact, the 88,7040-cell grid test case will show the advantage of the pseudo-FMG FAS method in the case of a big degres of freedom number.

The 88,7040-cell Grid

Again, the reference computation is confronted to the standard method computation. We give the CPU time and the elapsed time, more representative of the performances seen by the user. For the standard method computation the CPU time and the elapsed one are similar (no communication time). For the FAS methods, speed-up measurements based on the elapsed time lead to low bound values of the CPU time speed-up.

Table 2 gives a comparison between the results and the Figure 6 shows the convergences of the fine grid variables. For this high space resolution computation, the pseudo-FAS method performance is higher than the nested iteration method one (speed-up values about 2.0). As a whole, the FAS algorithm performance is very good, even larger than in the case of the 22,400-cell grid: about 2.5 instead of 1.5. Explanation is to be found in a relative cheaper coarse grid pseudo-time step cost and a lower efficiency of the standard method in term of convergence: 3,887 instead of 1,472 fine grid pseudo-time steps. After 25 multigrid cycles ( 375 fine grid pseudo-time steps), the dynamic multigrid cycle cut-off criterion Table 2: Clotaire mock-up simulations (steady-state: 10 ✝ 3 s ✝ 1 ); sequential two-grid pseudo-FMG FAS method: 88,704 cells, Figure 6: Clotaire mock-up simulations: comparison of the convergence histories of the fine grid mass flux and pressure for the FAS method, the nested iteration method and the standard solver in L 2 norm. Pseudo-FMG FAS methods: 88,704 cells,

cp 0 ✆ 15, cp 1 ✆ 60, α ✆ 0✂ 7, ε MG 1 ✆ 10 ✝ 3 , CG smoother. Standard Pseudo-FMG FAS Ω 0 Ω 1 Ω 0 Pseudo-
cp 0 ✆ 15, cp 1 ✆ 60, cp 2 ✆ 120, ε MG 1 ✆ 10 ✝ 3 , ε MG 2 ✆ 10 ✝ 3 , α ✆ 0✂ 7; CG smoother. (ε MG 1 ✆ 10 ✝ 3
) is reached for the mass flux. And after 35 multigrid cycles ( 525 fine grid pseudo-time steps), the dynamic multigrid cycle cut-off criterion is reached for the specific enthalpy. Then, the FAS method is stopped and the standard one is run. The slope change of the specific enthalpy convergence at this time is related to this solver change, see Figure 6.

Using a higher value of the ε MG 1 criterion leads to reduce the number of fine grid pseudo-time steps. With ε MG 1 ✆ 10 ✝ 4 , 50 multigrid cycles are performed and the fine grid pseudo-time step number is reduced to 1,080 (instead of 1,170). But the increase of the coarse grid pseudo-time step number limits the increase of the speed-up.

Clotaire Mock-up Parallel Red-black Three-grid Pseudo-FMG FAS Simulations

The simulation of the Clotaire mock-up using the 88,704-cell grid is now performed with the parallel red-black three-grid pseudo-FMG FAS algorithm introduced above. The computation parameters are similar to the Section 4.1.2 ones. We confront our reference computation with the sequential two-grid pseudo-FMG FAS method computation of Section 4.1.2.

Table 3 shows the reference computation results. Comparison with the Table 2 points out a fine grid pseudo-time step number reduction: 1,005 instead of 1,170, and an other one for the intermediate grid pseudo-time steps: 1,773 instead of 2,435. This enhanced convergence is summarized in the pseudotime step number speed-up: 3.1 instead of 2.7 for the sequential two-grid pseudo-FMG FAS method computation. Moreover, taking advantage of a parallel run (Ω 0 in parallel with Ω 2 ) and of the Ω 1 pseudo-time step number reduction, the elapsed time speed-up is close to the pseudo-time step speed-up (or the CPU time speed-up).

Table 3: Clotaire mock-up simulation (steady-state: 10 ✝ 3 s ✝ 1 ); parallel red-black three-grid pseudo-FMG FAS method: 88,704 cells, 3.0 Speed-Up (CPU time)

cp 0 ✆ 15, cp 1 ✆ 60 cp 2 ✆ 120, α ✆ 0✂ 7, ε MG 1 ✆ 10 ✝ 3 , ε MG 2 ✆ 10 ✝ 3 , CG smoother. Standard Pseudo-FMG FAS Ω 0 Ω 2 Ω 1 Ω 0 Pseudo-
3.1 Speed-Up (pseudo-time steps)

3.1

Concerning the dynamic multigrid cycle algorithm, the cut-off criterion (ε MG 1 ✆ 10 ✝ 3 ) is reached for the mass flux after 23 multigrid cycles and after 25 multigrid cycles for the specific enthalpy (more quickly than for the two-grid pseudo-FMG FAS case).

Choice Of The Fine Grid Correction Relaxation

Relaxation is a crucial point for the Clotaire mock-up simulation with the FAS method. Figure 7 shows the influence of the choice of the relaxation parameter α for computations involving the 22,400cell grid (two-grid method). We look at the mixture mass flux convergence because generally it is the variable that has the slowest convergence. Clearly, choosing α in the range [0.4; 0.9] leads to similar mass flux convergence histories (even if [0.7; 0.9] seems the best interval). Using a value of 1 (i. e. no relaxation) induces the divergence of the pseudo-FMG FAS method.

Choice Of The Coupling Periods

Figure 8 shows the influence of the choice of the coupling periods for computations involving the 22,400-cell grid (two-grid method). The dynamic multigrid cycle cut-off criterion is fixed to 10 ✝ 5 . This very low value allows to catch all the convergence behaviors induced by the multigrid corrections. Let us notice that the stalled regime is reached after about 1,000 fine grid pseudo-time steps (except for the (15; 60) computation, and after this point, a standard method is run.

As a whole, all the convergence histories are spread near the coarse grid standard method convergence history. But, particular choices of the coupling periods bring the best convergences. It is the case of the couple (15; 60) (understand cp 0 ✆ 15 and cp 1 ✆ 60). Moreover, some similarities in the convergence histories can be found. For instance, (15; 30) and (30; 60) lead to the same convergence. For these two couples, the coarse grid / fine grid ratio is the same. Similar convergence properties can also be found for the couples (15; 60) and (15; 120). This last point underlines that 60 coarse grid pseudotime steps is enough to solve the coarse grid problem. At the opposite, if the two coupling periods are too close each other, as (20; 20) or (15; 30), the coarse grid problem is not solve with enough accuracy, leading to a slower convergence on the fine grid.

Pseudo-FMG FAS Method Scalability

We can define the scalability as the following. Suppose a problem involving nd degrees of freedom solves in nt pseudo-time steps with a multigrid method using a given number of 3D grids. If we multiply nd by height and use one additional grid, then the number of pseudo-time steps should be nt again. Figure 9 shows the comparison of fine grid mass flux convergences for 88,704 cells (three grids), 22,400 cells (two grids) and 2,800 cells (one grid; 22,400/8 cells). Also it shows the convergences related to the one-grid standard method for 88,704 cells, 22,400 cells and 11,088 cells (88,704/8). During the first 500 pseudo-time steps (active multigrid solver), it indicates a pretty nice scalability of the pseudo-FMG FAS method. The cell number ratio between 88,704 and 22,400 is only about four and not eight, but standard method computations show similar mass flux convergence behaviors for the 22,400-cell grid and the 11,088-cell grid test cases (about the same space dicretization following the mean flow direction) that suggests roughly similar convergences for the 11,088-cell and 22,400-cell two-grid pseudo-FMG FAS method.

CONCLUDING REMARKS

The FAS multigrid method has been successfully implemented and tested in the GENEPI code. The high efficiency of this scheme as been proved in the case of an industrial simulation as the Clotaire Benchmark one. The FAS algorithm performances are very good, but need a large amount of computational cells to really be efficient. For a 22,400-cell test case, the CPU time speed-up is relatively low: roughly 1.5 (two grids). But for the 88,704-cell test case, high CPU time speed-up is obtained: about 2.5 (two grids) and 3.0 (three grids). Larger is the number of grid cells (here, one hundred thousand), bigger is the speed-up. Moreover, parameter studies has been performed providing the determination of a set of parameter values leading to the best results. In term of CPU time, the true version of the FMG method is not the best one to maximize the speed-up. Pseudo-FMG algorithms reach a better score. The relaxation of the fine grid error correction and a dynamic manadgement of the multigrid cycles are crucial points to get stable and fast convergences and to save CPU time. Concerning the coupling periods, a good choice is: cp 0 ✆ 15, cp 1 ✆ 60 and cp 2 ✆ 120 (if any).

NOMENCLATURE

Figure

  Figure 1: Forcing Terms

Figure 2 :

 2 Figure 2: Fine (22,400 cells) and coarse grid (2,800 cells).

Figure 3 :

 3 Figure 3: Fine (22,400 cells) and coarse grid (2,800 cells): cross-sections.

Figure 5 :

 5 Figure5: Clotaire mock-up simulations: convergence histories of the fine grid variable corrections (αP 1 0 e 1 ) and non-linear residuals in L 2 norm. Pseudo-FMG FAS method: 22,400 cells, cp 0 ✆ 15, cp 1 ✆ 60, ε MG

Figure 7 :

 7 Figure 7: Clotaire mock-up simulations: comparison of the convergence histories of the fine grid mass flux (L 2 norm) function of the relaxation parameter α. Pseudo-FMG FAS method: 22,400 cells, cp 0 ✆ 15, cp 1 ✆ 60, ε MG 1 ✆ 10 ✝ 3 ; CG smoother.

Figure 8 :

 8 Figure 8: Clotaire mock-up simulations: comparison of the convergence histories of the fine grid mass flux (L 2 norm) function of the coupling periods cp l . Pseudo-FMG FAS method: 22,400 cells, ε MG 1 ✆

Figure 9 :

 9 Figure 9: Clotaire mock-up simulations: comparison of the convergence histories of the fine grid mass flux (L 2 norm). Pseudo-FMG FAS method:cp 0 ✆ 15, cp 1 ✆ 60, cp 2 ✆ 1 ✁ 200, α ✆ 0✂ 7, ε MG 1 ✆ 10 ✝ 3 ; CG smoother. Scalability test.

  cp l : coupling period for grid Ω l cp 0 l : first coupling period for grid Ω l ✁ G : mixture mass flux (= ρ ✁ v ) ✁ g : gravity (m s ✝ 2 ) H : mixture specific enthalpy (J kg ✝ 1 ) H ls : saturated liquid specific enthalpy (J kg ✝ 1 ) l max : maximal number of computational grids L : latent heat (J kg ✝ 1 ) L v : typical vortex length (m) m : multigrid cycle counter n : iteration (or pseudo-time step) counter P : pressure (Pa) P k l : grid Ω k to grid Ω l transfer operator Q : heat source (W m ✝ 3 ) S l : Coarse grid (Ω l ) balance equation RHS (FAS) t : time (s) ✁ v : mixture velocity (m s ✝ 1 ) ✁ v R : relative velocity (gas minus liquid, m s ✝ 1 ) x : static quality ( H ✝ H ls L ) α : multigrid relaxation parameter β : porosity (:= ω m ω) δt : Pseudo-time step (s) χ T : turbulent diffusion coefficient for the mixture energy equation (kg m ✝ 1 s ✝ 1 ) µ T : two-phase turbulent dynamic viscosity (kg m ✝ 1 s ✝ 1 ) ρ : mixture density (kg m ✝ 3 ) Λ : two-phase friction tensor (s ✝ 1 ) ω : volume of the homogenization cell (m 3 ) ω m : mixture volume of the homogenization cell (m 3 ) Ω l : computation domain grids (0 ☎ l ☎ l max ) φ i : nodal function ψ e : element function


Isabelle Ramière
email: isabelle.ramiere@cea.fr

Michel Belliard
email: michel.belliard]@cea.fr

Philippe Angot
email: angot]@cmi.univ-mrs.fr

ON THE SIMULATION OF NUCLEAR POWER PLANT COMPONENTS USING A FICTITIOUS DOMAIN APPROACH

Keywords: Fictitious domain method, Immersed interface, Neptune, Steam Generator

This paper describes the first step of the introduction of fictitious domain methods (FDM) in the simulation of nuclear component fields. We consider the resolution of a diffusion problem ✁ div ✂ ã✄ ✆☎ ũ✝ ✟✞ f , in a physical domain ✠ , with different kinds of boundary conditions (B.C.). The shape of this domain is of the Steam Generator's type. The aim of this work is to present several ways to impose B.C. on an immersed interface, without modifying locally the numerical scheme and without using Lagrange multipliers. The numerical resolution is implemented using either a finite element or a finite volume method. Error estimates are performed in order to evaluate the capability of this method. A comparison with unstructured mesh methods is also presented.

INTRODUCTION

Within the framework of the federal program NEPTUNE launched by the CEA and EDF concerning the industrial simulation of two-phase flows in the nuclear power plant components as steam generators (SG) or cores, we propose investigations of new numerical methods based on structured meshes. The advantages of numerical methods for partial differential equations (PDE) involving structured meshes are well known : easy implementation for fast solvers (based for instance on finite volume methods with structured grids) and for multi-scale methods [START_REF] Angot | A comparison of locally adaptive multigrid methods[END_REF], good convergence properties, easy simulation of moving boundaries (for example moving obstacles in a flow, fluid/structure interaction...) ... However their main drawback is the difficulty to take into account boundaries or obstacles in complex industrial situations. As a fact, numerical methods involving unstructured meshes are classical in the industrial software world.

Fictitious domain method

Since a few years, fictitious domain methods (FDM) (Saul'ev, 1963. (in Russian), [START_REF] Marchuk | Methods of Numerical Mathematics[END_REF]Marchuk, (1rst ed. 1975)) have been arising for Computational Fluid Dynamic [START_REF] Khadra | Fictitious domain approach for numerical modelling of Navier-Stokes equations[END_REF]. In this approach, the initial computation domain ✠ (called physical domain) is immersed in a geometrical bigger and simpler other one ✠ , called fictitious domain (see Figure 1). The spatial discretization is now performed in the fictitious domain The main difficulty lies in the choice of the PDE solved in the fictitious domain ✠ and in the numerical scheme used for the resolution. These two choices have to be coupled in order to take into account the conditions on the initial boundary ✁ ✠ . Therefore the restriction of the solution on the physical domain is expected to be the initial problem's one.

Numerically, there are mainly three approaches to deal with the boundary conditions located on the immersed boundary: ✂ "Thin" interface approaches : the physical boundary is approximated without being enlarged in the normal direction. The physical boundary and the approximated one lie in the same ✄ n space. For example, in this group, we can find truncated domains methods [START_REF] Mccorquodale | A cartesian grid embedded boundary method for the heat equation on irregular domains[END_REF][START_REF] Ye | An accurate cartesian grid method for viscous incompressible flows with complex immersed boundaries[END_REF], immersed interface methods (I.I.M.) [START_REF] Leveque | The immersed interface method for elliptic equations with discontinuous coefficients and singular sources[END_REF], as well as an adapted Galerkin method proposed by Latché (submitted (october 2003)).

✂ "Spread" interface approaches : the support of the approximated interface is larger than the physical one. The approximated interface has one dimension more than the physical one. For example, the boundary can be estimated as the set of the elements of discretization crossed by itself. In particular, this kind of approach is used in the interface boundary methods (I.B.M.) [START_REF] Cortez | The blob projection method for immersed boundary problems[END_REF].

✂

Lagrange multipliers : in FDM using Lagrange multipliers [START_REF] Glowinski | On the solution of the Dirichlet problem for linear elliptic operators by a distributed Lagrange multiplier method[END_REF], the initial B.C. are imposed weakly as side constraints. This approach doesn't use a continuous approximated boundary but only needs discretization points lying on the initial boundary.

Moreover, penalization methods [START_REF] Angot | Finite volume methods for non smooth solution of diffusion models; Application to imperfect contact problems[END_REF][START_REF] Khadra | Fictitious domain approach for numerical modelling of Navier-Stokes equations[END_REF] have proved their efficiency to impose values on an immersed interface or to solve suitable equations in different domains.

In this paper, one is interested in a spread interface approach computed with a usual F.E. Scheme as well as a thin interface approach computed using a F.V. Scheme with flux and solution jumps [START_REF] Angot | A model of fracture for elliptic problems with flux and solution jumps[END_REF] are exposed.

A lot of papers are dedicated to embedded Dirichlet or Neumann B.C., e.g. (Saul'ev, 1963. (in Russian), [START_REF] Marchuk | Methods of Numerical Mathematics[END_REF]Marchuk, (1rst ed. 1975)), [START_REF] Glowinski | On the solution of the Dirichlet problem for linear elliptic operators by a distributed Lagrange multiplier method[END_REF][START_REF] Glowinski | Wavelet and finite element solutions for the Neumann problem using fictitious domains[END_REF][START_REF] Leveque | The immersed interface method for elliptic equations with discontinuous coefficients and singular sources[END_REF] and the reference herein... Only a few studies are devoted to other embedded Fourier B.C. (Kopčenov, 1974. (in Russian), [START_REF] Marchuk | Methods of Numerical Mathematics[END_REF](1rst ed. 1975[START_REF] Angot | Finite volume methods for non smooth solution of diffusion models; Application to imperfect contact problems[END_REF] ✝ isn't a saddle-point problem, the inf-sup condition must not be verified [START_REF] Glowinski | On the solution of the Dirichlet problem for linear elliptic operators by a distributed Lagrange multiplier method[END_REF]. Moreover, only one discretization grid is used (the structured regular mesh over the fictitious domain).

Presentation of the test problem

For sake of simplicity we choose to focus on 2D problems, but the formulations can be extended to 3D problems without more difficulty. We study the resolution of a diffusion problem in the unit disk ✠ . For symmetry reasons, we can solve this problem only on a quarter of this disk.

Let us consider the following model problem: ✠✄✂

✄ 2 open-bounded domain For ã ☎ ✂ L✆ ✂ ✠ ✝ ✝ 2 ✝ 2 and f ☎ L 2 ✂ ✠ ✝ , find ũ defined on ✠ such that ✞ ✟ ✡✠ ✁ div ✂ ã ✄ ✆☎ ũ✝ ✞ f in ✠ ✂ ã ✄ ✆☎ ũ✝ ✄ n ✞ 0 on ☛ B.C. on ☞ (1)
where B.C. represents :

✂ a Dirichlet condition : ũ ✞ u D , u D ☎ H 1✌ 2 ✂ ☞ ✝ z x 1 1 Ω Γ ~Σ ✂ a Robin (or Fourier) condition : ✁ ✂ ã ✄ ✆☎ ũ✝ ✄ n ✞ ✎✍ R ✂ ũ ✁ u R ✝ ✑✏ g R , ✍ R ☎ L✆ ✂ ☞ ✝ ✓✒ ✔✍ R ✕ 0✖ u R and g R ☎ L 2 ✂ ☞ ✝
(with n the outward normal unit vector on the circle) Remark: a Neumann condition, ✁ ✂ ã ✄ ✆☎ ũ✝ ✄ n ✞ g, is considerate as a particular Robin condition where ✍ R ✗ 0 and g R ✗ g.

Moreover, the symmetric tensor of diffusion ã ✗ ✂ ãij ✝ 1✘ i j✘ 2 verifies the classical ellipticity assumption:

✙ ã0 ✚ 0✖ ✜✛ ✣✢ ✤☎ ✄ 2 ✖ ã✂ x✝ ✄ ✡✢ ✄ ✡✢ ✕ ã0 ✥ ✢ ✥ 2 a✄ e✄ in ✠ ✂ A1✝
where ✥ ✄ ✥ is the Euclidean norm in ✄ 2 . In that case, classical variational techniques [START_REF] Dautray | Analyse mathématique et calcul numérique pour les sciences et les techniques[END_REF] prove the existence and the uniqueness of the solution ũ ☎ H 1 ✂ ✠ ✝ of this problem.

For sake of simplicity, we suppose here that u D and u R are constant. The non constant cases can be treated without more difficulty (e.g. with cell-centered approximation of these variables).

In a fictitious domain approach, the quarter disk is immersed in the unit square In the next section, we present in detail the numerical methods implemented on the fictitious domain ✠ . For each method, the equation solved on the fictitious domain is introduced. We expose several ways to impose the different B.C. on an immersed interface for a diffusion problem. Then, in Section 3, the accuracy of this fictitious domain approach will be illustrated with numerical results. To conclude, future works and perspectives are pointed out.

✠ ✞ ✧✦ 0✖ 1★ ✪✩ ✫✦ 0✖ 1★ . z x 1 1 Ω Γ Γ ∼ ∼ Σ Ω e Ω e

NUMERICAL METHODS

Finite Element (F.E.) Method

Numerical scheme

The F.E. numerical resolution is implemented using a 1 F.E. Scheme, where k stands for the space of polynomials of degree of each variable less than or equal to k. For example, 1 ✞ span✁ 1✖ x✖ y✖ xy✂ in ✄ 2 . The 1 discretization nodes are located on the vertex of the discrete elements.

A "spread" interface approach has been chosen. The set of the discretization elements K crossed by the immersed interface ☞ forms the approximated interface ✄ h ☎ .

Discretization of the fictitious domain

✠ : ✠ h , ✄ h ☎ and ✠ e h
In comparison with the initial formulation, the equation solved in

Find u h ☎ V h ✞ ✝✁ v h ☎ C 0 ✂ ✠ ✝ ✓✖ v h ✥ K ☎ ✞ 1 ✛ K ☎ ✟✆ h ✂ ✂ H 1 ✂ ✠ ✝ such that ✞ ✟ ✠ ✁ div ✂ a ✄ ✆☎ u h ✝ ✑✏ b ✂ u h ✁ u E ✝ ✞ f ✏ ✠ in ✠ h ✂ a ✄ ✆☎ u h ✝ ✄ n ✞ 0 on ☛ suitable B.C. for u h on ☛ e (2)
where a ☎

✂ L✆ ✂ ✠ h ✝ ✝ 2 ✝ 2 (symmetric tensor), f ✖ ✡✠ ✖ u E ☎ L 2 ✂ ✠ h ✝ , and b ☎ ✂ L✆ ✂ ✠ h ✝ ✝ such that: a ✥ ☛ h ✞ ã ✥ ☛ h ✖ f ✥ ☛ h ✞ f ✥ ☛ h
a satisfies the assumption (A1) in ✠ h and b the following one: [START_REF] Dautray | Analyse mathématique et calcul numérique pour les sciences et les techniques[END_REF] allow to conclude that the approximated solution u h of the problem (2) exists and is unique.

✙ b 0 ✕ 0✖ b ✂ x✝ ✕ b 0 a✄ e✄ in ✠ h (A2)
The additional terms b✖ u E and ✠ enable to take into consideration the immersed boundary [START_REF] Khadra | Fictitious domain approach for numerical modelling of Navier-Stokes equations[END_REF], Angot et al. , 1999). We expect thus that u h ✥ ☛ h ✞ ũ ✥ ☛ h . For each kind of boundary conditions lying on the immersed boundary ☞ , different possibilities to impose these conditions without modifying the numerical scheme and without using local unknowns are presented in Section 2.1.2.

Treatment of the initial B.C.

✂

Dirichlet case:

The Dirichlet B.C. are treated by penalization [START_REF] Angot | A penalization method to take into account obstacle in incompressible viscous flows[END_REF]. Setting b ✞ 1 ☞ (where 0 ✌ ✎✍ ✏✌ ✑✌ 1 is called penalty coefficient) and u E ✞ u D enables to impose u h ✞ u D . It's the L 2 -penalty. When the coefficient a is also equal to 1 ☞ (a ✞ 1 ☞ Id), we obtain the H 1 -penalty. The other coefficients of the equation ( 2 on ✠ (Ramière et al. , in preparation). In Section 3, we will compare the penalization of the spread interface to the penalization of the exterior domain.

Concerning

Find u ☎ H 1 ✂ ✠ ✝ such that ✛ v ☎ H 1 ✂ ✠ ✝ ✄ ☛ a✄ ✆☎ u✄ ✆☎ v dV ✞ ✄ ☛ f v dV ✁ ✏✌ ★ ★ ✂ a✄ ✆☎ u✝ ✄ n✦ ✦ ☎ ✆☎ ☎ ✖ v ✚ (3) where ✝ ★ ★ ✂ a✄ ✆☎ u✝ ✄ n✦ ✦ ☎ ✞ ✂ a✄ ✆☎ u✝ ✄ n ✥ ✟✞ ☎ ✁ ✂ a✄ ✆☎ u✝ ✄ n ✥ ✡✠ ☎ ☎ ☎ is a measure carried by the hypersurface ☞ ✖ ☎ ☎ ☞☛ v ✌ ✌ ☎ ☎ ✖ v ✚ ✞ ✎✍ ☎ v ✥ ☎ dS
In the distribution sense, we obtain the following equation:

✁ div ✂ a✄ ✆☎ u✝ ✞ f ✁ ★ ★ ✂ a✄ ✆☎ u✝ ✄ n✦ ✦ ☎ ☎ ☎ (4)
The jump of flux across ☞ can be interpreted as a source term carried by ☞ .

In our case, we want

✁ ✂ a✄ ✆☎ u✝ ✄ n ✥ ✡✠ ☎ ✞ ✄✍ R ✂ u ✁ u R ✝ ✑✏ g R and ✁ ✂ a✄ ✆☎ u✝ ✄ n ✥ ✟✞ ☎ ✞ 0 so that ★ ★ ✂ a✄ ✆☎ u✝ ✄ n✦ ✦ ☎ ✞ ✄✍ R ✂ u ✁ u R ✝ ✑✏ g R On ✠
, then we have:

✁ div ✂ a✄ ✆☎ u✝ ✞ f ✁ ★ ✡✍ R ✂ u ✁ u R ✝ ✑✏ g R ✦ ☎ ☎ (5)
However, with a structured mesh on ✠ , the support of ☞ is not exactly defined. We introduce a characteristic parameter ✏ in order to approximate the measure ☎ ☎ supported by ☞ by mollifiers [START_REF] Brezis | Analyse fonctionnelle, Théorie et applications[END_REF] on the spread interface ✄ h ☎ .

The principle is the following:

✄ ☎ ✁ ✍ R ✂ u ✁ u R ✝ ✏ g R ✂ v dS ✞ ✄ ✒✑ h✓ ✔ ✍ R ✂ u ✁ u R ✝ ✏ g R ✏ v dV ✛ v ☎ V h (6)
The equation solved on

✠ h is: ✁ div ✂ a✄ ✆☎ u h ✝ ✞ f ✁ ★ ✍ R ✏ ✂ u h ✁ u R ✝ ✏ g R ✏ ✦ ☎ ✑ h✓ ✔ (7)
where ☎ ✑ h✓ ✔ is a measure carried by the spread interface ✄ h ☎ .

With the equation ( 7), the parameter of the generic formulation (2) can be easily set.

✝ ✂ a✄ ✆☎ u h ✝ ✄ n ✞ 0 on ✄ h ☎ ☛ e u h ✞ u ext otherwise 6/17
The 11 th International Topical Meeting on Nuclear Thermal-Hydraulics (NURETH-11) Popes' Palace Conference Center, Avignon, France, October 2-6, 2005

where u ext is unspecified.

The parameter ✏ can be estimated by several ways (Ramière et al. , in preparation). We present some of them in Section 3. If h is the discretization step, [START_REF] Angot | Contribution à l'étude des transferts thermiques dans des systèmes complexes; Application aux composants électroniques[END_REF] showed that ✏ is in ✂ h✝ .

Local adaptive mesh refinement (A.M.R.)

Using a spread interface approach, A.M.R. techniques are necessary to improve the accuracy of the solution near the immersed interface and by the way on the whole physical domain. Most of these techniques are derived from multi-grid methods [START_REF] Hackbusch | Multi-grid Methods and applications[END_REF].

A Local Defect Correction (L.D.C.) method [START_REF] Hackbusch | Local Defect Correction Method and Domain Decomposition Techniques[END_REF] is performed. This method is a multi-grid method with a defect restriction in the ascent step. The local refinement zone is composed by all the elements of the spread interface ✄ h ☎ and their neighbors. This choice enables to correct the values of all the nodes located on the spread interface. As a consequence, the immersed B.C. are imposed more precisely. As the solution of an elliptic problem strictly depends on the B.C., we expect an improvement of the solution on the physical domain.

Finite Volume (F.V.) Method

Numerical scheme

The F.V. numerical resolution is performed using a model of fracture with flux and solution jumps [START_REF] Angot | A model of fracture for elliptic problems with flux and solution jumps[END_REF].

Given an "admissible" finite volume mesh ✆ h of ✠ [START_REF] Angot | A model of fracture for elliptic problems with flux and solution jumps[END_REF], with cell-centered discretization nodes, let S be the family of edges or sides ✁ of the "control volumes" K ✂

★ ★ ✂ a✄ ✆☎ u h ✝ ✄ n✂ ✦ ✦ ✂ ✞ ✍ ✂ ✂ u h ✥ ✂ ✁ U✂ ✝ ✁ h✂ (8) ✂ a✄ ✆☎ u h ✝ ✄ n✂ ✥ ✂ ✞ ✟ ✂ ★ ★ u h ✦ ✦ ✂ ✁ g✂ (9)
where the transfer coefficients ✍ ✂ ✖ ✝✟ ✂ satisfy the following ellipticity assumptions:

✍ ✠✂ ✖ ✝✟ ✡✂ ☎ L✆ ✂ ✁ ✝ ✓✒ ✙ ✟ 0 ✕ 0 ✍ ☛✂ ✂ x✝ ✓✖ ✝✟ ✡✂ ✂ x✝ ✕ ✟ 0 on ✁ and g✂ ✖ h✂ and U✂ are given in L 2 ✂ ✁ ✝ .
The finite volume scheme built on ✆ h doesn't introduce unknowns located on the side of the volume control. Indeed, the equations ( 8) and ( 9) enable to express these unknowns according to the cell-centered ones. This F.V. scheme is based on a "thin" interface approach.

In this F.V. approximation, we suppose that the thin approximated interface ☞ h lies on sides of control volumes crossing the physical interface ☞ . n☎ h denotes the outward unit normal vector of ☞ h (n☎ h ✄ n☎ ✚ 0). Equations ( 8) and ( 9) enable to impose the immersed B.C. on ☞ h . The problem solved on ✠ as the same generic form as in (2) with the above conditions of jump moreover.

✍ ✂ ✞ h✂ ✞ g✂ ✞ 0 and ✟ ✂ ✌ ✆☎ ✂ ✞✝ ✂ ✞ 1 ✟ ✂ ✞ 0✖ fracture resistance✝ (10)
Remark: Locally, on the physical domain ✠ h and on the exterior domain ✠ e h , we get a F.V. scheme without jump of flux and without jump of the solution as in [START_REF] Eymard | Finite Volume methods[END_REF].

We focus on the transfer coefficients lying on

✁ ☎ ☞ h ✁ ✂ ☛✄✂ ☛ e ✝ .

✂

Dirichlet case:

We present two manners to deal with a Dirichlet condition. The first one consists in penalizing the exterior domain at u D . To impose u D on the interior side of the approximated interface ☞ h (side located in the physical domain), we set

✟ ✂ ✞ 1 ☞ ✛ ✁ ☎ ☞ h ✁ ✂ ☛✟✂ ☛ e ✝ .
The other transfer coefficients are equal to zero. In the second approach, u D is imposed directly by surface penalization on the approximated interface ☞ h using the jump equations ( 8) and (9). We set

✟ ✂ ✞ 1 ☞ in order to have u ✞ ✞ u ✠ on ✁ and then u ✥ ✂ ✞ u ✞ ✞ u ✠ ✞ u ✥ ✂ . In this case, penalizing also ✍ ✂ ✞ 1 ☞ while setting U✂ ✞ u D induces u✂ ✞ u D ✛ ✁ ☎ ☞ h ✁ ✂ ☛✠✂ ☛ e ✝ .
The exterior domain has no influence on the solution obtained in the physical domain.

In these two approaches, the B.C. on ☞ h ☛ e are the only one of interest :

u h ✥ ☎ h ✡ ✂ e ✞ u D . ✂ Robin case: We want that ✁ ✂ a✄ ✆☎ u h ✝ ✄ n☎ h ✥ ✠ ☎ h ✞ ✄✍ R ✂ u ✥ ☎ h ✁ u R ✝ ✑✏ g R , where u ✥ ☎ h ✞ u ✠ ✥ ☎ h .
There are many ways to deal with the transfer coefficients in order to impose a Robin B.C. on ☞ h (there is five coefficients for 2 equations). One of particular interest is the one without control on the exterior domain. In this case, we eliminate ✁ ✂ a✄ ✆☎ u h ✝ ✄ n☎ h ✥ ✞☎ h and u ✞ ✥ ☎ h in the equations ( 8) and ( 9):

✛ ✁ ☎ ☞ h ✁ ✂ ☛✟✂ ☛ e ✝ ✍ ✂ ✞ 2✍ R ✖ ✟ ✂ ✞ ✍ R 2 ✖ u✂ ✞ u R ✖ h✂ 2 ✁ sg ✂ n☎ h ✄ n✂ ✝ g✂ ✞ g R ✂ e✄ g✄ h✂ ✞ ✁ 2g R and g✂ ✞ 0✝
As in the spread interface approach, a characteristic parameter ✏ is required to respect the equality between the flux carried by the approximated interface and the flux carried by the physical immersed interface. Remark: We can also use a spread interface approach by the same way as exposed Section 2.1.2. In this case, ✛ ✁ ☎ S, the transfer coefficients are set as in (10).

8/17

The 11 th International Topical Meeting on Nuclear Thermal-Hydraulics Popes' Palace Conference Center, Avignon, France, October 2-6, 2005

NUMERICAL RESULTS

To solve the following test problems, the different approaches have been computed on six meshes: 4x4, 8x8, 16x16, 32x32, 64x64 and 128x128 cells. The discretization steps are respectively

h ✞ 1 4 ✖ 1 8 ✖ 1 16 ✖ 1 32 ✖ 1 64 ✖ 1
128 . All these simulations have been computed on a 3.2 GHz Intel Xeon bi-processor.

The linear system can be solved using a direct solver or an iterative one. As the matrix system is symmetric positive definite, an iterative solver of conjugate gradient's type seems to be suitable. Moreover, the penalization coefficients render the matrix ill-conditioned. Thus, a diagonal preconditioner is necessary in order to speed up the resolution of the system. The F.E. simulations have been computed thanks to the finite element industrial code PYGENE [START_REF] Grandotto | Calculs des écoulements diphasiques dans les échangeurs par une méthode aux éléments finis[END_REF][START_REF] Grandotto | A 3D finite element analysis for solving two phase flow problems in PWR steam generators[END_REF] of the Neptune project. This project, codeveloped by the CEA and EDF, is dedicated to the simulation of two-phase flows in Nuclear Power Plants. A F.V. code using a structured approach had been implemented for the scheme [START_REF] Angot | A model of fracture for elliptic problems with flux and solution jumps[END_REF] presented here.

Dirichlet case

The problem solved on the physical domain ✠ is :

✞ ✟ ✠ ✁ ✁ ũ ✞ 4 in ✠ ✂ ˜u ✂ n ✞ 0 on ☛ ũ ✞ 0 on ☞
With the generic formulation (1), we set:

ã ✗ 1 ✂ ã ✗ Id✝ ✓✖ f ✗ 4 and ũ ✞ u D ✞ 0 on ☞ (Dirichlet B.C.)
The analytic solution of this problem is :

ũ ✞ 1 ✁ r 2 in ✠
We compute this problem on the fictitious square domain ✠ using the two approaches described Section 2.1.2 for the F.E. method. For the F.V. method, we test the two approaches exposed Section 2.2.2 for an exterior thin interface, and we also implement the second approach (surface penalization) with a "cut" interface approach.

As in the F.E. method, the discretization nodes are located on the vertex of the discrete elements the H 1 -penalty leads to the same results as the L 2 -penalty as long as the penalization parameter ✍ is small enough. However investigating the behavior of the error according to the penalization parameter (modelization error), we observe that the error obtained for the H 1 -penalty converges faster toward the discretization error as the L 2 one. We set ✍ ✞ 10 ✠ 12 to get a negligible modelization error. In the F.V. method, the discretization nodes are cell-centered. A H 1 -penalty is necessary to impose the Dirichlet value on the approximated interface (located on sides of control volume).

For each approach, Figure 4 represents the relative discrete L 2 error norms versus the discretization step of the fictitious domain.

As excepted, the order of the method is approximatively equal to one for all the approaches since the approximation of the interface is in ✂ h✝ (for the spread as for the thin one). The interface discretization error leads to a "global" discretization error in ✂ h✝ even if the numerical scheme error is in In the F.E. approach, the penalization of the exterior domain is more accurate than the penalization of the spread interface. These results strictly depend on the geometry of the physical domain. Since the discretization nodes are located on the elements vertex, nodes inside the physical domain are penalized with a spread interface penalization. In the quarter disk case, interior penalized nodes are globally farther of the physical interface than the exterior nodes.

In the F.V. approach, we can observe that, with an exterior approximated interface, the H 1 -penalty of the exterior domain and the surface penalization lead to the same errors. Indeed, with the H 1 -penalty, the solution and its gradient is penalized. So u ✞ u D on the exterior domain until the approximated interface. Moreover, with a thin interface approach, the errors obtained are better for a "cut" approximated interface than for an exterior interface. The cut approximation of the interface is more precise. The last conclusion to draw is that the effect of the H 1 exterior penalization is similar with the two schemes.

As drawn on the Figure 5 the main differences between the approximated solution and the analytic one are located around the spread interface. Using the F.E. scheme, an adaptive mesh refinement is performed in this zone. At each level an exterior penalization is performed. A three-grid LDC algorithm (two refinement levels) is applied on each initial mesh. At each level, the local refinement zone is composed by all the elements of the spread interface ✄ h ☎ and their neighbors. This algorithm converges by three V-cycles. The errors computed using a three-grid mesh refinement and the ones obtained without refinement are represented on Figure 6. The results obtained with a local refinement are in ✂ h✝ where h is the discretization step of the finest refinement grid. With A.M.R., on the initial "coarse" mesh, the error obtained after correction is similar to the one obtained without A.M.R. on a mesh with a discretization step equal to the finest local grid's one. Even if the order of the method doesn't increase with a local refinement, the errors on the initial meshes are reduced.

Robin case

We consider the following problem:

✞ ✟ ✠ ✁ ✁ ũ ✞ 16r 2 in ✠ ✂ ˜u ✂ n ✞ 0 on ☛ ✁ ✂ ˜u ✂ n ✞ u ✏ 3 on ☞
Identifying with the generic formulation (1), we get:

ã ✗ 1 ✂ ã ✗ Id✝ ✓✖ f ✞ 16r 2 ✖ ✍ R ✞ 1✖ u R ✞ 0 and g R ✞ 3 ✂ Robin B✄ C✄ ✝
The solution of this problem is : ũ ✞ 2 ✁ r 4 in ✠

In the F.E. approach, the characteristic parameter ✏ can be estimated be several ways. In this section we present three of them. From the equation ( 6), supposing that the flux carried by the spread interface and the flux carried by the physical interface are rather similar (this approximation is as much justified as the discretization step becomes small), we obtain :

✄ ☎ dS ✄ ✑ h ✓ ✔ 1 ✏ dV (11)
The three estimations of ✏ can be deduced of the equation ( 11):

✂ In a first approach, ✏ is supposed constant on ✄ h ☎ : ✏ ✞ meas ✂ ✄ h ☎ ✝ meas ✂ ☞ ✝ ✂
In the second approach, in the equation ( 11), the integration on ✄ h ☎ is weighted by a coefficient ✁ . This coefficient estimates the presence rate of the physical domain in each element K crossed by the boundary ☞ (K By this way, the right hand side of ( 11) is integrated only on the physical domain included in ✄ h ☎ . We get :

✂ ✄ h ☎ ). By construction, ✁ is constant on each K : ✁ K ✞ volume of ✠ included in K volume of
✏ K ✞ e ★ ✁ e ✄ meas ✂ K ✝ ✦ ✁ e ✄ meas ✂ ☞ ✝
The value of ✏ is given by element. This approximation of ✏ is called volume weighted approximation.

✂

In the last approach, the boundary ☞ is piecewise linear approximated by a segment ☞ l K in each element K included in ✄ h ☎ . The equation ( 11) is written for each of these elements K (respect to equation ( 6), we use test functions v lying only on K). Using the linear approximation

☞ l ✞ ✂ K ✁ ✑ h✓ ✔ ☞ l K of ☞ , we get: ✄ ☎ l✓ K dS ✄ K 1 ✏ dV Finally, ✏ K ✞ meas ✂ K✝ meas ✂ ☞ l K ✝
As in the second case, the value of ✏ depends on the element K. This approach induces a local reconstruction of the interface.

For these three estimations of ✏ , the errors obtained according to the methodology described Section 2.1.2 are reported Figure 7(a).

In the F.V. approach, the characteristic parameter is estimated by two manners.

✂

Global correction : the characteristic parameter ✏ is deduced from the same equation as in ( 11) with a surface integral on the thin approximated interface ☞ h instead of the volume integral on the spread interface ✄ h ☎ :

✏ ✟✞ meas ✂ ☞ h ✝ meas ✂ ☞ ✝ ✂
Local correction : A local estimation of ✏ is made in each cell K of the mesh crossed by the immersed interface ☞ . In this case ✏ K corresponds to the ratio between the sum of the measures of the normal projections of the physical boundary ☞ on each mesh axis (x and y in 2D), and the measure of the physical immersed interface ☞ itself. In 2D, with a piecewise linear approximation ☞ l of ☞ composed by a segment ☞ l K in each K crossed by ☞ , we get:

✏ K ✞ cos ✂ K ✏ sin ✂ K
where ✂ K indicates the angle between the normal direction n☎ l✓ K of ☞ l K and the horizontal axis x.

Figure 7(b) shows the errors obtained for these two estimations of ✏ with the transfer parameters introduced Section 2.2.2 for two kinds of thin approximation of the interface (exterior and cut approximated interface).

All the F.E. variants are approximatively of order one (a little bit under). The method using a piecewise linear approximation of the interface leads to little better results but not as excepted. Some improvements are presented Section 3.3. Concerning the F.V. variants (using a thin interface approximation), the global correction leads to an asymptotic stagnation of the error with a cut interface and then the first-order precision is lost. This stagnation disappears with the local correction (see also Figure 8(b)). In Fig. 7(b), for an exterior interface, the two estimations of ✏ seem to lead to the same errors and then to the same first-order method. However, in Fig. 8, we can observe that for finer meshes (256 ✩ 256✖ 512 ✩ 512 cells) a stagnation appears also for a global correction since the local correction keeps the first-order precision. As in the Dirichlet case, the errors obtained with the cut interface are better as the ones obtained with an exterior interface.

For the F.E. approach, a local adaptive mesh refinement is performed on the method involving a constant epsilon (✏ ✞ mes ✁ ✑ h✓ ✔ ✄✂ mes ✁ ☎ ✂ ) and a local epsilon (linear approximation ✏ K ✞ meas ✁ K ✂ meas ✁ ☎ l✓ K ✂ ). As in the Dirichlet case, we computed a three-grid LDC algorithm, which converges by three V-cycles. The results are presented Figure 9.

Here again the refinement method has a discretization error in ✂ h✝ with h the discretization step of the local finest grid. However, performing a constant ✏ , we observe a stagnation of the error for fine meshes and the first-order is lost. In this case, the modelization error is reached. This modelization error seems to be due to the global estimation of ✏ (equation ( 11)) since no stagnation appears with a local estimation of ✏ . The modelization error for a local ✏ decreases with the discretization step.

For the two interface modelling (spread and thin), a local correction is required to keep the first-order method. With a global correction, the rough estimation of ✏ leads to a stagnation of the error.

Some improvements

We present briefly some improvements under investigations for a Robin B.C.. Those improvements are presented for the F.E. approach.

✂ Firstly, we try to correct the equation coefficients on the spread interface. The aim of this correction is to better approximate the physical domain.

A first step is to set ✁ e f instead of f on ✄ h ☎ . In this case the source term is approximatively integrated on the physical domain (the integration is exact for a constant f ). Figure 10(a) presents the errors obtained for the three ✏ introduced Section 3.2 with a source term weighted on the spread interface. Figure 10(b) enables to appreciate the improvement of this approach compared to the "classic" one of Section 3.2. The correction of the source term f leads to better results for each ✏ approach. Moreover, for a constant ✏ , the stagnation to the modelization error is reached for a 32 ✩ 32 mesh, without using A.M.R.

In a second step, for a local ✏ (✏ linear approximated), the diffusion coefficient is also weighted by ✁ e . This approach is implemented in order to get the diffusion coefficient compatible with the This approach has a much better order as the other ones. Indeed, the order reaches almost 2. However, even with these corrections of the source term and of the diffusion coefficient, the discretization error stagnates around a modelization error. With a local ✏ , this modelization error is reduced compared to the constant ✏ 's one.

✂

Secondly, beginning from the equation (3), we computed on the spread interface the surface integral of the jump of flux :

✌ ★ ★ ✂ a✄ ✆☎ u h ✝ ✄ n✦ ✦ ☎ ☎ ☎ ✖ v ✚ ✞ ✄ ☎ ★ ★ ✂ a✄ ✆☎ u h ✝ ✄ n✦ ✦ ☎ v dS ✛ v ☎ V h ✂ H 1 ✂ ✠ ✝
On each element of the spread interface, we choose the base functions associated at each node as test functions. If we call x B the barycenter of the immersed interface intersected in an element, the local surface integral is approximated by an one point Gauss formula on x B . This approach needs a local reconstruction of the interface but avoids the estimation of a characteristic parameter ✏ in order to obtain a volume integral. Following figures (Figure 12(a) and 12(b)) compare this surface integral approach (weighting the source term f and the diffusion coefficient a) with an unstructured mesh approach in term of error/discretization step and error/CPU time. In the approximated surface integral approach, a slope break prevents from determining the order of the method. However, the errors obtained with this approach are comparable to those obtained for unstructured meshes, in term of discretization error as in term of in CPU time. Moreover, in

The 11 th International Topical Meeting on Nuclear Thermal-Hydraulics (NURETH-11) Popes' Palace Conference Center, Avignon, France, October 2-6, 2005 15/17 term of CPU time, the results presented here for the surface integral approach aren't optimal as the software used is dedicated to unstructured meshes.

CONCLUSIONS

The first step of the introduction of a general fictitious domain approach to deal with structured meshes for complex industrial applications has been introduced. The fictitious domain approach introduced here don't use local unknowns and don't modify locally the numerical scheme. Moreover, all the usual embedded B.C. (Dirichlet, Neumann, Fourier) can be treated with these methods.

For a diffusion problem in a domain of SG's type, we exposed : ✂ the methodology to take into account usual B.C. lying on a immersed interface using two numerical schemes : a F.E. scheme and a F.V. scheme. The two main approximated interface approaches have been tested : spread and thin interface. ✂ some recent improvements in order to get better results (comparable to unstructured meshes).

The results presented here give confidence for the use of fictitious domain methods in more complex industrial applications. Even if the discretization error is in ✂ h✝ , the computed relative errors are rather small. Moreover the use of structured meshes is full of promise : meshing costs are weak, local adaptive mesh refinements are easily implemented, moving boundaries can be simulated without mesh reconstruction, fast solvers can be used... The resolution of an advection-diffusion problem with the FDM presented here is on hand. Future works will focus on homogeneous Navier-Stokes equations for the simulation of two-phase flows in Nuclear Power Plants.

NOMENCLATURE

˜✄

: subscript invoking data of the physical problem ✠ : physical domain

✠Figure 1 :

 1 Figure 1: Example of a physical domain ˜ immersed in a fi ctitious domain

  . The fictitious domain approaches presented here deal with Dirichlet, Robin or Neumann B.C. on an immersed interface without requiring neither the modification of the numerical scheme near the The 11 th International Topical Meeting on Nuclear Thermal-Hydraulics (NURETH-11) Popes' Palace Conference Center, Avignon, France, October 2-6, 2005 3/17 immersed interface nor the use of Lagrange multipliers. Since the fictitious problem ✂ ✁

Figure 2 :

 2 Figure 2: Immersion of the unit quarter disk in the unit square

✠h

  contains additional terms that impose the initial B.C.. Given a rectangular mesh ✆ h of ✠ , the problem solved on the discrete fictitious domain ✠ h has the following generic form :

  ) are arbitrary extensions in ✠ h of the initial coefficients lying The 11 th International Topical Meeting on Nuclear Thermal-Hydraulics (NURETH-11) Popes' Palace Conference Center, Avignon, France, October 2

  the B.C. on ☛ e , for a penalization of the spread interface, the Dirichlet B.C. u h ✞ u D must only be imposed on the nodes lying on ✄ h ☎ ✁ ☛ e . The B.C. on the rest of ☛ e have no influence on the solution obtained in the physical domain. To penalize the exterior domain, the B.C. on the whole exterior boundary ☛ e must be Dirichlet B.C. u h ✥ ✂ e ✞ u D . ✂ Robin case: We consider the continuous transmission problem between ✠ and ✠ e . The addition of the two weak formulations of the two diffusion subproblems defined on respectively ✠ and ✠ e leads to the following problem :

  Only the B.C. of the nodes lying on ✄ h ☎ ✕ ☛ e has an effect on the solution obtained in the physical domain. These B.C. must be homogeneous Neumann B.C. in order to have an external flux equal to zero. However the B.C. on the whole boundary ☛ e can't be homogeneous Neumann B.C. since the problem solved on ✠ becomes in this case bad-posed. The solution obtained on the fictitious domain won't be unique. So the B.C. on ☛ e are separated into two B.C.:

  ✆ h . Let n✂ be the unit normal vector on ✁ oriented from K to the exterior. For a function ✄ in L 2 ✂ ✠ ✝ , let ✄ ✞ and ✄ ✠ be the traces of ✄ on each side of ✁ . We define the arithmetic mean of traces of ✄ as ✄ ✥ ✂ ✞ ✂ the jump of traces of ✄ on ✁ , oriented by n✂ , as ★ ★ ✄ ✦ ✦ ✂ ✞ ✂ ✄ ✞ ✁ ✞✄ ✠ ✝ . Then, the following fracture conditions lie on each side ✁ ✂ S :

  The 11 th International Topical Meeting on Nuclear Thermal-Hydraulics (NURETH-11) Popes' Palace Conference Center, Avignon, France, October 2-6, 2005 7/17 (a) Exterior thin interface (b) Interior thin interface (c) "Cut" thin interface Figure 3: Examples of thin approximated interface 2.2.2 Treatment of the initial B.C. We are interesting in the B.C. lying on the immersed approximated interface ☞ h . On the sides ✁ ✆ ☎ ☞ h , ✁ ☎ S ✁ ✂ ☛✄✂ ☛ e ✝ , we set:

  The B.C. on ☛ e ✁ ☞ h are free (without effect on the solution on the physical domain). On ☞ h ☛ e , we set the Robin B.C. of the initial problem.

Figure 4 :

 4 Figure 4: Estimation of the discretization error for a Dirichlet B.C.

Figure 5 :

 5 Figure 5: Absolute error for the exterior penalization -F.E Scheme -Dirichlet case -32x32 mesh

Figure 6 :

 6 Figure 6: Relative error norms with or without refi nement -F.E. Scheme -Exterior penalization -Dirichlet case

(Figure 7 :(Figure 8 :

 78 Figure 7: Estimation of the discretization error for a Robin B.C.

(Figure 9 :

 9 Figure 9: Relative error norms with and without refi nement -F.E. Scheme -Robin case

Figure 10 :

 10 Figure 10: Relative error norms for a weighted source term f

Figure 12 :

 12 Figure 12: Relative error norms for a surface integral approximation

✂

  numerical results for Dirichlet and Robin B.C. on the immersed interface. These two cases provide satisfactory results with errors in ✂ h✝ .

  used to impose the immersed flux

  the element K

	The 11 th International Topical Meeting on Nuclear Thermal-Hydraulics (NURETH-11)	11/17
	Popes' Palace Conference Center, Avignon, France, October 2-6, 2005

  The 11 th International Topical Meeting on Nuclear Thermal-Hydraulics Popes' Palace Conference Center, Avignon, France, October 2-6, 2005 source term. Locally, on each element of the spread interface, the balance of flux of the initial problem is approximatively respected. However, we are conscious that this correction isn't equivalent to a physical domain integration of the diffusion term. Results are presented on Figure11.

	14/17
		Relative error norms (discrete L2 norm)
	10 0 ||u-uh||/||u||		Robin B.C. -F.E. Scheme -a and f weighted or not
	10 -1
	10 -2
			Variable epsilon (lin.approx.) => order =0.86
			Lin. approx. -f weighted => order =1.
			Lin. Approx. -a and f weighted => order = 1.84
	10 -3	10 -2	10 -1	discretization step

Figure 11: Relative error norms for a weighted diffusion coeffi cient a