C Chateau

L Gélébart

M Bornert

J Crépin

E Boller

C Sauder

W Ludwig

In-situ X-ray microtomography characterization of damage in SiC f /SiC minicomposites

Keywords: A. Ceramic-matrix composites (CMCs), B. Matrix cracking, C. Damage mechanics, D. X-ray computed microtomography (µ-CT)

The purpose of the present study is to characterize matrix crack propagation and fiber breaking occurrences within SiC/SiC minicomposite in order to validate later on a multiscale damage model at the local scale. An in-situ X-ray microtomography tensile test was performed at the European Synchrotron Radiation Facility (ESRF, ID19 beamline) in order to obtain 3-dimensional (3D) images at six successive loading levels.

Results reveal a slow and discontinuous propagation of matrix cracks, even after the occurrence of matrix crack saturation. A few fiber failures were also observed.

However, radiographs of the whole length (14 mm) of the minicomposites under a load and after the failure were more appropriate to get statistical data about fiber breaking. Thus, observations before the ultimate failure revealed only a few fibers breaking homogenously along the minicomposite. In addition, an increase in fiber breaking density in the vicinity of the fatal matrix crack was observed after failure. These experimental results are discussed in regards to assumptions used in usual 1dimensional (1D) models for minicomposites.

Introduction

SiCf/SiC composites are prospective candidates for functional uses in future nuclear reactors -such as gas cooled fast reactors (GFR) -because of their favorable mechanical properties at high temperatures and after irradiation. The composites under investigation are made from a 2D fibrous preform composed of the new near-stoechiometric SiC fibers (Hi-Nicalon type S or Tyranno SA3 fibers), using the Chemical Vapor Infiltration (CVI) process. The material exhibits a nonlinear behavior due to the accumulated damage occurring between and inside the woven tows, such as through matrix cracking, fiber/matrix debonding as well as fiber breaking. Thus, a characterization of damage mechanisms within the tow is required to build and validate at local scale a multiscale predictive model. Due to their simple geometry, minicomposites (unidirectional composites containing a single bundle of fibers) are well suited to study these mechanisms. They are also frequently used to optimize the fiber/matrix interphase which dictates the matrix crack deflection along the fibers and consequently the nonlinear behavior of the composite [START_REF] Morscher | Tensile stress rupture of SiCf/SiCm minicomposites with carbon and boron nitride interphases at elevated temperatures in air[END_REF][START_REF] Bertrand | SiC/SiC minicomposites with nanoscale multilayered fibre coatings[END_REF][START_REF] Sauder | Influence of interface characteristics on the mechanical properties of Hi-Nicalon type-S or Tyranno-SA3 fiber-reinforced SiC/SiC minicomposites[END_REF].

Several 1D statistical models of the evolutional damage have previously been reported [START_REF] Baxevanakis | Fracture statistics of a unidirectional composite[END_REF][START_REF] Hui | An exact closed-form solution for fragmentation of weibull fibers in a single filament composite with applications to fiber-reinforced ceramics[END_REF][START_REF] Lissart | Damage and failure in ceramic matrix minicomposites: Experimental study and model[END_REF][START_REF] Dibenedetto | Statistical simulation of fiber fragmentation in a single-fiber composite[END_REF][START_REF] Curtin | Modeling brittle and tough stress-strain behavior in unidirectional ceramic matrix composites[END_REF][START_REF] Lamon | Stochastic approach to multiple cracking in composite systems based on the extreme-values theory[END_REF][START_REF] Castelier | Three consistent approaches of the multiple cracking process in 1D composites[END_REF][START_REF] Phoenix | Statistical issues in the fracture of brittle-matrix fibrous composites[END_REF]. They are based on matrix and fiber failure probability laws (such as the Weibull law) and are complemented by a stress redistribution assumption in the vicinity of matrix cracks. These models may lead to satisfactory predictions of the macroscopic response. However, microscopic phenomena cannot be fully validated because of the lack of experimental damage characterization. In fact, if the qualitative damage evolution is accepted, then observations reported in other literature were limited to the sample surface and were mostly collected after the ultimate failure [START_REF] Morscher | Tensile stress rupture of SiCf/SiCm minicomposites with carbon and boron nitride interphases at elevated temperatures in air[END_REF][START_REF] Lissart | Damage and failure in ceramic matrix minicomposites: Experimental study and model[END_REF][START_REF] Bertrand | Hi-Nicalon/SiC minicomposites with (pyrocarbon/SiC)(n) nanoscale multilayered interphases[END_REF][START_REF] Martinez-Fernandez | Room and elevated temperature tensile properties of single tow Hi-Nicalon, carbon interphase, CVI SiC matrix minicomposites[END_REF].

The purpose of this article is to present an experimental characterization of damage in SiCf/SiC minicomposites under a tensile load using X-ray microtomography. As reviewed by Stock [START_REF] Stock | Recent advances in X-ray microtomography applied to materials[END_REF], microtomography has been successfully used in material research. In particular, it has been used to study damage or fatigue cracks [START_REF] Buffière | Three-dimensional visualisation of fatigue cracks in metals using high resolution synchrotron X-ray micro-tomography[END_REF] in several materials such as fiber reinforced metal [START_REF] Mcdonald | X-ray tomographic imaging of Ti/SiC composites[END_REF][START_REF] Breunig | MicroCT (microtomography) quantification of microstructure related to macroscopic behaviour -Part 2 -Damage in SiC-Al monofilament composites tested in monotonic tension and fatigue[END_REF][START_REF] Hung | Fatigue crack growth and load redistribution in Ti/SiC composites observed in situ[END_REF] or polymer [START_REF] Schilling | X-ray computed microtomography of internal damage in fiber reinforced polymer matrix composites[END_REF][START_REF] Aroush | A study of fracture of unidirectional composites using in situ high-resolution synchrotron X-ray microtomography[END_REF][START_REF] Moffat | Micromechanisms of damage in 0 degrees splits in a [90/0](s) composite material using synchrotron radiation computed tomography[END_REF] matrix composites, aluminium alloys [START_REF] Guvenilir | New direct observations of crack closure processes in Al-Li 2090 T8E41[END_REF][START_REF] Ferrié | Fatigue crack propagation: In situ visualization using X-ray microtomography and 3D simulation using the extended finite element method[END_REF] or polymers (PMMA) [START_REF] Muller | Visualising complex morphology of fatigue cracks in voxel based 3D datasets[END_REF]. However, tomography applied to SiCf/SiC composites has been limited to porosity observations [START_REF] Lee | Pore geometry in woven fiber structures: 0 degrees/90 degrees plain-weave cloth layup preform[END_REF][START_REF] Morales-Rodriguez | Porosity analysis of long-fiber-reinforced ceramic matrix composites using X-ray tomography[END_REF][START_REF] Gélébart | X-ray tomographic characterization of the macroscopic porosity of chemical vapor infiltration SiC/SiC composites: Effects on the elastic behavior[END_REF], crack observations requiring a very high resolution because openings are smaller than 1 µm. In order to investigate matrix crack propagation through SiCf/SiC minicomposite sections, 3D images of a minicomposite under several tensile loading levels were acquired using the X-ray synchrotron source provided by the European Synchrotron Radiation Facility (ESRF). These images, reconstructed from a large number of radiographs, have been analyzed to detect matrix cracks within a small volume. Therefore, the morphology and kinetics of crack propagation through the minicomposite section can be described. In order to get statistical data on the fiber failure locations, fiber breaking has also been directly observed using a single radiograph of the entire sample. This was done for a single tensile loading level (about 80 % of the stress to failure) and after failure.

Material and methods

Minicomposites

The studied minicomposites were made [START_REF] Sauder | Influence of interface characteristics on the mechanical properties of Hi-Nicalon type-S or Tyranno-SA3 fiber-reinforced SiC/SiC minicomposites[END_REF] from a fibrous yarn constructed from 500 Hi-Nicalon type S fibers, with an average diameter of 13 m. The 100 nm interphase (pyrocarbon) and the SiC matrix were deposited on the fibers using the CVI process.

The residual porosity due to the CVI process, and fiber fractions were estimated at 0.12 (±0.04) and 0.58 (±0.09) from polished cross-section pictures (figure 1a).

Like the woven composite, minicomposites have an elastic, damageable behavior. Their macroscopic behavior (figure 2) follows typical successive steps in accordance with damage evolution [START_REF] Sauder | Influence of interface characteristics on the mechanical properties of Hi-Nicalon type-S or Tyranno-SA3 fiber-reinforced SiC/SiC minicomposites[END_REF][START_REF] Bertrand | Hi-Nicalon/SiC minicomposites with (pyrocarbon/SiC)(n) nanoscale multilayered interphases[END_REF]. The first one is an elastic domain of the minicomposite behavior: no cracking occurs. A second nonlinear step is associated with the matrix cracking (figure 1b) until saturation of crack number density (reached for a total strain of about 0.3%). A second linear domain associated with the additional elastic deformation of fibers is then observed after matrix crack saturation. The final step is characterized by a slight nonlinearity associated with fiber breaking just before the ultimate failure (close to 0.7%).

Experimental procedure

The in-situ microtomography tensile test was carried out on the ID19 beamline at the ESRF, in Grenoble, France. A specific in-situ tensile testing machine dedicated to ID19 was used to manually load the specimen (called specimen #1). The minicomposite was glued onto aluminium tabs and had a gauge length of 14 mm. Only the load was monitored using a 500 N load cell. The test was interrupted at six successive loading levels (50, 68, 74, 86, and 92 N) to record the tomography images. The corresponding load levels are reported on the typical load-strain curve presented in figure 2, and were obtained from another sample of the same batch with a classical macroscopic device.

Microtomography consists in recording a set of radiographs acquired at high resolution and at various angular positions of the sample with respect to the X-ray beam [START_REF] Salvo | X-ray micro-tomography an attractive characterisation technique in materials science[END_REF].

Appropriate algorithms are then used to reconstruct the 3D image from this set of radiographs. A high resolution was necessary to observe both cracks and the microstructure of the tow. Such a resolution with moderate acquisition times could only be reached through a synchrotron radiation, which gives a monochromatic, parallel and high intensity beam. The highest resolution provided on the ID19 beamline at the ESRF is a voxel (volumetric pixel) size of 0.28 m, identical in all three directions. To observe a 1.65 mm total length, three acquisitions at three successive axial positions were required. In addition, these acquisitions, or scans, were performed at four distinct distances (8, 14, 26 and 36 mm) from the sample to the camera (a Fast Readout Low Noise -FreLoN -14 bit CCD camera with a resolution of 2048 × 2048 pixels [START_REF] Labiche | The fast readout low noise camera as a versatile Xray detector for time resolved dispersive extended X-ray absorption fine structure and diffraction studies of dynamic problems in materials science, chemistry, and catalysis[END_REF]), using a 20.5 keV energy beam. From these four scans, two 3D images were reconstructed from the radiographs recorded at the two shorter distances, using the standard absorption mode (filtered back projection reconstruction) and a third one using the holotomographic mode (based on phase contrast). 3D holotomographic images were reconstructed by combining the four scans [START_REF] Cloetens | Holotomography: Quantitative phase tomography with micrometer resolution using hard synchrotron radiation X rays[END_REF], using the specific algorithm proposed by Langer [START_REF] Langer | Quantitative comparison of direct phase retrieval algorithms in in-line phase tomography[END_REF][START_REF] Langer | Regularization of phase retrieval with phaseattenuation duality prior for 3-D holotomography[END_REF]. Finally, these scans were performed for each loading level (50, 68, 74, 80, 86 and 92 N). In total, the entire experiment required 72 scans (3 fields of view x 4 distances x 6 loading levels) with each scan consisting of 1500 radiographs.

Due to the length of acquisition time for each scan (25 min), it was not possible to observe the entire gauge length of the specimen. However, a simple radiograph for the whole length was acquired at 92 N. Finally, another sample (called specimen #2) was loaded until failure (115 N) with the same tensile device. A radiograph of the whole specimen was acquired after failure. To summarize, three distinct kinds of observation have been made: A. the 3D images of a 1.65 mm long region of the specimen #1 observed at six successive loading levels, B. the radiograph of the whole specimen #1 (about 14 mm long) under tensile loading at 92 N, C. the radiograph of the longest part of the broken specimen #2 (about 12 mm).

Damage detection

Matrix cracking (obs.A)

One goal of this study is to characterize the matrix crack morphology and its evolution with a load increase. Indeed, as presented in figure 3 (the detection procedure is detailed below), matrix cracks have a specific morphology which requires accurate locating.

Such a result was obtained by spotting grey level variations in 3D images caused by the crack (figure 4). These variations could be visually noticed by routine observation of transverse slices of the tow. In the following, a "crack width" (along the axial direction) is defined from the number of transverse slices crossed by the crack (see figure 3). Note that crack openings cannot be characterized with accuracy because image resolution is not sufficient.

Moreover, the crack effects are more noticeable on the images reconstructed from the second distance radiographs than from the first. In fact, it is now well known that when a synchrotron source is used, an additional phase contrast due to diffraction effects appears on projections that improves the efficiency of damage detection [START_REF] Buffière | Three-dimensional visualisation of fatigue cracks in metals using high resolution synchrotron X-ray micro-tomography[END_REF][START_REF] Cloetens | Observation of microstructure and damage in materials by phase sensitive radiography and tomography[END_REF]. This contrast, invisible at the zero distance, is enhanced as the sample/detector distance increases. However, an increase in distance also emphasizes diffraction fringes on the crack edges, making a quantitative observation of the crack more difficult. Although the crack is less noticeable on holotomographic reconstructions, diffraction fringes do not appear on these reconstructed images (figure 4) as the reconstruction procedure takes into account their physical origin [START_REF] Buffière | Three-dimensional visualisation of fatigue cracks in metals using high resolution synchrotron X-ray micro-tomography[END_REF][START_REF] Cloetens | Holotomography: Quantitative phase tomography with micrometer resolution using hard synchrotron radiation X rays[END_REF]. In order to precisely locate the crack within the section, a specific procedure was established to automatically detect the cracks in the volume by combining the three distinct reconstructions.

Considering a sub-volume containing a matrix crack (figure 1), the detection procedure is based on the processing of the grey level variation in the axial (Z) direction for all (X,Y) positions on the transverse plane. An example of such an evolution for a given (X,Y) position is provided in figure 5a. The crack leads to an important variation of the grey levels on a relatively short height among the three referenced reconstructions. The aim of the automatic procedure is to spot and locate this fluctuation in the axial direction, while reducing the detection of artifacts such as microstructure changes (matrix/pore interface) or image noise effects. Thus, a specific filter has been developed (figure 5b). It corresponds to the ratio of the norm of the grey level gradient along the Zdirection to the average derivative in a vicinity along Z (defined in figure 5a). In addition, the derivative was smoothed prior to the filtering in order to reduce the noise effect, using a centered moving average (10 voxels window along the Z-direction). A crack was then detected when the three global maxima (corresponding to the three reconstructions) of the filtered grey levels were separated by less than 10 voxels, ie 2.8 µm, which specifies the accuracy of the crack location. Otherwise, the material is considered healthy for the (X,Y) position and for all Z slices of the sub-volume. All (X,Y) positions are processed independently this way.

A combination of the three reconstructed volumes allows for reducing the detection of artifacts, like diffraction fringes at other interfaces. Moreover, the few number of detections located in the porous phase were removed.

Besides the main cracks that are detectable using the automatic procedure, there are small matrix cracking zones located on the periphery of the sample (isolated fiber group) or within the minicomposite. These are called minor cracks. They lead to a far smaller grey level variation along the Z axis, similar in intensity to image noise. As a result, the automatic detection procedure could not be used. However, because of their small size (about 5% of the global section) and their limited number, a manual axial location based on a direct observation of the images was sufficient to characterize them.

The tomographic images resulting from observation A also show a few fiber failures which were easy to detect visually. Nevertheless, the observed field of view was far too small and the fiber breaking density too low to get statistically representative data about fiber breaking. Fiber/matrix debonding was not observable on these images.

Fiber breaking (obs.B and C)

In order to study fiber break distribution under load (92 N) and after failure, simple radiographs were useful to locate both fiber failures and main matrix cracks along the entire minicomposite length, as shown on figure 6a. While a few minor cracks were detected (they were harder to find on those images), they were not reported here. The detection was limited to visual observation leading to an inventory of fiber break and main matrix crack positions. Fiber crack openings and crack widths, as defined in section 1.3.1, can also be estimated (figure 6).

Results and discussion

Matrix cracking

Matrix crack distribution

The axial distribution of main matrix cracks and fiber breaks (as defined in section 1.3) are presented in figure 7 from radiograph observation performed at 92 N of the entire minicomposite (obs. B). At this loading, the matrix cracking can be assumed to be completed (see figure 2) so the average intercrack distance at saturation of crack number density was estimated to be equal to 250 µm. This value is consistent with crack spacing distances measured on similar minicomposites after classical tensile tests [START_REF] Sauder | Influence of interface characteristics on the mechanical properties of Hi-Nicalon type-S or Tyranno-SA3 fiber-reinforced SiC/SiC minicomposites[END_REF] (173 µm to 350 µm) .The field of view of the tomography observations (obs. A) is also reported in figure 7. Before studying matrix crack morphology and propagation in detail, the evolution of the damage axial distribution was investigated by visual detection in the CT-reconstructions for the six different loading levels (figure 8). A first observation is that most cracks evolve when the load is increasing, in the sense that their width, as defined before, increases. Moreover, some of the minor cracks developed into main cracks, but the majority remained isolated.

The damage distribution is then compared with the effective section and surface porosity fraction variations in the axial direction (figure 8). The purpose is to study a potential link between damage location and these global characterizations of the microstructure. Both were estimated from a threshold of the holotomographic image.

The porosity of a section is here defined as the ratio of the closed pore area over the largest connex part of the composite section. Open pores on the periphery are not taken into account. These measurements are sensitive to the threshold choice. Its impact (relative error) was estimated at ±2% for the effective section and ±10% for the porosity fraction. The emphasis here is on the porosity variation along the minicomposite which does not significantly depend on the threshold. Note also that the sharp porosity variations correspond to openings or closings of peripheral pores. These observations did not evidence any clear link between crack locations and these global microstructure quantifications. Crack locations do generally not coincide with local extrema of effective section or porosity. Crack 2 is an exception, as it coincides with a local minimum of the effective section and thus a local maximum of average axial stress, but this is far from a general features.

Matrix crack propagation

The automatic procedure described in section 1.3 was used to detect the six main cracks, numbered in figure 8, for all of the six loading levels applied during observation A. Six sub-volumes, centered on the six main cracks observed at 92 N, were defined to compute the detection. Figure 9 describes the 3D geometry of the six detected cracks, in the form of the height Z of each crack as a function of the (X,Y) positions on the transverse plan. Only loading levels indicating crack propagation are presented (for example crack 2 did not propagate between 68 N and 74 N).

Two preliminary comments can be made before going into detail. Firstly, projections of cracks 1 and 6 are less clear than others. Contrast variations (like the one presented in figure 5) were far less pronounced for these two cracks. This is probably a consequence of a smaller crack opening. Secondly, the minor crack 3b presented in figure 10 and main crack 3a are complementary. As presented in figure 8, the two cracked zones were independent before the last loading step.

Concerning the crack morphology at 92 N, all six main cracks spread across the entire minicomposite section. They are not flat and typically follow a spiral shape around the fiber direction, except for crack 4 which presents an axial symmetry. Note also that the crack elongation in the fiber direction -from 40 µm up to more than 100 µm -is quite large with respect to the intercrack distance (250 µm at saturation).

Pertaining to crack propagation, two or three distinct propagation states were observed depending on the crack. For example, cracks 2, 3, 4 and 6 propagation states followed three typical steps: crack initiation in a localized peripheral zone of the minicomposite, crack propagation on the entire periphery, and finally, propagation towards the section center. Crack 5 also initialized on the periphery, but the intermediate propagation state was not caught during observations. Lastly, although it was not revealed by the detection process, crack 1 initiation was visually detected and occurred inside the minicomposite section at 86 N.

Moreover, cracks 2, 3 and 4 did not propagate between 68 and 74 N. This revealed a discontinuous propagation. Furthermore, even if the matrix cracking saturation was supposed to be reached (ie no additional crack appears), the propagation was not achieved at 74 N. These results contradict the usual 1D modeling of minicomposite [START_REF] Lissart | Damage and failure in ceramic matrix minicomposites: Experimental study and model[END_REF][START_REF] Lamon | Stochastic approach to multiple cracking in composite systems based on the extreme-values theory[END_REF], which assumes that the matrix cracks suddenly spread across the entire section with a uniform crack opening (ie propagation stage is ignored).

Fiber breaking

Radiographs resulting from observations B and C lead to statistical data about fiber breaking along the minicomposite under loading (92 N) and after failure. Following figure 7, fiber breaks appeared homogeneously along the minicomposite, mostly located in the vicinity of the matrix cracks. The fiber break density was rather low, with 4.8 broken fibers per mm (ie less than 1% fibers per mm).

Considering post-failure observations without counting final breaks (figure 11), approximately 220 fiber breaks were observed in the 12 mm long half minicomposite.

The break density is far less homogeneous (figure 12). In fact, most fiber breaks are located within 2 mm of the failure zone. On the other hand, the fiber break distribution is quite homogeneous at other locations. It has a density of approximately 5 breaks per mm (calculated on region III in figure 12), very similar to the density at 92 N.

Additionally, observation A shows that the first fiber breaks observed with 3D tomography appeared at 80 N (figure 8). This load matches the final stage of matrix crack propagation inside the minicomposite.

To summarize, it appears that fiber breaking occurs at the completion of matrix cracking propagation with a rather low density and a uniform distribution with preferred locations near matrix cracks. At the load resulting in ultimate failure, fiber breaking was centralized within an area (about 4 mm long) surrounding the ultimate failure location.

Finally, fiber crack openings reported on figure 12 are 7 µm long in average far from the failure zone. Fiber breaks located near the failure zone exhibit larger openings (up to 50 µm).

In agreement with these observations, two different scenarios could be proposed. On the one hand, a random uniform fiber breaking develops slowly as the stress increases and the localization of the fiber density observed on figure 12 is the consequence of dynamic effects during failure. On the other hand, the random fiber breaking starts with a uniform distribution, then a localization of the break density occurs and leads to the ultimate failure. This second scenario is in better agreement with the final slight nonlinearity observed on macroscopic curves, but additional experiments are required to conclude. However, in both cases, this characterization shows that one fiber can fail twice along the minicomposite. This has also been observed on a fiber in obs.A.

These experimental results are inconsistent with the assumption of the classical fiber bundle model [START_REF] Daniels | The statistical theory of the strength of bundles of threads .1[END_REF] used in some 1D models [START_REF] Lissart | Damage and failure in ceramic matrix minicomposites: Experimental study and model[END_REF][START_REF] Lamon | A micromechanics-based approach to the mechanical behavior of brittlematrix composites[END_REF] which supposes that once broken, a fiber does not participate anymore to the load transfer throughout the whole composite (and as a consequence breaks only once). Fiber break observations are rather in line with works which take into account a sliding length around a break necessary for the fiber to recover back its previously carried stress [START_REF] Curtin | Theory of mechanical-properties of ceramic-matrix composites[END_REF][START_REF] Phoenix | Scalings in fracture probabilities for a brittle matrix fiber composite[END_REF][START_REF] Phoenix | Statistical issues in the fracture of brittle-matrix fibrous composites[END_REF]. Moreover, those works also emphasize a localization region around some plane weaker than the others along the composite where fiber breakage will continue once final breakdown occurs there [START_REF] Curtin | Fiber pull-out and strain localization in ceramic matrix composites[END_REF][START_REF] Curtin | The tough to brittle transition in brittle matrix composites[END_REF][START_REF] Phoenix | Size effects in the distribution for strength of brittle matrix fibrous composites[END_REF].

Concluding remarks

Standard 2D techniques are not sufficient to characterize damage mechanisms such as matrix crack propagation and fiber breaking occurring within the material. Therefore, an in-situ tensile test was performed on a SiCf/SiC minicomposite under X-ray microtomography to observe damage evolution. The matrix crack detection required the use of a high resolution equipment available through synchrotron radiation at the ESRF. Two observation techniques were used and resulted in the following damage characterization:

-The tomography (absorption and holotomographic mode) performed on a minicomposite at 6 successive loading levels was necessary to observe matrix crack evolutions within the sample. A few fiber breaks may also be observed.

-Simple radiograph acquisition was much faster than tomography acquisition.

Therefore, the whole minicomposite observations could be performed under loading and after failure. Matrix crack and fiber break characterizations are limited to their spatial location along the sample, but it continues to be a useful way to get statistical data about fiber breaking. Moreover, the resolution was sufficient to estimate fiber break openings.

These observations provided a deeper insight on damage mechanisms:

-Besides well known transverse cracks and fiber breaks, 3D images revealed minor matrix cracks. These small matrix cracking zones located within the minicomposite could become a transverse crack or could remain localized. However, the cracked sections are so small that the influence on the macroscopic damageable behavior is likely to be insignificant compared to main cracks.

-A specific automatic procedure was developed to detect and locate main cracks along the minicomposite. Most of them are not flat but typically follow a spiral shape which is quite elongated in the fiber direction. In addition, the matrix crack propagates slowly and discontinuously within the bundle of fibers, even after the cracking saturation. They mostly appear at first on the peripheral zone, and then propagate towards the center of the section. This experimental evidence is contrary to common assumption of 1D models which ignore crack propagation. Such assumptions that are appropriate for microcomposites (containing a single fiber) [START_REF] Castelier | Three consistent approaches of the multiple cracking process in 1D composites[END_REF], may be insufficient for minicomposites.

-Fiber breaking seems to begin immediately after matrix crack propagation ends. Fibers fail, at first homogeneously, and are typically located in the vicinity of matrix cracks. Fiber break density stabilizes around 5 fiber failures per mm, except in a short zone (a few millimeters) surrounding the ultimate failure where the break number density is much greater. Fiber openings are also much larger in this area. Two scenarios have been proposed to explain these observations, involving dynamic effects induced by the ultimate failure or fiber breaking localization leading to the ultimate failure. Moreover, experimental results suggest selecting models based on frictional load-sharing rather than neglecting the contribution of broken fibers to the load transfer.

Figure captions

Figure 1 :

 1 Figure 1: SEM micrographs of (a) minicomposite polished cross section; (b) cracked

Figure 2 :

 2 Figure 2: Macroscopic tensile behavior obtained with a macroscopic device. The

Figure 3 :

 3 Figure 3: (a) Sub-volume containing a matrix crack; (b-c) Detected matrix crack from

Figure 4 :

 4 Figure 4: Reconstruction of a transverse slice in absorption contrast (first two distances)

Figure 5 :

 5 Figure 5: Filtering of grey levels :(a) grey level evolution in the fiber direction through a

Figure 6 :

 6 Figure 6: (a) Radiograph of the minicomposite at a tensile load of 92 N (b) Zoom on a

Figure 7 :

 7 Figure 7: Damage location along the entire minicomposite at 92 N as detected from the

Figure 8 :

 8 Figure 8: Damage location along the minicomposite in the area observed by

Figure 9 :

 9 Figure 9: Results of main crack detections (projections onto the transverse plane). Color

Figure 10 :

 10 Figure 10: Projection onto the transverse plane of crack 3b, which is an extension of the

Figure 11 :

 11 Figure 11: Radiographs after failure located near (I.) and far from the ultimate failure

Figure 12 :

 12 Figure 12: Fiber break local density and opening as a function of the distance to the

  Figure 1

Figure 2 Figure 3

 23 Figure 2

Figure 5 Figure 6 0

 56 Figure 5

Figure

  Figure 7

Figure 8 Figure 9

 89 Figure 8

Acknowledgement

The authors are very grateful to M.Langer for his help in reconstructing holotomographic images. The ID19 team is also acknowledged for its help during the experiment at the ESRF.


Stéphane Brull
email: stephane.brull@math.u-bordeaux1.fr

Fabrice Deluzet
email: fabrice.deluzet@math.univ-toulouse.fr

Alexandre Mouton
email: alexandre.mouton@math.univ-lille1.fr

Numerical resolution of an anisotropic non-linear diffusion problem *

Keywords: Anisotropic diffusion problems, Singular perturbation, Asymptotic-Preserving schemes AMS Subject Classification 35J60, 35J62, 65M06, 65M12, 65N06, 65N12

This paper is devoted to the numerical resolution of an anisotropic non-linear diffusion problem involving a small parameter ǫ, defined as the anisotropy strength reciprocal. In this work, the anisotropy is carried by a variable vector function b. The equation being supplemented with Neumann boundary conditions, the limit ǫ → 0 is demonstrated to be a singular perturbation of the original diffusion equation. To address efficiently this problem, an Asymptotic-Preserving scheme is derived. This numerical method does not require the use of coordinates adapted to the anisotropy direction and exhibits an accuracy as well as a computational cost independent of the anisotropy strength.

Introduction

The aim of this paper is to build an efficient numerical method for solving an anisotropic diffusion problem where the anisotropy is carried by a vector b. This work is motivated by investigations of strongly magnetized plasmas, more specifically the study of the Euler-Lorentz model in a low Mach number regime and in the presence of a large magnetic field. This framework is characteristic of the magnetically confined plasma fusion [START_REF] Degond | Asympotic-preserving schemes for fluid models of plasmas, in Collection "Panorama et Synthèses[END_REF][START_REF] Hazeltine | Plasma Confinement[END_REF][START_REF] Miyamoto | Controlled fusion and plasma physics[END_REF]. In this context, the asymptotic parameter ǫ represents the gyro-period of particles as well as the square root of the Mach number, the vector field b being the magnetic field direction. Therefore the ǫ values can be very small in some sub-regions of the computational domain where the magnetic field is large, inducing then a severe anisotropy of the medium, while being large in other sub-domains for intermediate and small strength of the magnetic field. Another important property of this system is the time dependence of the magnetic field defining the anisotropy direction. These two main characteristics define the framework of the present paper whose purpose is to design a numerical scheme for anisotropy ratios ranging from ǫ ≪ 1 to ǫ ∼ O(1) and for a time varying anisotropy direction. In order to address efficiently these requirements, the numerical method should not rely on a coordinate system adapted to this anisotropy direction. The use of adapted coordinates would imply mesh modifications accordingly to the evolution of b, an intricate and expensive procedure we wish to avoid. Thus, the numerical method introduced here will carry out the anisotropic non-linear diffusion problem on a mesh independent of the anisotropy direction. This scheme will be detailed on the following model problem

           -∇ x • H ǫ (b ⊗ b) ∇ x p ǫ -S ǫ ǫ + g ǫ (p ǫ ) = f ǫ , in Ω, H ǫ (b ⊗ b) ∇ x p ǫ -S ǫ ǫ
• ν ≡ 0, on ∂Ω.

(1.1)

In this system, Ω is a bounded subset of R d (d = 1,2,3), ǫ ≥ 0 is a fixed constant parameter and, for any x ∈ ∂Ω, ν = ν(x) stands for the unit outward normal vector. ∇ x and ∇ x • stand for the gradient and the divergence operators with respect to the space variable x. We assume that H ǫ : Ω → R * + , b : Ω → R d -{0}, f ǫ : Ω → R, S ǫ : Ω → R d are given and the unknown of the problem is the function p ǫ : Ω → R. The tensor product of two vectors u and v is denoted u ⊗ v. Finally, we assume that, for any ǫ, the function p → g ǫ (p) is strictly increasing and can be non-linear. This equation is well suited for the plasma fusion context above depicted. It allows the computation of the plasma pressure in order to guarantee that the forces vanish in the low Mach regime for strongly magnetized plasma. The function denoted g ǫ defines the internal energy of the fluid with respect to the pressure. This relation may be non-linear, which motivates the investigation of non-linear anisotropic problems. However, the derivation of this equation is out of the scope of the present paper and we refer to related works (see [START_REF] Brull | Asymptotic-Preserving scheme for a two-fluid Euler-Lorentz model[END_REF][START_REF] Brull | Numerical degenerate elliptic problems and their applications to magnetized plasma simulations[END_REF][START_REF] Degond | Asympotic-preserving schemes for fluid models of plasmas, in Collection "Panorama et Synthèses[END_REF]) for detailed explanations. Furthermore, we wish to present the numerical method in a context wider than the strict plasma context, since anisotropic diffusion problem are encountered in many applications. Good examples of these applications are, for instance, image noise filtering, convection dominated diffusion equations and more generally diffusion problem with strong medium anisotropies. The model equation (1.1) is representative of a large enough variety of problems, up to slight changes, and will be considered to detail the numerical method.

Developing an efficient numerical method to compute the solution of this diffusion problem, regardless to ǫ values, is a difficult task. Indeed, the limit ǫ → 0 is a singular limit for the problem (1.1), the diffusion equation degenerating into the following one

-∇ x • (H 0 (b ⊗ b)(∇ x p0 -S 0 )) = 0, in Ω, (H 0 (b ⊗ b)(∇ x p0 -S 0 )) • ν ≡ 0, on ∂Ω. (1.2)
The system (1.2) is ill-posed, its solution being non-unique. More precisely, if p0 is a solution of (1.2) and c : Ω → R is a function verifying b • ∇ x c = 0 on Ω, then p0 + c defines a new solution of (1.2). However, p 0 the limit of p ǫ solution of (1.1) is uniquely defined by the limit problem as demonstrated in Section 2, but a direct discretization of the diffusion problem (1.1) gives rise to a linear system with a conditioning number that blows up for vanishing ǫ. This property has been outlined for numerical studies of elliptic equation singular perturbations (see [START_REF] Besse | Efficient numerical method for strongly anisotropic elliptic equations[END_REF][START_REF] Degond | An Asymptotic-Preserving scheme for an anisotropic elliptic problem[END_REF]).

To tackle this difficulty, an Asymptotic-Preserving (AP) scheme is introduced to compute the solution of the anisotropic diffusion problem for ǫ = O(1) and to capture p 0 , the solution of the limit problem, for small ǫ values. This property should be provided without any limitations on the discretization parameters related to the value of ǫ. These requirements are compliant with the properties of AP-schemes originally introduced in [START_REF] Jin | Efficient Asymptotic-Preserving (AP) schemes for some multiscale kinetic equations[END_REF] and developed in [START_REF] Lemou | A new Asymptotic-Preserving scheme based on micro-macro decomposition for linear kinetic equations in the diffusion limit[END_REF] for diffusive regimes of transport equations. These techniques have received numerous extensions to other singular perturbation problems: relaxation limits of kinetic plasma descriptions [START_REF] Crouseilles | An Asymptotic-Preserving scheme based on a micro-macro decomposition for collisional Vlasov equations: diffusion and high-field scaling limits[END_REF][START_REF] Filbet | A class of Asymptotic-Preserving schemes for kinetic equations and related problems with stiff sources[END_REF][START_REF] Filbet | An Asymptotic-Preserving Scheme for the ES-BGK model of the Boltzmann equation[END_REF], quasi-neutral limit of fluid and kinetic plasma models [START_REF] Belaouar | An asymptotically stable semi-lagrangian scheme in the quasi-neutral limit[END_REF][START_REF] Crispel | Quasi-neutral fluid models for current-carrying plasmas[END_REF][START_REF] Crispel | An Asymptotic-Preserving scheme for the two-fluid Euler-Poisson model in the quasi-neutral limit[END_REF][START_REF] Crispel | A plasma expansion model based on the full Euler-Poisson system[END_REF][START_REF] Crouseilles | Two-Scale Macro-Micro decomposition of the Vlasov equation with a strong magnetic field[END_REF][START_REF] Degond | Asymptotic-Preserving Particle-In-Cell method for the Vlasov-Poisson system near quasi-neutrality[END_REF][START_REF] Degond | Numerical approximation of the Euler-Poisson-Boltzmann model in the quasi-neutral limit[END_REF][START_REF] Degond | Analysis of an Asymptotic-Preserving scheme for the Euler-Poisson system in the quasi-neutral limit[END_REF], hydrodynamic low Mach number limit [START_REF] Degond | All speed scheme for the low Mach number limit of the isentropic Euler equation[END_REF][START_REF] Klar | An Asymptotic-Preserving numerical scheme for kinetic equations in the low Mach number limit[END_REF], radiative hydrodynamics [START_REF] Buet | Diffusion limit of the Lorentz model: Asymptotic-Preserving schemes[END_REF][START_REF] Buet | Asymptotic-Preserving and positive schemes for radiation hydrodynamics[END_REF], fluid and particle flows [START_REF] Carrillo | Simulation of fluid and particles flows: Asymptotic-Preserving schemes for bubbling and flowing regimes[END_REF] and strongly magnetized plasmas as well as heterogeneous media [START_REF] Besse | Efficient numerical method for strongly anisotropic elliptic equations[END_REF][START_REF] Brull | Numerical degenerate elliptic problems and their applications to magnetized plasma simulations[END_REF][START_REF] Brull | Asymptotic-Preserving scheme for a two-fluid Euler-Lorentz model[END_REF][START_REF] Degond | Dualitybased Asymptotic-Preserving method for highly anisotropic diffusion equations[END_REF][START_REF] Degond | An Asymptotic-Preserving method for highly anisotropic elliptic equations based on a micromacro decomposition[END_REF][START_REF] Degond | Asympotic-preserving schemes for fluid models of plasmas, in Collection "Panorama et Synthèses[END_REF][START_REF] Degond | An Asymptotic-Preserving scheme for an anisotropic elliptic problem[END_REF][START_REF] Degond | An Asymptotic-Preserving scheme for the Euler equations in a strong magnetic field[END_REF].

The Asymptotic-Preserving property of the presented method is obtained thanks to a decomposition of the solution, introduced in [START_REF] Degond | An Asymptotic-Preserving scheme for an anisotropic elliptic problem[END_REF] and also used in [START_REF] Besse | Efficient numerical method for strongly anisotropic elliptic equations[END_REF][START_REF] Brull | Numerical degenerate elliptic problems and their applications to magnetized plasma simulations[END_REF][START_REF] Degond | Dualitybased Asymptotic-Preserving method for highly anisotropic diffusion equations[END_REF]. It consists of the following identity p ǫ = π ǫ + q ǫ , π ǫ being the solution mean part, with respect to the anisotropy (b) direction, q ǫ the fluctuating part. These two components verify π ǫ ∈ K and q ǫ ∈ K ⊥ , K defining the functions constant along the b-direction, K ⊥ the functions of zero mean value along b. This decomposition was first developed for meshes adapted to the anisotropy direction [START_REF] Besse | Efficient numerical method for strongly anisotropic elliptic equations[END_REF][START_REF] Degond | An Asymptotic-Preserving scheme for an anisotropic elliptic problem[END_REF], for which, the discretization of K is straightforward. A direct discretization of the sub-space K ⊥ is, on the other side, much more intricate. This difficulty is overcome thanks to the introduction of a Lagrangian multiplier, in order to penalize the zero mean value property of the functions belonging to K ⊥ . The method is extended in [START_REF] Degond | Dualitybased Asymptotic-Preserving method for highly anisotropic diffusion equations[END_REF] for computations with meshes independent of the anisotropy direction. This is achieved by introducing two more Lagrangian multipliers to discretize the sub-spaces. The size of the linear system providing the problem solution is then significantly enlarged. However, this drawback may be corrected thanks to a slightly different decomposition. In [START_REF] Degond | An Asymptotic-Preserving method for highly anisotropic elliptic equations based on a micromacro decomposition[END_REF] the solution is decomposed in two non-orthogonal parts which allows the definition of two sub-spaces whose direct discretization is readily obtained without any Lagrangian multipliers. The size of the linear system obtained with this approach is considerably lowered compared to the previous method [START_REF] Degond | Dualitybased Asymptotic-Preserving method for highly anisotropic diffusion equations[END_REF]. This method has been extended in [START_REF] Mentrelli | Negulescu Asymptotic-Preserving scheme for highly anisotropic non-linear diffusion equations[END_REF][START_REF] Lozinski | Negulescu Highly anisotropic temperature balance equation and its asymptotic-preserving resolution[END_REF] to non-linear diffusion equations. The path followed in the present paper still relies on the decomposition in K and K ⊥ . However, the discretization of these sub-spaces is achieved using a differential characterization, similar to the one introduced in [START_REF] Brull | Numerical degenerate elliptic problems and their applications to magnetized plasma simulations[END_REF]. This finally allows the computation of the solution thanks to a second-order problem for π ǫ and a fourth-order problem for q ǫ . This latter problem, in the framework of Neumann boundary conditions considered in this paper, can be recast into two elliptic problems.

The method proposed finally reduces to the computation of three standard elliptic problems for which very efficient solvers can be used (for instance multigrid solvers). For the former approaches [START_REF] Degond | Dualitybased Asymptotic-Preserving method for highly anisotropic diffusion equations[END_REF][START_REF] Degond | An Asymptotic-Preserving method for highly anisotropic elliptic equations based on a micromacro decomposition[END_REF], the equations providing both components are not classical elliptic equations and the resolution of the linear system requires more sophisticated solvers. This complexity is resource demanding and may be challenging for realistic three-dimensional computations. Finally this paper also presents an extension to non-linear reaction diffusion problems, a class of problems that has never been investigated in the previous works.

The paper is organized as follows: in Section 2, the decomposition methodology is presented. The linear case, i.e. with g ǫ (p)(x) = G ǫ (x)p(x) where G ǫ : Ω → R * + is a given function sequence, is first investigated: more precisely, we describe the decomposition procedure in the specific case where G ǫ is a strictly positive constant denoted λ ǫ , then we generalize this procedure to any function G ǫ : Ω → R * + by using well-chosen Sobolev spaces. Finally the non-linear problems are addressed by invoking Gummel's iterative algorithm. Section 3 is devoted to presentation of the discretization. Finally, the efficiency of the numerical method is demonstrated in Section 4.

Scale separation and solution decomposition

In this section a scale separation is introduced to ensure the Asymptotic-Preserving property of the scheme. This is achieved by transforming the singular perturbation problem (1.1) into an equivalent system for which the limit ǫ → 0 is regular. For simplicity reasons, the linear case with constant G ǫ is first considered for detailing the decomposition method. In this framework, the singular nature of the limit ǫ → 0 is outlined and the limit problem, providing p 0 = lim ǫ→0 p ǫ , is stated. A development to linear cases with variable positive functions G ǫ is then presented and finally, thanks to Gummel's iterative method [START_REF] Gummel | A self-consistent iterative scheme for one-dimensional steady state transistor calculations[END_REF], the non-linear case is addressed by using a sequence of linear problems.

2.1 AP-scheme derivation for linear problems

2.1.1 A simplified framework: constant G ǫ
We assume here that the given sequence (g ǫ ) ǫ ≥ 0 is of the form

g ǫ (p)(x) = λ ǫ p(x),
where λ ǫ > 0 is a known constant for any ǫ ≥ 0. Then the diffusion problem (1.1)

writes -∇ x • (H ǫ (b ⊗ b)(∇ x p ǫ -S ǫ )) + ǫλ ǫ p ǫ = ǫf ǫ , in Ω, (H ǫ (b ⊗ b)(∇ x p ǫ -S ǫ )) • ν ≡ 0, on ∂Ω . (2.1)
The limit solution p 0 of the singular perturbation problem (2.1) verifies the limit problem

     -lim ǫ → 0 ∇ x • H ǫ (b ⊗ b) ∇ x p ǫ -S ǫ ǫ + λ 0 p 0 = f 0 , in Ω, (H 0 (b ⊗ b)(∇ x p 0 -S 0 )) • ν ≡ 0, on ∂Ω. (2.2)
This algebraic equation admits a unique solution under the assumption

H ǫ (b ⊗ b)(∇ x p ǫ -S ǫ ) = O(ǫ), (2.3)
a requirement that must be fulfilled by the numerical method. To ensure this property, the methodology consists in using a decomposition similar to that of [START_REF] Brull | Numerical degenerate elliptic problems and their applications to magnetized plasma simulations[END_REF][START_REF] Brull | Asymptotic-Preserving scheme for a two-fluid Euler-Lorentz model[END_REF][START_REF] Crouseilles | Two-Scale Macro-Micro decomposition of the Vlasov equation with a strong magnetic field[END_REF][START_REF] Degond | Dualitybased Asymptotic-Preserving method for highly anisotropic diffusion equations[END_REF][START_REF] Degond | An Asymptotic-Preserving scheme for an anisotropic elliptic problem[END_REF]. The solution is decomposed into π ǫ its mean part with respect to the anisotropy direction and q ǫ the fluctuating part, which exhibits the property to have a zero mean value along the anisotropy direction. These two functions verify π ǫ ∈ K and q ǫ ∈ K ⊥ , K being the kernel of the elliptic operator defined by equation 1.2. These properties are capitalized on, to isolate in the problem 2.1 the macro scale (providing π ǫ ) from the micro scale (giving q ǫ ) and thereby, build the Asymptotic-Preserving scheme. The main difficulty of the procedure lies in the characterization of the sub-spaces associated to the different scales. In [START_REF] Besse | Efficient numerical method for strongly anisotropic elliptic equations[END_REF][START_REF] Crouseilles | Two-Scale Macro-Micro decomposition of the Vlasov equation with a strong magnetic field[END_REF][START_REF] Degond | Dualitybased Asymptotic-Preserving method for highly anisotropic diffusion equations[END_REF][START_REF] Degond | An Asymptotic-Preserving scheme for an anisotropic elliptic problem[END_REF] the property of the functions populating K or K ⊥ are imposed by a penalization technique. The methodology developed in this paper operates a similar decomposition on to K and K ⊥ , but with a different characterization of these sub-spaces. Here, we shape the technique introduced in [START_REF] Brull | Numerical degenerate elliptic problems and their applications to magnetized plasma simulations[END_REF] for a very specific framework, in order to discriminate the functions in K and K ⊥ thanks to differential properties, providing thus, an easy discretization.

With this aim, we introduce the following Sobolev spaces:

V = p ∈ L 2 (Ω) : b • ∇ x p ∈ L 2 (Ω) , W = q ∈ L 2 (Ω) : ∇ x • (bq) ∈ L 2 (Ω) , W 0 = {q ∈ W : (bq) • ν ≡ 0 on ∂Ω} ,
and we define K ⊂ V as

K = {π ∈ V : b • ∇ x π = 0 on Ω} .
The goal is to reproduce the function decomposition into its mean and fluctuating parts. The functions of K correspond to the mean part and the complementary part is demonstrated to belong to K ⊥ . This is the purpose of the following theorem:

Theorem 2.1 We denote by ∇ x • (bW 0 ) the subspace of functions θ ∈ L 2 (Ω) such that ∃χ ∈ W 0 , θ = ∇ x • (bχ), (2.4)
and we equip it with the usual norm on L 2 (Ω). Then:

• W 0 equipped with the norm p W0 = ∇ x • (bp) L 2 is a Hilbert space, • ∇ x • (bW 0 ) is a closed subspace in L 2 (Ω),
• K is a closed subspace in L 2 (Ω).

• We have the orthogonal decomposition

L 2 (Ω) = K ⊕ K ⊥ , with K ⊥ = ∇ x • (bW 0 ). (2.5)
The demonstration of this theorem will be omitted. It can be readily adapted from that of Theorem 2.1 from [START_REF] Brull | Numerical degenerate elliptic problems and their applications to magnetized plasma simulations[END_REF]. As a consequence of this theorem, the decomposition

p ǫ = π ǫ + q ǫ , π ǫ ∈ K, q ǫ ∈ K ⊥ , (2.6)
exists and is unique for any ǫ ≥ 0. Therefore finding the particular solution p 0 which is exactly the limit of (p ǫ ) ǫ > 0 is equivalent to find π 0 and q 0 as the respective limits of (π ǫ ) ǫ > 0 and (q ǫ ) ǫ > 0 . Then, our goal is now to find some equations for π ǫ and q ǫ which are well-posed for any value of ǫ, including ǫ = 0. For this purpose, the decomposition (2.6) is introduced into (2.1), yielding

-∇ x • (H ǫ (b ⊗ b)(∇ x q ǫ -S ǫ )) + ǫλ ǫ (π ǫ + q ǫ ) = ǫf ǫ , in Ω, (H ǫ (b ⊗ b)(∇ x q ǫ -S ǫ )) • ν ≡ 0, on ∂Ω. (2.7)
The variational formulation on V writes

Ω H ǫ (b • ∇ x q ǫ )(b • ∇ x θ)dx + ǫλ ǫ Ω (π ǫ + q ǫ )θ dx = ǫ Ω f ǫ θ dx + Ω H ǫ (b • S ǫ )(b • ∇ x θ)dx, (2.8)
for any test function θ ∈ V .

In order to exhibit the equation providing

π ǫ ∈ K, the variational formulation (2.8) is tested against θ ∈ K giving Ω (λ ǫ π ǫ -f ǫ )θ dx = 0, which means that λ ǫ π ǫ -f ǫ ∈ K ⊥ for any ǫ ≥ 0. According to Theorem 2.1, there exists a function h ǫ ∈ W 0 such that λ ǫ π ǫ -f ǫ = ∇ x • (bh ǫ ).
(2.9)

This equation furnishes a means of computation for π ǫ . Firstly, applying the differential operator b • ∇ x onto (2.9) leads to an equation for h ǫ

-b • ∇ x (∇ x • (bh ǫ )) = b • ∇ x f ǫ , in Ω, (bh ǫ ) • ν ≡ 0, on ∂Ω, (2.10)
then, π ǫ is retrieved thanks to To derive an equation for q ǫ ∈ K ⊥ , we now assume that the test function θ in (2.8) is in K ⊥ . According to Theorem 2.1, there exist two functions χ and l ǫ in W 0 such that

π ǫ = 1 λ ǫ [f ǫ + ∇ x • (bh ǫ )] . ( 2
q ǫ = ∇ x • (bl ǫ ), (2.12)
and

θ = ∇ x • (bχ).
As a consequence, the variational formulation of (2.7) can be rewritten as follows:

Ω H ǫ (b • ∇ x (∇ x • (bl ǫ ))) (b • ∇ x (∇ x • (bχ))) dx + ǫλ ǫ Ω (∇ x • (bl ǫ )) (∇ x • (bχ)) dx = ǫ Ω f ǫ (∇ x • (bχ)) dx + Ω H ǫ (b • S ǫ )(b • ∇ x (∇ x • (bχ))) dx.
We recognize the variational formulation of

       b • ∇ x (∇ x • (H ǫ (b ⊗ b)∇ x (∇ x • (bl ǫ )))) -ǫλ ǫ b • ∇ x (∇ x • (bl ǫ )) = -b • ∇ x (ǫf ǫ -∇ x • (H ǫ (b ⊗ b)S ǫ )) , in Ω, (H ǫ (b ⊗ b)∇ x (∇ x • (bl ǫ ))) • ν ≡ (H ǫ (b ⊗ b)S ǫ ) • ν , on ∂Ω, (bl ǫ ) • ν ≡ 0, on ∂Ω.
(2.13) Therefore, coupling this system with (2.12), we recognize a complete definition of q ǫ which is well-posed for any ǫ ≥ 0, including ǫ = 0. Moreover this computation of q ǫ is totally compliant with the condition (2.3) and guarantees the Asymptotic-Preserving property of the scheme.

At this point, we have established a system of equations for π ǫ and q ǫ which is well-posed for any ǫ > 0 but also for ǫ = 0. Then, solving the well-posed equations (2.10), (2.13), (2.11) and (2.12) provides π 0 and q 0 as the respective limits of (π ǫ ) ǫ > 0 and (q ǫ ) ǫ > 0 when ǫ → 0. As a consequence, the sum π 0 + q 0 is exactly the solution p 0 of (2.2). Furthermore, we can remark that the limit ǫ → 0 is regular for the reformulated model (2.10)-(2.13)-(2.11)-(2.12).

Case with variable G ǫ

In this paragraph, we extend the method we have presented to the general linear case, i.e. to cases where g ǫ is of the form

g ǫ (p)(x) = G ǫ (x)p(x).
G ǫ : Ω → R is given for any ǫ ≥ 0, and is supposed to be strictly positive on Ω. In such a case, the diffusion problem (1.1) writes

-∇ x • (H ǫ (b ⊗ b)(∇ x p ǫ -S ǫ )) + ǫG ǫ p ǫ = ǫf ǫ , in Ω, (H ǫ (b ⊗ b)(∇ x p ǫ -S ǫ )) • ν ≡ 0, on ∂Ω. (2.14)
The study of these cases is motivated by the fact that the use of Gummel's algorithm on the non-linear case leads to the resolution of a sequence of linearized problems which are similar to (2.14). We refer to Section 2.2 for more details about the linearization procedure.

In order to solve the linear problem (2.14) for any value of ǫ, we use the method presented in the previous paragraph. Firstly, we define L 2 (Ω;G ǫ ) by

L 2 (Ω;G ǫ ) = p : Ω → R, p 2 L 2 (Ω;Gǫ) = Ω G ǫ (x) |p(x)| 2 dx < +∞ .
Then we introduce the following weighted Sobolev spaces:

V ǫ = p ∈ L 2 (Ω;G ǫ ) : b • ∇ x p ∈ L 2 (Ω;G ǫ ) , W ǫ = q ∈ L 2 (Ω;G ǫ ) : ∇ x • (G ǫ bq) ∈ L 2 (Ω;G ǫ ) , W 0,ǫ = {q ∈ W ǫ : (G ǫ bq) • ν ≡ 0 on ∂Ω} ,
and the set representing the functions constant along the magnetic field lines

K ǫ = {π ∈ V ǫ : b • ∇ x π = 0 on Ω} .
Following the methodology presented in the previous paragraph and in [START_REF] Brull | Numerical degenerate elliptic problems and their applications to magnetized plasma simulations[END_REF], we deduce

Corollary 2.2 W 0,ǫ equipped with the norm p W0,ǫ = ∇ x • (G ǫ bp) L 2 (Ω;Gǫ) is a Hilbert space and ∇ x • (G ǫ bW 0,ǫ ) is a closed space in L 2 (Ω;G ǫ ).
Furthermore, K ǫ is also a closed space in L 2 (Ω;G ǫ ) and we have the orthogonal decomposition

L 2 (Ω;G ǫ ) = K ǫ ⊕ K ⊥ ǫ , with K ⊥ ǫ = 1 G ǫ ∇ x • (G ǫ bW 0,ǫ ). (2.15)
From the orthogonal decomposition (2.15), the solution p ǫ of (2.14) can be uniquely decomposed as

p ǫ = π ǫ + q ǫ , π ǫ ∈ K ǫ , q ǫ ∈ K ⊥ ǫ . (2.16)
Then, if we identify the limits π 0 and q 0 of the sequences (π ǫ ) ǫ > 0 and (q ǫ ) ǫ > 0 , we will find the limit p 0 of (p ǫ ) ǫ > 0 by taking p 0 = π 0 + q 0 . In order to identify a set of equations satisfied by π ǫ and q ǫ , we follow the same procedure as in the previous paragraph: we multiply (2.14) by a test function θ ∈ V ǫ and we integrate over Ω. By choosing θ in K ǫ or in K ⊥ ǫ , we prove that π ǫ and q ǫ are respectively of the form

π ǫ = 1 G ǫ [f ǫ + ∇ x • (G ǫ bh ǫ )] , q ǫ = 1 G ǫ ∇ x • (G ǫ bl ǫ ), (2.17)
where h ǫ and l ǫ are solutions of

     -b • ∇ x 1 G ǫ ∇ x • (G ǫ bh ǫ ) = b • ∇ x f ǫ G ǫ , in Ω, (G ǫ bh ǫ ) • ν ≡ 0, on ∂Ω, (2.18)
and

                                   b • ∇ x 1 G ǫ ∇ x • H ǫ (b ⊗ b)∇ x 1 G ǫ ∇ x • (G ǫ bl ǫ ) -ǫb • ∇ x 1 G ǫ ∇ x • (G ǫ bl ǫ ) = -b • ∇ x 1 G ǫ (ǫf ǫ -∇ x • (H ǫ (b ⊗ b)S ǫ )) , in Ω, H ǫ (b ⊗ b)∇ x 1 G ǫ ∇ x • (G ǫ bl ǫ ) • ν ≡ (H ǫ (b ⊗ b)S ǫ ) • ν , on ∂Ω, (G ǫ bl ǫ ) • ν ≡ 0, on ∂Ω. (2.19)
As in the previous paragraph, we observe that the equations (2.17)-(2.18)-(2.19) remain well-posed for any ǫ ≥ 0. As a consequence, the particular solution p 0 of the limit problem we are looking for is exactly the sum π 0 + q 0 where π 0 and q 0 are computed by solving (2.17)-(2.18)-(2.19) with ǫ = 0. Furthermore, the resolution of the fourth order problem (2.19) can be replaced by the successive resolution of two homogeneous Dirichlet type problems which are

     -b • ∇ x 1 G ǫ ∇ x • (H ǫ bL ǫ ) + ǫL ǫ = -ǫb • ∇ x f ǫ G ǫ -S ǫ , in Ω, (H ǫ bL ǫ ) • ν ≡ 0, on ∂Ω, (2.20)
and

     -b • ∇ x 1 G ǫ ∇ x • (G ǫ bl ǫ ) = L ǫ -b • S ǫ , in Ω, (G ǫ bl ǫ ) • ν ≡ 0, on ∂Ω.
(2.21)

AP-scheme derivation for non-linear problems

Finally, we consider the general model (1.1) given in the introduction when the function p → g ǫ (p) is non-linear. When ǫ goes to 0, the model becomes

-∇ x • (H 0 (b ⊗ b)(∇ x p0 -S 0 )) = 0, in Ω, (H 0 (b ⊗ b)(∇ x p0 -S 0 )) • ν ≡ 0, on ∂Ω. (2.22)
Due to the non-linearity of the function p → g ǫ (p) the orthogonal decomposition method cannot be used. Then we choose to linearize the diffusion equation (1.1) by using Gummel's algorithm developed in [START_REF] Gummel | A self-consistent iterative scheme for one-dimensional steady state transistor calculations[END_REF]. This iterative method consists in the approximation of the solution p ǫ by a sequence (p ǫ,N ) N ≥ 0 defined by

p ǫ,N +1 = p ǫ,N + δ ǫ,N , (2.23)
and initialized with an arbitrary p ǫ,0 . In this method, each δ ǫ,N is viewed as a small correction of p ǫ,N in order to obtain p ǫ,N +1 . Then, assuming that p ǫ,N +1 is a solution of (1.1), it holds that

               -∇ x • H ǫ (b ⊗ b) ∇ x p ǫ,N + ∇ x δ ǫ,N -S ǫ ǫ +g ǫ (p ǫ,N ) + δ ǫ,N g ′ ǫ (p ǫ,N ) + O(δ 2 ǫ,N ) = f ǫ , in Ω, H ǫ (b ⊗ b) ∇ x p ǫ,N + ∇ x δ ǫ,N -S ǫ ǫ • ν ≡ 0, on ∂Ω.
(2.24)

Then, neglecting second order terms in δ ǫ,N , we obtain a linear diffusion problem for δ ǫ,N which writes

-∇ x • (H ǫ (b ⊗ b)(∇ x δ ǫ,N -S ǫ,N )) + ǫG ǫ,N δ ǫ,N = ǫf ǫ,N , in Ω, (H ǫ (b ⊗ b)(∇ x δ ǫ,N -S ǫ,N )) • ν ≡ 0, on ∂Ω, (2.25)
where G ǫ,N , f ǫ,N and S ǫ,N are defined by

G ǫ,N = g ′ ǫ (p ǫ,N ), f ǫ,N = f ǫ -g ǫ (p ǫ,N ), S ǫ,N = S ǫ -∇ x p ǫ,N .
For each value of N , the problem (2.25) is of the same kind as (2.14). So we can solve it by applying the method described in the paragraph 2.1.2. This sequence of linearized problems can also be obtained from Newton's iterative method to solve

F ǫ (p ǫ ) = 0, (2.26)
where the differential operator F ǫ is defined as

F ǫ (p) = -∇ x • H ǫ (b ⊗ b) ∇ x p -S ǫ ǫ + g ǫ (p) -f ǫ .
Indeed, Newton's method for solving (2.26) writes

DF ǫ (p ǫ,N )(p ǫ,N +1 -p ǫ,N ) = -F ǫ (p ǫ,N ),
where DF ǫ (p) is the derivative in p of the differential operator F ǫ (p) and is of the form

DF ǫ (p)(δ) = -∇ x • H ǫ (b ⊗ b) ∇ x δ ǫ + g ′ ǫ (p) × δ .

Numerical method

In this section, we present a numerical method which allows to solve the diffusion problems (2.14) and (1.1) by using the decomposition approaches we have presented. First, we introduce some notations which will be used for the construction of the scheme, then we present the scheme itself for the general linear case (2.14). Finally, we present the discretized version of Gummel's algorithm for the non-linear case.

Notations and definitions

We consider a uniform mesh (x i ,y j ) defined by

x i = x min + i∆x, y j = y min + j ∆y , ∆x = x max -x min N x + 1 , ∆y = y max -y min N y + 1 ,
and we assume that the simulation domain is

Ω = [x -1/2 ,x Nx+1/2 ] × [y -1/2 ,y Ny+1/2
]. We also consider the following subsets of Z 2 : I = {0,...,N x } × {0,...,N y },

I = {-1,...,N x + 1} × {-1,...,N y + 1}, I * = {0,...,N x -1} × {0,...,N y -1}, I * = {-1,...,N x } × {-1,...,N y },
and we consider the notation h = max(∆x,∆y).

(3.1)

Since the decomposition method we have presented in paragraph 2.2 is based on variational formulations of the diffusion problem for p ǫ and uses the duality between the operators p → b • ∇ x p and p → ∇ x • (bp), we choose to approach these differential operators by ∂ h and ∂ h, * respectively such that the duality property is preserved at the discrete level. For this purpose, we define ∂ h and ∂ h, * such that

(∂ h θ) i+1/2,j+1/2 = b i+1/2,j+1/2 •     θ i+1,j+1 -θ i,j+1 + θ i+1,j -θ i,j 2∆x θ i+1,j+1 -θ i+1,j + θ i,j+1 -θ i,j 2∆y
    , (3.2)
for all θ = (θ i,j ) (i,j) ∈ I , and

(∂ h, * χ) i,j = α∈{±1} (b x χ) i+1/2,j+α/2 -(b x χ) i-1/2,j+α/2 2∆x + α∈{±1} (b y χ) i+α/2,j+1/2 -(b y χ) i+α/2,j-1/2 2∆y , (3.3)
for all χ = (χ i+1/2,j+1/2 ) (i,j) ∈ I * .

Linear problems

We assume that the function p → g ǫ (p) is given by

g ǫ (p)(x,y) = G ǫ (x,y)p(x,y),
where G ǫ is analytically known. We also assume that the functions b, H ǫ , f ǫ and S ǫ are analytically known and we consider the following notations:

f ǫ,i,j = f ǫ (x i ,y j ) , G ǫ,i,j = G ǫ (x i ,y j ) , G ǫ,i+1/2,j+1/2 = G ǫ (x i+1/2 ,y j+1/2 ) , H ǫ,i+1/2,j+1/2 = H ǫ (x i+1/2 ,y j+1/2 ) , S ǫ,i+1/2,j+1/2 = S ǫ (x i+1/2 ,y j+1/2 ) , b i+1/2,j+1/2 = b(x i+1/2 ,y j+1/2 ).
Then, the diffusion problem (2.14) can be discretized under the following form:

(-∂ h, * (H ǫ (∂ h p ǫ,app -b • S ǫ ) + ǫG ǫ p ǫ,app ) i,j = ǫf ǫ,i,j , ∀(i,j) ∈ I , (H ǫ (∂ h p ǫ,app -b • S ǫ )(b • ν)) i+1/2,j+1/2 = 0, ∀(i,j) ∈ I * \I * , (3.4)
the approximation of p ǫ,app of p ǫ is computed at the points (x i ,y j ) ∈ Ω. Since ∂ h and ∂ h, * have been chosen to be dual operators, we follow the decomposition approach we have presented in Section 2.1.2 at a discrete level by using some discrete variational formulations of (3.4). Writing p ǫ,app,i,j = π ǫ,app,i,j + q ǫ,app,i,j , with π ǫ,app satisfying (∂ h π ǫ,app ) i+1/2,j+1/2 = 0, ∀(i,j) ∈ I * , π ǫ,app and q ǫ,app are completely defined by

π ǫ,app,i,j = 1 G ǫ,i,j f ǫ,i,j + (∂ h, * (G ǫ h ǫ,app )) i,j , (3.5)
q ǫ,app,i,j = 1 G ǫ,i,j (∂ h, * (G ǫ l ǫ,app )) i,j ,
where h ǫ,app = (h ǫ,app,i+1/2,j+1/2 ) (i,j) ∈ I * and l ǫ,app = (l ǫ,app,i+1/2,j+1/2 ) (i,j) ∈ I * are computed by inverting the following systems:

              -∂ h 1 G ǫ ∂ h, * (G ǫ h ǫ,app ) i+1/2,j+1/2 = ∂ h f ǫ G ǫ i+1/2,j+1/2 , ∀(i,j) ∈ I * , h ǫ,app,i+1/2,j+1/2 = 0, ∀(i,j) ∈ I * \I * , (3.6)
               -∂ h 1 G ǫ ∂ h, * (H ǫ L ǫ,app ) + ǫL ǫ,app i+1/2,j+1/2 = -ǫ ∂ h f ǫ G ǫ -b • S ǫ i+1/2,j+1/2
, ∀(i,j) ∈ I * , L ǫ,app,i+1/2,j+1/2 = 0, ∀(i,j) ∈ I * \I * ,

and

         -∂ h 1 G ǫ ∂ h, * (G ǫ l ǫ,app ) i+1/2,j+1/2 = (L ǫ,app -b • S ǫ ) i+1/2,j+1/2 , ∀(i,j) ∈ I * , l ǫ,app,i+1/2,j+1/2 = 0, ∀(i,j) ∈ I * \I * .
(3.8)

Non-linear problems

In this paragraph, we detail the discretized version of Gummel's algorithm presented in Section 2.2. In order to initialize the loop, we compute the following initial datas:

p ǫ,0,i,j = p ǫ,0 (x i ,y j ), ∀(i,j) ∈ I , S ǫ,i+1/2,j+1/2 = S ǫ (x i+1/2 ,y j+1/2 ), ∀(i,j) ∈ I * , f ǫ,i,j = f ǫ (x i ,y j ), ∀(i,j) ∈ I , H ǫ,i+1/2,j+1/2 = H ǫ (x i+1/2 ,y j+1/2 ), ∀(i,j) ∈ I * .
Then, the N -th iteration of Gummel's algorithm is set as follows:

• Step 1: assuming that p ǫ,N,i,j , ∀(i,j) ∈ I , are known, we have

(b • S ǫ,N ) i+1/2,j+1/2 = (b • S ǫ -∂ h p ǫ,N ) i+1/2,j+1/2 , ∀(i,j) ∈ I * , f ǫ,N,i,j = f ǫ,i,j -g ǫ (p ǫ,N,i,j ), ∀(i,j) ∈ I , G ǫ,N,i,j = g ′ ǫ (p ǫ,N,i,j ), ∀(i,j) ∈ I , G ǫ,N,i+1/2,j+1/2 = g ′ ǫ (p ǫ,N,i+1/2,j+1/2 ), ∀(i,j) ∈ I * ,
where p ǫ,N,i+1/2,j+1/2 = 1 4 (p ǫ,N,i+1,j+1 + p ǫ,N,i+1,j + p ǫ,N,i,j+1 + p ǫ,N,i,j ).

•

Step 2: we compute h ǫ,N,i+1/2,j+1/2 and l ǫ,N,i+1/2,j+1/2 for all (i,j) ∈ I * by solving

               -∂ h 1 G ǫ,N ∂ h, * (G ǫ,N h ǫ,N ) i+1/2,j+1/2 = ∂ h f ǫ,N G ǫ,N i+1/2,j+1/2 , ∀(i,j) ∈ I * , h ǫ,N,i+1/2,j+1/2 = 0, ∀(i,j) ∈ I * \I * ,                -∂ h 1 G ǫ,N ∂ h, * (H ǫ L ǫ,N ) + ǫL ǫ,N i+1/2,j+1/2 = -ǫ ∂ h f ǫ,N G ǫ,N -b • S ǫ,N i+1/2,j+1/2
, ∀(i,j) ∈ I * , L ǫ,N,i+1/2,j+1/2 = 0, ∀(i,j) ∈ I * \I * , and

         -∂ h 1 G ǫ,N ∂ h, * (G ǫ,N l ǫ,N ) i+1/2,j+1/2 = (L ǫ,N -b • S ǫ,N ) i+1/2,j+1/2 , ∀(i,j) ∈ I * , l ǫ,N,i+1/2,j+1/2 = 0, ∀(i,j) ∈ I * \I * .
• Step 3: we compute δ ǫ,N,i,j for all (i,j) ∈ I by using

δ ǫ,N,i,j = 1 G ǫ,N,i,j f ǫ,N,i,j + (∂ h, * (G ǫ,N (h ǫ,N + l ǫ,N ))) i,j ,
and we obtain p ǫ,N +1,i,j for all (i,j) ∈ I.

• Step 4: we compute p ǫ,N +1,i,j for all (i,j) ∈ I\I by using the boundary condition

(∂ h p ǫ,N +1 ) i+1/2,j+1/2 -(b • S ǫ ) i+1/2,j+1/2 = 0, ∀(i,j) ∈ I * \I * .

Numerical investigations of the AP-scheme

This section is devoted to numerical investigations of the Asymptotic-Preserving scheme derived in Sections 2 and 3. The validation procedure consists in manufacturing p ǫ , an analytic solution of the model problem (1.1) which is compared to the numerical approximation p ǫ,app carried out thanks to the AP-scheme. These experiments are performed in two dimensions using a uniform Cartesian mesh independent of the anisotropy direction. For simplicity purpose, the first numerical experiments are performed in the framework on the linear model, but the conclusions drawn from these investigations apply to the general non-linear problem.

Numerical convergence of the scheme

The first numerical tests aim at demonstrating the convergence of the APscheme regardless to the asymptotic parameter values. With this aim, an analytic solution is manufactured for the problem (1.1) in the linear case, i.e. with g ǫ (p) = G ǫ (x)p(x). First, the expression for the anisotropy direction b and the functions G ǫ and H ǫ are defined on

Ω = [1,2] × [1,2] thanks to G ǫ (x,y) = H ǫ (x,y) = 1 + sin 2 (x) sin 2 (y), (4.1) b = (sinθ,-cosθ) with θ(x,y) = arctan y x . (4.2)
Then the expression of p ǫ , as defined by

p ǫ (x,y) = 1 1 + x 2 + y 2 ,
is used to analytically compute f ǫ and S ǫ with

f ǫ = G ǫ p ǫ , S ǫ = ∇ x p ǫ . (4.3)
These definitions are inserted in the numerical method described in Section 3.2 to compute the numerical approximation p ǫ,app finally compared to the exact solution p ǫ . The relative errors denoted e p , p ∈ {2;∞}, are defined by

e 2 = p ǫ -p ǫ,app ℓ 2 (I) p ǫ ℓ 2 (I) , e ∞ = p ǫ -p ǫ,app ℓ ∞ (I) p ǫ ℓ ∞ (I) .
These quantities are displayed on Figure 1(a) as functions of the space step h and for different anisotropy strengths ǫ = 10 -1 , ǫ = 10 -9 and ǫ = 0. A linear decrease of the errors is observed with the mesh refinement, the slope being equal to 2, which is consistent with the definitions (3. including ǫ = 0. This demonstrates the ǫ-invariance of the numerical scheme second order accuracy with respect to the space step h.

The ability of the scheme to compute a solution component π ǫ with no gradient in the anisotropy direction is also investigated. The numerical approximation of π ǫ , π ǫ,app provided by (3.4)- (3.5), should verify a discrete analogous of the property b • ∇ x π ǫ = 0. This is analyzed thanks to Figure 1(b), where the evolution of ∂ h π ǫ,app ℓ p (I * ) / p ǫ ℓ p (I) as a function of the space step is displayed for p = 2,∞, ǫ = 10 -1 , ǫ = 10 -9 and ǫ = 0. Note that the quantity ∂ h π ǫ,app is the residual of the linear system solved to compute the solution of (3.4), and consequently characterizes the precision of the linear system solver. For these test cases, a sparse direct solver being used [START_REF] Amestoy | A fully asynchronous multifrontal solver using distributed dynamic scheduling[END_REF], the accuracy is very close to the computer arithmetic precision, at least for small linear system sizes. This precision is observed to deteriorate moderately with the increase of the system size which explains the growth of the error with vanishing mesh sizes. However this does not affect the precision of the scheme, as demonstrated by the results of Figure 1(a).

Anisotropy angle influence on the method accuracy

In this section, we quantify the sensitivity of the numerical method with respect to the anisotropy direction variations. More precisely, we wish to analyze the accuracy of the method as a function of α, the angle measured between the anisotropy direction and the first direction (associated to the first coordinate). The anisotropy direction b is assumed to be uniform and defined as

b = (sinα,-cosα), α ∈ [0, π 2 ].
In order to manufacture an analytic solution for the problem, we introduce a system of coordinates which is adapted to b. These coordinates are denoted (X,Y ) and are deduced from (x,y) by the relations

X = x cosα + y sinα, Y = x sinα -y cosα. (4.4)
In these coordinates, the linear diffusion problem (2.14) writes

-∂ Y (H ǫ ∂ Y P ǫ ) + ǫG ǫ P ǫ = ǫF ǫ -∂ Y (H ǫ B • S ǫ ), in Ω, ∂ Y P ǫ = B • S ǫ , on ∂Ω, (4.5)
with

P ǫ (X,Y ) = p ǫ (x,y), H ǫ (X,Y ) = H ǫ (x,y), G ǫ (X,Y ) = G ǫ (x,y), F ǫ (X,Y ) = f ǫ (x,y), S ǫ (X,Y ) = S ǫ (x,y) and B(X,Y ) = b(x,y).
It is straightforward to verify that the function P ǫ given by

P ǫ (X,Y ) = sin(X) + 1 G ǫ (X,Y ) ∂ Y (G ǫ L ǫ ),
is the solution of (4.5) provided that F ǫ and S ǫ satisfy

F ǫ = G ǫ P ǫ , B • S ǫ = ∂ Y P ǫ ,
and where L ǫ (X,Y ) = l ǫ (x,y) with l ǫ| ∂Ω = 0. This requirement is met by the following definition

l ǫ (x,y) = sin 2π (x -x -1/2 ) x Nx+1/2 -x -1/2 sin 2π (y -y -1/2 ) y Ny+1/2 -y -1/2 ,
which ensures l ǫ (x i+1/2 ,y j+1/2 ) = 0 for any (i,j) ∈ I * \I * . The problem is stated in Cartesian coordinates thanks to the change of variables (4.4) yielding to p ǫ (x,y) = π ǫ (x,y) + q ǫ (x,y) with

π ǫ (x,y) = sin(x cosα + y sinα) , q ǫ (x,y) = 1 G ǫ (x,y) ∇ x • (G ǫ (x,y)b(x,y)l ǫ (x,y)) ,
the other coefficients being manufactured similarly with G ǫ and H ǫ given by (4.1).

In the following tests, the computation domain Ω = [1,2] × [1,2] is discretized thanks to a uniform mesh constituted of 200 × 200 cells. The relative approximation error as a function of the angle α is displayed on Figure 2 for different norms. The numerical method accuracy is observed to remain almost unaltered by the anisotropy direction changes. More precisely, we observe a variation of the relative errors in norms ℓ 1 and ℓ 2 lower than 4% and a variation of ℓ ∞ lower than 7%. Furthermore, these observations are redundant for several values of ǫ: in Figure 2, we have considered ǫ = 10 -3 and ǫ = 10 -8 and the obtained error curves are very close. Note that other experiments have been carried out for anisotropy strengths ranging from ǫ = 0 to ǫ = 1 and with other definitions of G ǫ and of H ǫ , with comparable results. The curves being very similar to that of Figure 2, these plots are omitted. These observations confirm one of the main ideas of the present paper: the accuracy of the method is almost independent of the anisotropy direction relatively to the grid, i.e. the mesh over Ω can be constructed whatever the anisotropy direction and strength, without a significant loss of accuracy.

Convergence of Gummel's loop

The third test sequence is devoted to the convergence of the linear problems sequence defined in Sections 2.2 and 3.3 for solving the non-linear model (1.1). The process detailed in the preceding sections is again implemented to manufacture an analytic solution for the non-linear problem.

The computational domain remains Ω = [1,2] × [1,2], the anisotropy direction b is a function of the space variables whose expression is given by equation (4.2), H ǫ and g ǫ being defined as

H ǫ (x,y) = 1 + cos 2 (x) cos 2 (y), g ǫ (p) = p 6 . (4.6)
Note that this choice of g ǫ (p) introduces a severe non-linearity in the problem. Several tests have also been performed with other definitions of g ǫ (p), for instance

g ǫ (p) = -p(1 -p) p - 1 2 ,
which defines an anisotropic diffusion-reaction equation similar to the steadystate Allen-Cahn equation (see [START_REF] Beněs | Computational studies of non-local anisotropic Allen-Cahn equation[END_REF][START_REF] Elliot | The limit of the fully anisotropic doubleobstacle Allen-Cahn equation in the non-smooth case[END_REF]) used in phase transition problems. These tests produce results almost identical to the results which are obtained when g ǫ (p) = p 6 so we only consider the strongly non-linear reaction term defined in (4.6) within the presentation of the numerical results in the next lines.

The solution p ǫ is constructed thanks to a cubic spline S, precisely

p ǫ (x,y) = 1 + S x -x mid L x S y -y mid L y , (4.7)
with S(z) = 0 for |z| / ∈ [0,2] and

S(z) = 1 6 (2 -|z| 3 ), if 1 ≤ |z| ≤ 2, 2 3 -|z| 2 + 1 2 |z| 3 , if 0 ≤ |z| < 1, (4.8)
with (x mid ,y mid ) = ( 3 2 , 3 2 ) and L x = L y = 1/10. To analyze the convergence with respect to the number of Gummel's iterations, the sequence is initiated with p ǫ,0 , a perturbation of the non-linear problem solution, reading p ǫ,0 (x,y) = p ǫ (x,y) + η max 0,1 -µ(x -x mid ) 2 -µ(y -y mid ) 2 , (4.9)

where µ and η µ are parameters controlling the support and the magnitude of the perturbation. Since Gummel's method is constructed on a linearization of the problem its convergence cannot be guaranteed with a poor estimation of the solution as initial guess. It means that the parameters µ and η cannot be chosen completely arbitrarily: indeed, several simulations have been performed, all with the same parameters except η ranging in {0,10,20,40,60,100,1000} and µ ranging in {1,10,60,100} and it has been observed that Gummel's method does not converge as N → ∞ when η is larger than 10 2 . Concerning the parameter µ, the simulation sequence reveals that the convergence of Gummel's method is almost not affected by the amplitude of µ.

The successive relative errors measured between the iterates of the Gummel's loop and the exact solution are plotted on Figure 3(a). The computations are carried out on two different meshes, M 100 and M 1000 with 100 × 100 and 1000 × 1000 cells, with η = 0.1 and µ = 60 and for anisotropy strengths including ǫ = 0. Along with the graphical representation of the solution approximation error, the evolution of the corrector norm relative to that of the solution, namely the quantity δ ǫ,N,app ℓ 2 (I) / p ǫ ℓ 2 (I) , is also plotted in Figure 3(b). These last results being almost identical for both meshes, the plot related to the finest mesh is omitted in this figure .  In spite of the large perturbation amplitude, Gummel's iterative method converges in a small number of iterations, for both meshes and for all ǫ-values. The corrector term δ ǫ,N,app rapidly decreases to reach the computer precision threshold (10 -15 ) after 4 iterations. In the same time, the relative error also decreases but the approximation is not improved by subsequent iterations, the error remaining constant for iteration numbers greater than 4. At this stage, the precision of the approximation is not limited by the linearization process of the Gummel's loop anymore, but by the discretization error of the linearized problem, explaining the plateau described by the error. To document this analyzis further, we summarize in Table 1 the values of the relative error measured 1: Relative error E p,N f = p ǫ -p ǫ,N f ,app ℓ p (I) / p ǫ ℓ p (I) (p = 1,2,∞) for the non-linear problem defined by g ǫ (p) = p 6 . The computations are carried out on uniform meshes M k constituted of k × k cells (k = 100,200,500,1000) with several values of ǫ and after a number of iteration of Gummel's loop N f large enough for the convergence to be effective. between the exact solution and the approximation obtained after N f iterations of the Gummel's loop. This quantity is referred to as E p,N f (p = 1,2,∞) and computed for N f large enough to ensure that the plateau above mentioned is reached. For the investigations carried out, this requirement is met as soon as N f ≥ 4. The approximation error E p,N f is observed to quadratically decrease with the space mesh: the error norms related to the computations performed on a 1000 × 1000 mesh are for instance 10 2 times as small as those carried out on a mesh with 100 × 100 cells. This is a consequence of the second order accurate discretization of the spatial operator already outlined in section 4.1. Finally, the results of Table 1 also demonstrate the independence of the numerical method precision with respect to the anisotropy intensity.

Highlight of the scheme Asymptotic-Preserving property

These last experiments are devoted to illustrate the Asymptotic-Preserving property of the numerical method, i.e. its ability to compute an accurate approximation of p 0 , the solution of the limit problem (2.2). The solution of the problem is constructed as a sequence (p ǫ ) ǫ > 0 defined by The functions g ǫ , H ǫ and b are defined as in the previous test sequence, the initial guess for Gummel's loop being constructed following (4.9) using the same perturbation. We now wish to evaluate the error measured between the exact solution of the limit problem p 0 and the approximation computed thanks to the AP-scheme for vanishing ǫ. This error, denoted E ǫ and defined as E ǫ = p ǫ,app -p 0 ℓ 2 (I) / p 0 ℓ 2 (I) , is plotted on Figure 4(a) as a function of ǫ. The data represented on this figure are obtained after convergence of the Gummel's loop. Two regimes can be identified. The first one is related to the largest values of ǫ for which a linear decrease of the error is observed. The second one is a plateau whose value depends on the mesh step h, this value being lower for refined meshes. Precisely we note a quadratic decrease of this value with the mesh size. To explain these features, we use the following identity p ǫ,app -p 0 = p ǫ,app -p 0,app + p 0,app -p 0 . This yields E ǫ ≤ E ǫ,app + e 0 where e 0 = p 0,app -p 0 ℓ 2 (I) / p 0 ℓ 2 (I) represents the approximation error of p 0 , p 0,app being the numerical approximation of p 0 provided by the AP-scheme with ǫ = 0, and E ǫ,app = p ǫ,app -p 0,app ℓ 2 (I) / p 0 ℓ 2 (I) .

p ǫ = p 0 + ǫ p1 ǫ , with p 0 (x,y) = 1 + S x -x mid L x S y -y mid L y , ( 4
The error E ǫ linearly decreases with ǫ as long as the approximation error e 0 is negligible compared to E ǫ,app (see Figure 4(b)). Below a given ǫ-value, varying with the mesh size, the total error can be assimilated to e 0 and the decrease of ǫ is ineffective. The discrete operators being second order accurate e 0 is quadratically decreasing with the mesh step h.

As a consequence, we can conclude that p ǫ,app converges to p 0 when ǫ converges to 0 alongwith h. This is exactly the Asymptotic-Preserving property of the scheme we intended to validate.

Conclusions and perspectives

In this paper we have presented an Asymptotic-Preserving numerical method for singular perturbation of non-linear anisotropic reaction-diffusion problems. The Asymptotic-Preserving property of the scheme is ensured thanks to a solution decomposition explained in full details in the most simple framework of a linear problem. This method is then generalized to non-linear problems thanks to Gummel's linearization method.

In a second part, several two-dimensional numerical investigations of the AP-scheme are performed. These tests reveal a very weak dependence of the scheme accuracy with respect to the anisotropy direction, demonstrating the relevance of the use of non-adapted coordinates. The Asymptotic-Preserving property of the scheme is also validated for vanishing ǫ on linear as well as non-linear problems. The solution of the limit problem is accurately captured with no restrictions on the anisotropy strength. Furthermore, the computational efficiency of the method, in terms of memory as well as CPU usage, does not depend on this anisotropy strength.

Several applications of the present work can be investigated: at present time, the method has been used for the resolution of linear anisotropic diffusion problems for a two-fluid Euler-Lorentz model (see [START_REF] Brull | Asymptotic-Preserving scheme for a two-fluid Euler-Lorentz model[END_REF]) and the non-linear version of the method will be coupled to an Asymptotic-Preserving scheme for a onefluid full Euler-Lorentz model (see [START_REF] Brull | A numerical investigation of the full Euler-Lorentz model with a large magnetic field[END_REF]).

. 11 )

 11 Note that the system (2.10)-(2.11) is well-posed and does not degenerate for any value of ǫ ≥ 0, including ǫ = 0. It provides a means of computing the macro component of the solution regardless to ǫ values.

  2) and (3.3) of ∂ h and ∂ h, * as second order accurate approximations of the differential operators b • ∇ x and ∇ x • (b•). Furthermore, this property holds for all considered values of ǫ, ep = pǫ -pǫ,app ℓ p (I) / pǫ ℓ p (I) . ep = ∂ h πǫ,app ℓ p (I * ) / pǫ ℓ p (I) .

Figure 1 :

 1 Figure 1: Relative error p ǫ -p ǫ,app ℓ p (I) / p ǫ ℓ p (I) (left) and parallel gradient ∂ h π ǫ,app ℓ p (I * ) / p ǫ ℓ p (I) (right) in ℓ 2 and ℓ ∞ norms as functions of the space step h for non-uniform data G ǫ , H ǫ and b defined by (4.1)-(4.2), and different anisotropy strengths ǫ = 1, ǫ = 10 -9 and ǫ = 0.

  pǫ -pǫ,app ℓ p (I) / pǫ ℓ p (I) , ǫ = 10 -3 . pǫ -pǫ,app ℓ p (I) / pǫ ℓ p (I) , ǫ = 10 -8 .

Figure 2 :

 2 Figure 2: Relative error p ǫ -p ǫ,app ℓ p (I) / p ǫ ℓ p (I) (p = 1,2,∞) as a function of α: case with G ǫ (x,y) = H ǫ (x,y) = 1 + sin 2 (x) sin 2 (y) and ǫ = 10 -3 (left) and ǫ = 10 -8 (right).

  E 2,N , M 100 and M 1000 meshes. δ ǫ,N,app ℓ 2 (I) / pǫ ℓ 2 (I) , M 100 mesh.

Figure 3 :

 3 Figure 3:Gummel's iteration convergence: evolution in log 10 -scales of the relative error E 2,N = p ǫ -p ǫ,N f ,app ℓ 2 (I) / p ǫ ℓ 2 (I) (left) and of δ ǫ,N,app ℓ 2 (I) / p ǫ ℓ 2 (I) (right) as functions of th iteration number N for the non-linear problem. The simulations are performed with the uniform meshes M 100 and M 1000 composed of 100 × 100 and 1000 × 1000 cells (for the corrector norm, the values being very similar for both meshes, only those of the coarsest one are displayed).

2 -

 2 norm of the error p ε,app -p 0 in log 10 -scale p Eǫ, M 200 and M 1000 meshes.

2 -

 2 norm of the error p ε,app -p 0,app in log 10 -scale p Eǫ,app, M 200 and M 1000 meshes.

Figure 4 :

 4 Figure 4: Evolution in log 10 -scales of E ǫ = p ǫ,app -p 0 ℓ 2 (I) / p 0 ℓ 2 (I) (left) and of E ǫ,app = p ǫ,app -p 0,app ℓ 2 (I) / p 0 ℓ 2 (I) (right) as functions of ǫ computed on uniform meshes constituted of 200 × 200 cells (M 200 ) and 1000 × 1000 cells (M 1000 ) are considered. The simulations are performed with b, H ǫ , g ǫ defined by (4.6) and (4.2).

  10 -5 9.8116 × 10 -6 1.5673 × 10 -6 3.9166 × 10 -7 10 -1 E 2,N f 1.0446 × 10 -4 2.6188 × 10 -5 4.1988 × 10 -6 1.0505 × 10 -6 E ∞,N f 6.0730 × 10 -4 1.5793 × 10 -4 2.5942 × 10 -5 6.5451 × 10 -6 E 1,N f 3.9796 × 10 -5 9.8969 × 10 -6 1.5808 × 10 -6 3.9504 × 10 -7 10 -12 E 2,N f 1.0496 × 10 -4 2.6311 × 10 -5 4.2184 × 10 -6 1.0554 × 10 -6 E ∞,N f 6.1098 × 10 -4 1.5885 × 10 -4 2.6087 × 10 -5 6.5815 × 10 -6 E 1,N f 3.9796 × 10 -5 9.8969 × 10 -6 1.5808 × 10 -6 3.9504 × 10 -7 0 E 2,N f 1.0496 × 10 -4 2.6311 × 10 -5 4.2184 × 10 -6 1.0554 × 10 -6 E ∞,N f 6.1098 × 10 -4 1.5885 × 10 -4 2.6087 × 10 -5 6.5815 × 10 -6

	ǫ	Error	M 100	M 200	M 500	M 1000
	E 1,N f 3.9452 × Table

Acknowledgement. This work has been supported by the french magnetic fusion programme FR-FCM, by the INRIA large-scale initiative 'FUSION', by 'BOOST' and 'IODISSEE' ANR projects and by the CEA-Cadarache in the frame of the contract 'APPLA' (# V3629.001 av. 2). The authors wish to thank P. Degond for suggesting this problem and for very fruitful discussions on the topic.