Ablative and catalytic behavior of carbon-based porous thermal protection materials in nitrogen plasmas

Accepted Manuscript Ablative and catalytic behavior of carbon-based porous thermal protection materials in nitrogen plasmas Gerard L. Vignoles, Alessandro Turchi, Daniele Bianchi, Pierre Blaineau, Xavier Lamboley, Damien Le Quang Huy, Cyril Levet, Olivier Caty, Olivier Chazot PII:

S0008-6223(18)30338-5

DOI:

10.1016/j.carbon.2018.03.087

Reference:

CARBON 13030

To appear in:

Carbon

Received Date: 18 July 2017 Revised Date:

27 March 2018

Accepted Date: 28 March 2018

Please cite this article as: G.L. Vignoles, A. Turchi, D. Bianchi, P. Blaineau, X. Lamboley, D. Le Quang Huy, C. Levet, O. Caty, O. Chazot, Ablative and catalytic behavior of carbon-based porous thermal protection materials in nitrogen plasmas, Carbon (2018), doi: 10.1016/j.carbon.2018.03.087. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Rebuilt flow conditions Measured recession rate

Fiber diameter (µm)

11

2

Affected depth δa

M AN U

10 9

TE D

8

Ablation rate (GS difference x frame rate )

4

SC

12

Avg. grayscale

0

Concentration (walkers/slice)

800

7 6

EP

0

200

Distance from surface (µm) 400

AC C

Press. (Pa) 1700 1665 1630 Temp. (K) 5500 4200 2900 1600

New analytical and numerical models

RI PT

Nitrogen plasma ablation testing on carbon preforms

600

800

Affected depth observation & measurement Effective rate constants

400

1000 0

0

400

1200 800 1600 Z Position in image (µm)

Intrinsic rate constants

2000

ACCEPTED MANUSCRIPT

Ablative and catalytic behavior of carbon-based porous thermal protection materials in nitrogen plasmas

RI PT

Gerard L. Vignolesa,∗, Alessandro Turchib , Daniele Bianchic , Pierre Blaineaua , Xavier Lamboleya,d , Damien Le Quang Huyb , Cyril Leveta,d , Olivier Catya , Olivier Chazotb a

M AN U

SC

University of Bordeaux, Laboratoire des Composites ThermoStructuraux (LCTS) UMR 5801: CNRS-Safran-CEA-Univ. Bordeaux 3, All´ee de La Bo´etie, 33600 Pessac, France Tel: (+33) 5 5684 4708 b von Karman Institute for Fluid Dynamics, Chauss´ee de Waterloo 72, B-1640 Rhode-Saint-Gen`ese, Belgium c Sapienza University of Rome, Department of Mechanical and Aerospace Engineering, Via Eudossiana 18, 00184 Rome, Italy d CEA,CESTA, 15 Avenue des Sabli`eres, 33114 Le Barp, France

Abstract

TE D

Heat shields used to protect space capsules during very high-speed atmospheric entry incorporate lightweight insulating refractories based on carbon-fiber preforms. These ablators, with up to 80% porosity, present exceptional thermal and chemical properties. A joint experimental and modeling approach to study such

EP

materials is presented, which contributes to improving their design in relevant operating conditions. Samples were tested in an inductively coupled plasma at 1.5

AC C

kPa and 25 kPa pressures, with surface temperatures ranging from 1500 to 2800 K. The recession rate was measured in-situ and the local flow conditions were reconstructed numerically from experimental data. The porous medium was imaged by X-ray Computerized Micro-Tomography (CMT), and the depth affected by the gas-solid interaction phenomena due to the plasma exposure was extracted from these data. A simple analytical model has been derived to relate observable ∗

Corresponding author. E-mail: [email protected]

Preprint submitted to Elsevier

April 5, 2018

ACCEPTED MANUSCRIPT

and reconstructed quantities and the fiber-scale reaction constants. Image-based numerical simulations of ablation by nitridation, with simultaneous catalytic re-

RI PT

combination of atomic nitrogen, were compared to the analytical model and used to extract the intrinsic, fiber-scale reaction rate constants from test data. Results show, among others, that the fiber-scale nitridation rate constant is close to the rate of graphite nitridation, and that it decreases strongly with pressure.

SC

Keywords: Ablation; nitridation; carbon fibers properties; modeling ; fibrous media; plasma test

M AN U

Highlights

• Affected depth was retrieved from X-ray CMT scans of plasma-tested fibrous carbon preforms by an original image processing technique. • Model and 3D numerical simulations of material recession linkrelate af-

TE D

fected depth to diffusion/reaction ratios.

• Fiber-scale recessionnitridation and catalytic recombination reactions rate constants were identified from experimental results, analytical model and

EP

3D numerical simulations.

• Fiber-scale nitridation rate constants are close to graphite values.

AC C

• The identified recombination rate constants were found higher than the nitridation rate constants.

• The effective (i.e. assuming that the surface is an equivalent flat, impermeable boundary) and intrinsic nitridation rate constants decreases strongly with pressure.

2

ACCEPTED MANUSCRIPT

1. Introduction The entry into a planetary atmosphere is one of the most critical phases of

RI PT

space missions. During this phase the Thermal Protection System (TPS) is a key

component to guarantee payload integrity against the severe heat load due to the

harsh post-shock environment. For spacecraft entering at a very high velocity (higher than 10 km·s−1 ), engineers have to employ ablative materials to build effi-

SC

cient heat shields. Design of their thicknesses is a critical issue: the best trade-off

between thermal protection efficiency and heat-shield mass minimization needs to

M AN U

be found.

Carbon-based ablative materials are extensively studied, as they represent the most performing subclass of ablative TPS materials; here, we will consider lightweight ablators, made of carbon fiber preforms consolidated with a few percent of pyrolyzed phenolic resin, considered for moderately severe re-entry scenarios, will be addressed. An important characteristic of lightweight ablators is the possibil-

TE D

ity to ablate in-depth, due to their high porosity. Therefore, the gaseous species coming from the outer flow field do not just diffuse inside the boundary layer before reaching the surface but also through the porous material. Hence, in the

EP

experiment, the heterogeneous gas-solid reactions are not taking place entirely at the surface but rather on a finite thickness of material. The depth of this

AC C

region is directly related to the competition between reactant diffusion through the fibers and surface reactivity of the fibers in the porous material. Additionally, another gas-surfacesolid interaction phenomenaon aremay be enhanced by the large surface/volume ratio: catalytic recombination of dissociated gasatomic gaseous species, particularly important in TPS design because they areit is a strongly exothermic process. This phenomenon cannot be neglected in real reentry conditions because the extent of dissociation of the post-shock fluid is important. Its consequences are first a marked rise of the surface temperature with 3

ACCEPTED MANUSCRIPT

respect to non-catalytic conditions, and consequently an increase of the gas diffusion rate and a much larger increase of the heterogeneous reaction rates, therefore

RI PT

of the recession velocity. Most state-of-the-art CFD codes do not take into account the in-depth diffusion

of reactants and hence are necessarily limited to a macroscopic description of the

reactivity at the surface reactivity [1]. However, in a recent study [2], surface and

SC

volume ablation have both been accounted for in a single model, in which in-depth ablation is treated only macroscopically, in a volume-averaged wayframework. In contrast, numerical methods likesuch as Monte-Carlo Random Walks (MCRW),

M AN U

capable of taking into account the reactant diffusion through the fibers and the reaction with the fibers, allow for a microscopic approach of the in-depth diffusion/reaction process [3–6].

Summarizing these points, it is evident that further efforts to finely tune the design of TPS containing porous ablators are needed. The AblaCat project, supported

TE D

by ESA,The work reported here aims to fulfill these tasks by addressing two objectives: (i) to investigate possible catalytic properties of porous ablators, and (ii) to improve the computational models of ablation and catalysis for TPS design. A past study hads led to the determination of fiber-related intrinsic reaction

EP

constants, in the case of oxidation of fiber bundles [7]; in thisat work, a simple model for diffusion-reaction competition throughout the fibersfibrous medium

AC C

hads been used to recover the fiber-related reaction rate constant from the effec-

tive rate constant obtained by comparing the overall oxidation flux to the oxygen partial pressure close to the sample surface. We will extend hHere, the ideas developed in this previous work will be extended with a more accurate model that accounts for the presence of an ablation affected layer and works directlysimultaneously with the couple of two competing heterogeneous chemical reactions that are nitridation and catalytic recombination.

4

ACCEPTED MANUSCRIPT

The emphasis is set here on pure nitrogen plasmas. This allows a determination of rate constants linked to this speciesrelated only to this element (i.e. nitri-

RI PT

dation and nitrogen recombination). When the atmosphere contains oxygen, the reactivity being much higher, it is expected that volumetric ablation is not related

tosignificant for oxygen; on the other hand, nitridation and nitrogen recombination still take place and can develop their extent in the volume as in the absence

SC

of oxygen. So, knowledge of the nitrogen-related reaction rate constants is rele-

vant even when studying air plasmas. Moreover, some atmospheres (e.g., Titan)

relevance to the present study.

M AN U

for which entries have to be studied are very poor in oxygen, which further adds

In this article, we will report an experimental campaign of ablation tests performed in the von Karman Institute (VKI) Plasmatron facility [8–10] on a modelrelevant TPS material and a posterior examination of its surface state by X-ray Computerized Micro-Tomography (CMT) and image processing will be reported. AThis

TE D

allowed for a direct method for the identification of the affected ablation depth has been designed and will be depicted here. Then, we will describe ourMoreover, a numerical identification method allowing the retrieval of the ablation and catalytic rate constants from the experimental data has been designed and will be de-

EP

scribed. For thisTo this goal, computations at different levels have been chained: CFD codes allowed rebuilding the flowfree-stream conditions, and were coupled

AC C

to material recession and catalytic recombination simulations, able to reproduce the material surface state while matching the correct heat and mass fluxes. To help identifying the chemical kineticreaction rate parameters, we will also introduce a simple analytical model relating the depth of the affected zone to the intrinsic reactivity of the constituents and to the test conditions has been developed; application of the method to foresee the values of the effective reaction rates will be described and used as a first guess, inserted in the numerical code for final

5

ACCEPTED MANUSCRIPT

identificationnumerical experiments were then performed, resulting in a slight update to the analytical results; finally, the fiber-related reaction rate constants have

RI PT

been identified. The results of this novel approach will then be discussed and some of its perspectives of the approach will be given in the conclusion. 2. Experimental techniquesExperiments

SC

In this section we describe successively the tested material, the Plasmatron

facility and its in-situ diagnostics, and finally the post-test material analysis. An

M AN U

experimental campaign was carried out with the final aim of extracting information to be used in the theoretical/computational analysis described later on. This campaign was composed by: i) a pre-test material characterization; ii) a plasma testing campaign; and iii) a post-test characterization. The pre-test material characterization provided an average value for the pristine fiber diameter before plasma exposure. Plasma testing exposed the carbon preform to specific

TE D

(target) conditions and allowed to characterize its response in terms of surface temperature and recession rate. The post-test characterization was used to determine the depth of the zone affected by the gas-surface interaction.

EP

The different components of the experimental campaign are described in the following after a short description of the tested material. A short summary of the

AC C

experimental results concludes the present section. 2.1. Tested material and samples The tested material is a commercial carbon fiber preform, Calcarb r CBCF

18-20001 : it contains a small amount of pyrolyzed phenolic resin, used to bind

fibers. TheIts porosity is ≈ 90%, and the measured average density (average on all 1

https://www.mersen.com/sites/default/files/publications-media/3-gs-calcarb-grade-cbcf-18-

2000-mersen.pdf

6

ACCEPTED MANUSCRIPT

samples) was 176 kg·m−3 and ≈ 5-6 µm, respectively..measured average density (average over all tested samples) and fiber diameter read 176 kg·m−3 and ≈ 5-610

RI PT

µm, respectively. The average fiber diameter is ≈ 5-6 µm. Figures 1(a)-1(b) areis a renderings of the fibrous architecture as obtained by X-ray microtomography.

Here, it is clear ; they show that clusters of fibers alternate with large void regions and regions in which the fibers and microporosity are more evenly distributed.

SC

The tested samples (Figure 1(c)) were machined in house at the von Karman Institute (VKI) and consisted of a hemispherical nose (2.5 cm radius) and a cylindrical aft body (2.5 cm height). Figure 4 shows that the surface fibers are thinned

M AN U

away after a test. Even though mechanical erosion cannot be ruled out as a competing mechanism for surface recession, the fact that pointed fibers are present

AC C

EP

TE D

confirms the presence of chemical ablation.

7

9 2.

SC

2.3 mm

RI PT

ACCEPTED MANUSCRIPT

m

M AN U

m

m

4m

(a) 3D rendering of subsample.

(b) Zoom on an extract located in

AC C

EP

TE D

the center of the subsample.

(c) Full sample.

Figure 1: (a) 3D rendering of a pristine carbon-preform subsample (1700 × 800 × 800 voxels) obtained by X-ray CMT with pixel size of 2 µm; (b) Zoom on an extract located in the center of

the sample; (c) Axonometric-like picture of a pristine plasma-test carbon-preform sample.

8

ACCEPTED MANUSCRIPT

2.2. Plasma wind tunnel facilitytesting The VKI Plasmatron is a high enthalpy wind tunnel in which plasma is gen-

RI PT

erated by electromagnetic induction. This type of facilities, also known as Inductively Coupled Plasma (ICP) torches, thanks to their electrodeless technol-

ogy for plasma generation, present the advantage of producing high-purity plasma

flows. The Plasmatron uses a high frequency, high power, high voltage (400 kHz,

SC

1.2 MW, 2 kV) solid state (MOS technology) generator, feeding the single-turn

inductor of the 160 mm diameter plasma torch. The torch is mounted on a 1.4 m

M AN U

diameter, 2.5 m long, water-cooled test chamber, fitted with different portholes that allow unrestricted optical access to the test section. Hot gas from the test chamber exits through a 700 kW heat exchanger to a group of three rotary-vane vacuum pumps and a roots pump, which are capable of extracting 3900 m3 ·h−1 , with a terminal vacuum capability of 0.04 mbar. Thanks to this pumping system, the plasma jet flows into the test chamber at sub-atmospheric pressure (ranging

TE D

from 10 mbar to 800 mbar). A 1050 kW cooling system using a closed loop deionized water circuit (2090 l·min−1 ) and fan-driven air coolers provide cooling to all facility components. The facility is computer-controlled using a 719 I/O lines

EP

PLC, and two PCs for controlling and monitoring the Plasmatron operation. 2.2.1. Diagnostics and measurements

AC C

The fully-instrumented Plasmatron test chamber allows monitoring of the plasma free-stream conditions before injection of the test article into the jet. A copper calorimeter is used for the cold-wall heat flux measurement and a Pitot probe, connected to a static pressure port, measures the dynamic pressure of the flow. Additionally, when ablation tests are performed, optical measurements of quantities of interest such as sample recession rate, surface temperature, and spectrally resolved boundary layer emission are usually carried out through the optical accesses of the chamber [11]. In the present test campaign the non-intrusive setup 9

ACCEPTED MANUSCRIPT

shown in Figure 2 was used. It included a high-speed camera (HSC) for timeresolved recession measurements, an ICCD camera equipped with a spectral filter

RI PT

for free-stream temperature determination, a two-color pyrometer and a broadband radiometer for surface temperature and emissivity measurements, respectively. Among these, surface recession rate and temperature measurements are the most relevant for the present analysis. The short exposure time of the HSC

SC

(Vision Research Phantom 7.1 CMOS) allows the monitoring of the sample recession (resolution of ∼ 0.3 mm) despite the high brightness of the surface. A low

sampling frequency was applied (100 Hz), which allowed to record the whole test

M AN U

duration without filling up the internal memory of the HSC. The two-color pyrometer (MR1SC Raytek Marathon Series) operates over two overlapping wavebands between 0.95 µm and 1.1 µm, providing a measurement of the surface temperature which is independent of surface emissivity. The instrument lower and upper

TE D

measurement bounds are 1273 K and 3273 K, respectively. 2.2.2. Test conditions

The steady-state stagnation-point surface temperature of the test sample defined the target condition for the tests. Three nominal values of surface tempera-

EP

tures, 1500 K (low Plasmatron power), 2000 K (mid Plasmatron power), and 2800 K (high Plasmatron power) were selected. Tests have been performed in nitrogen

AC C

atmosphere (18 g·s−1 for the low-power tests and 16 g·s−1 for the mid- and highpower tests) for two neatly distinct static chamber pressures (1500 Pa and 25000

Pa). The actual value of the surface temperature was measured during the test with the two-color pyrometer (see section 2.2.1). The surface temperature is set by the complex balance of the energy fluxes on the sample surface, and the Plasmatron power needed to reach the target temperature was, in principle, unknown a priori.

The different Plasmatron electric powers necessary to reach each target temperature were, in principle, unknown a priori, being the surface temperature defined 10

function generator

SC

2-D ICCD camera

RI PT

ACCEPTED MANUSCRIPT

2-color pyrometer

M AN U

spectral filter

test sample exhaust & heat exchanger

TE D

plasma torch

broadband radiometer

AC C

EP

high-speed camera

Figure 2: Plasmatron set-up for in-situ measurement for hemispherical samples.

11

ACCEPTED MANUSCRIPT

by a complex balance of several energy fluxes. Therefore, a specific test procedure was defined based on the following steps: i) adjust the Plasmatron power

RI PT

to target a first tentative value of the measured cold-wall heat flux based on previous test campaigns on similar materials, ii) inject the sample into the plasma flow and monitor the surface temperature, iii) adjust the power (if necessary) to reach the target surface temperature (±100 K), iv) eject the samples after a pre-

SC

defined exposure time, v) measure the cold-wall heat flux. It is worth noting that

the test campaign included also test in similar conditions in air atmosphere. Being ablation more severe in air, it was decided to perform those tests first, setting the

M AN U

acceptable target temperatures within ±100 K from the nominal target ones (i.e., 1500 K, 2000 K, and 2800 K). This allowed to reduce the exposure time during the trimming. After all the air tests were performed, updated target temperatures were defined for the nitrogen tests (i.e., air test temperature ±50 K). This procedure proved to be reliable andas only minor adjustments to the plasma power

TE D

(step iii) were necessary to reach the targetsatisfactory steady-state surface temperatures. Time profiles of measured surface temperatures are shown in Figure 3 for the 6 nitrogen tests analyzed in this work, as well as for the 6 corresponding air tests. The average value of the achieved steady-state surface temperatures (uncer-

EP

tainty ±15 K coming from a certified calibration of the two-color pyrometer) are

AC C

reported in Table 21 along with other quantities measured during the experiments.

12

3000

2500

2500

1500

1000

2000

1500

1000 0

50

100

150

0

Time, s

SC

2000

RI PT

3000

Surface temperature, K

Surface temperature, K

ACCEPTED MANUSCRIPT

50

100

150

Time, s

(b) 250 hPa.

M AN U

(a) 15 hPa.

Figure 3: Time profiles of surface temperature for the tests at chamber static pressures of 1500 Pa (a) and 25000 Pa (b). a) 25000 Pa ; b) 1500 Pa. The decrease of the Plasmatron power, following the explained procedure, to approach the target steady-state surface temperature is evident in the

TE D

low- and mid-power tests at 25000 Pa.

Table 1: Plasmatron test conditions: test ID, static pressure p s (static port in the Plasmatron chamber), dynamic pressure pd (Pitot probe), meanaverage surface temperature T s (two-color pyrometer), cold wall heat flux q˙ cw (copper calorimeter), sample exposure time τ. p s (kPa)

pd (Pa)

T s (K)

q˙ cw (MW·m−2 )

τ (s)

P-N-15-T1500

1.65

58

1490

0.30

541

P-N-15-T2000

1.5

147

2087

1.25

244

P-N-15-T2800

1.56

338

2711

3.97

123

P-N-250-T1500

25

3

1636

0.42

720

P-N-250-T2000

25

8

2139

1.37

247

P-N-250-T2800

25

32

2756

3.48

182

AC C

EP

Test ID

2.3. Post-test samples characterization Tested-sample microscopic inspection (Scanning Electron Microscopy) showed the presence of pointed fibers near the surface, confirming the occurrence of chem13

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

Figure 4: Scanning Electron Micrograph of ablated-sample surface. Fibers are thinned away after

TE D

testing due to chemical ablation.

ical ablation in a kinetic-limited regime (see Figure 4).

EP

After weight and height variation determinations on the samples, subsamples were cut in their apical part, around the stagnation point. They were then stud-

AC C

ied by X-ray CMT using a GE - VTomeX M laboratory tomograph, with a 180 keV X-ray generator operated at 60 keV and 240 µA with a 0.1 mm copper filter. 1500 transmission images were acquired rotating the sample onby 360°, averag-

ing 2 images and with an exposure time of 2000 ms. The resolution was 2 µm per voxel, enough to properly distinguish the fibers in the structure, as illustrated in Figure 1(a). The purpose of this post-test characterization was the determination of the thickness of the region chemically affected by the exposure to the plasma flow. However, finding the depth of the affected zone is not trivial, since 14

ACCEPTED MANUSCRIPT

the material porosity is highly variable throughout the samples, with characteristic lengths of variation often larger than the affected zone size. Another difficulty

RI PT

is that the solid volume fraction vanishes close to the surface, making it difficult to identify the precise location of the extreme surface.to precisely identify the material edge location. To bypass these difficulties, a specific image processing routine has been developed. It features the following steps: First, tThe image is

SC

firstthresholded based on the tomography grayscale distribution. The present fi-

brous material, composed by voids and fibers only, shows a bi-modal distribution. Setting the threshold value halfway between the two modes provides an accurate

M AN U

segmentation of the two phases. A porosity (or solid volume fraction) profile may be drawn from the resulting image. In every slice taken parallel to the surface, the intersections of the fibers with the slice planes should be ellipses. In practice, ellipses are fitted to the 2D objects or “particles,” of which weone can measure the area and smallest diameter. The built-in routines of the Fiji distribution [12] of

TE D

ImageJ [13] are used for this purpose. They obtain a maximal similarity of the first and second moments of the pixels distributions with respect to the ellipse center position and axes, additionally requiring an exactly equal area. However, it can happen that tThe measured values are not all acceptable, i.e., those which are too

EP

small are discarded (high-pass filtering: areas below one half of the expected fiber section are excluded), as well as those which correspond to fiber clusters (low-

AC C

pass filtering: diameters above twice the expected fiber diameter are excluded) are discarded. Then, aAn average fiber diameter is computed on the filtered particle population for each slice. The average diameter and is plotted against the

depth in the material. Finally, tThe affected depth is obtained visually from the plots of average diameter and porosity vs. depth. Figure 5 is an example of plots giving the average particle diameter and porosity vs. depth in the sample. The affected depth appears rather clearly: it is the

15

ACCEPTED MANUSCRIPT

depth after which the fiber diameter and solid volume fraction reach the charac-

AC C

EP

TE D

M AN U

SC

RI PT

teristic value of pristine samples and do not increase any more.

16

ACCEPTED MANUSCRIPT

SC

RI PT

Preform, N2, 1500 Pa, 1500 K

a)

M AN U

12

Affected depth δa

Fiber diameter (µm)

11 10 9 8

6 0 0.16 0.14

TE D

7

100

200

300 400 500 600 700 Distance from surface (µm)

b) 800

900

1000

0.1

0.08

AC C

Solid volume fraction (-)

EP

0.12

0.06 0.04

c)

0.02 0

0

100

200

300 400 500 600 700 Distance from surface (µm)

800

900

1000

Figure 5: Analysis of the preform sample P-N-15-T1500 scanned after a plasma exposure torch test under 15 hPa N2 at 1500◦ C close to the surface. a) Projected view; b) Plot of the average fiber diameter (see text for details) vs. depth in the porous material. c) Plot of the solid volume fraction profile vs. depth in the porous material.

17

ACCEPTED MANUSCRIPT

2.4. Experimental results Table 2 summarizes all tests in which the same above-described material was

RI PT

exposed to nitrogen plasmasthe experimental results. Values between parentheses have limited confidence because they are close to the detection limit. We can

seeIt is interesting to note that the recession rate increases and the affected depth decreases with temperature for both pressures. In addition, tThe recession rate is

SC

significantly higher at lower pressure.

Table 2: Plasmatron test conditions and results: test ID, static pressure p s , dynamic pressure pd ,

M AN U

meanaverage surface temperature T s , cold wall heat flux q˙ cw , sample exposure time τ, recession rate vr and measured affected zone depth δa .

Test ID

p s (kPa)

T s (K) vr (µm·s−1 ) δa (µm)

1.65

1490

2

464

P-N-15-T2000

1.5

2087

5

232

P-N-15-T2800

1.56

2711

17

266

P-N-250-T1500

25

1636

(0.6)

460

P-N-250-T2000

25

2139

(1.1)

326

P-N-250-T2800

25

2756

2

272

EP

TE D

P-N-15-T1500

AC C

3. Identification of the reaction rate constants The whole identification strategy involves the reconstruction of the flow con-

ditions outside the sample up to the sample surface, and the simulation of the chemical phenomena inside the porous sample. In order to identify the reaction rate constants from the material behavior and knownthe measured/reconstructed conditions at the material surface, we have used first a simple analytical model has

first been introduced for the production of crude estimates; then, the first-guess values are inserted in a detailed numerical model of the sample recession based 18

ACCEPTED MANUSCRIPT

on Monte-Carlo random walks has been applied for a more accurate identification. The obtained closed-form relationships allowed for a direct identification of the

RI PT

sought reaction rate constants. 3.1. A simple analytical model

The considered model is unidimensional, considering onlyonly accounting for the depth inside the sample, z. We will consider that the solid gasification reaction

SC

has a small heat of reaction, as compared to the recombination reaction.Heat trans-

fer effects will not be accounted for, considering a constant temperature over the

M AN U

relatively small region of interest (hundreds of µm at most). Let us consider the porous medium as homogenized with a porosity distribution along the depth below its surface, ε(z). Its effective internal surface area σv given in m2 ·m−3 , or m−1 , is the ratio between the specific surface area and the mass density and is a key parameter of the problem because it strongly influences strongly the importance

TE D

of heterogeneous reactions, as compared to transfer phenomena. At any position, σv is known, since it is related to the porosity by the following relationship: σv =

2p (1 − ε) (1 − ε0 ) r0

(1)

EP

This relationship has been obtained considering separate cylindrical fibers with a constant spacing (i.e. in a square array) and variable radius, starting with an initial

AC C

radius r0 and porosity ε0 . See Appendix A for details. On the other hand, it can be considered that the effective diffusion coefficient D p of the active species D p is

linked to the porosity and to its pure gas diffusion coefficient as: ε D p = Dg η

(2)

where η is the tortuosity of the porous medium, a quantity that does not vary much in the interval [ε0 ; 1]. For sake of simplicity, we will assume η = 1 will be assumed in the following. 19

ACCEPTED MANUSCRIPT

The gas mass balance inside the porous medium results from the competition between diffusion and the heterogeneous chemical reactions with rate constants

RI PT

knit and krec for nitridation and recombination, respectively:   ∂εC + ∇ · −D p ∇C = −σv (knit + krec ) C ∂t where C is the active gaseous species concentration.

(3)

C(z = 0) = Cb

SC

The boundary condition at the porous medium surface is:

(4)

M AN U

The solid mass balance yields an evolution equation for the porosity: ∂ε = Ω s σv knitC ∂t

(5)

where Ω s = M s /ρ s is the molar volume of the pure solid phase, i.e. excluding the contribution of porosity. In this equation, we see that the algebraic porosity

solid molar volume.

TE D

increase is directly linked to the ablation nitridation reaction molar flux, times the

The material recedes with a given velocity, vr , and we seek solutions are sought under the form of a traveling wave, that is, a steady profile of concentration and

EP

porosity traveling at constant speed vr . A reference frame linked to this traveling wave (z0 = z − vr t, t0 = t) will transform the time operator ∂/∂t into a space

AC C

operator multiplied by the velocity −vr d/dz0 . We therefore have tThe following couple of ODEs is obtained instead of eqs. (3) and (5): !   d dC      dz0 −D p dz0 − vr εC = −σv (knit + krec ) C     dε     −vr 0 = Ω s σv knitC dz

(6a) (6b)

Combining both equations one gets: ! ! d dC vr dε krec − 0 D p 0 + vr εC = 1+ dz dz Ω s dz0 knit 20

(7)

ACCEPTED MANUSCRIPT

Let us introduce κ = 1 +

krec . knit

A prime integral appears:

  dC vr ε κΩ−1 + C + D p 0 = I1 s dz

RI PT

(8)

The value of the prime integral is obtained by looking at the z0 → ∞ limit, for which the gas is completely depleted, so: I1 = vr ε0 κΩ−1 s

SC

(9)

Therefore eqs. (2), (8) and (9) yield a relationship between the porosity and

ε

M AN U

concentration profiles:

Dg dC −1 + εκΩ−1 s − ε0 κΩ s + εC = 0 0 vr dz

(10)

We will now introduce tThe variable φ = 1 − ε, i.e. the solid volume fraction, will now be introduced, since it which will be easier to handle than the porosity, and the following dimensionless variables:

TE D

 C   χ=    Cb     z0    ξ= δa

(11a) (11b)

EP

Equation (10) can be rewritten as:

AC C

! ρc dχ (1 − φ) ∆ + χ − φ + φ0 = 0 κ dξ

(12)

where two dimensionless quantities appear:  ρc = Ω s Cb     Dg     ∆= δa vr

(13a) (13b)

The first one, ρc , is a condensation ratio (i.e. the ratio between gas and solid concentrations) and is expected to be small with respect to unity; the second one, 21

ACCEPTED MANUSCRIPT

∆, compares the diffusion rate through the affected zone depth to the recession rate.

(eq. (1)), one gets finally:

(14a)

(14b)

SC

 dφ    =Aχ √     2 φdξ    dχ κ φ0 − φ    +χ=−  ∆ dξ ρc 1 − φ

RI PT

Combining with the solid mass balance (eq. (6b)) and the surface area law

M AN U

where we have defined some other dimensionless groups have been defined:

p  2  φ0 ρc A = ∆Sh λ  nit      k r    Shnit = nit 0  Dg      δa     λ= r0

(15a) (15b) (15c)

From eq. (14a) one gets a simple relation between the concentration profile

TE D

χ(ξ) and the solid volume fraction profile φ(ξ):

 ξ 2  Z    φ(ξ) = A χ(ξ)dξ  

(16)

0

EP

We therefore can approximate tThe solution can be approximated by choosing only a candidate function for χ(ξ). A simple choice is, in the spirit of the boundary

AC C

layer theory, a 3rd -order polynomial over [0; 1]: χ(ξ) = 1 −

3 1 ξ + ξ3 2 2

(17)

By inserting the condition that φ(1) = φ0 , one gets: A=

8p φ0 3

(18)

or : Shnit λ2 ∆ρc = 22

8 3

(19)

ACCEPTED MANUSCRIPT

Substituting the profiles of χ and φ in eq. (14b) and evaluating at ξ = 0, a new relation appears:

κφ0 3 ∆−1= 2 ρc

RI PT

(20)

Noticing that ∆ >> 1, eq. (20) rewrites: κφ0 3 = ∆ρc 2

(21)

SC

The full problem solution is contained in eqs. (19) and (21). The actual values

of the numerical coefficients in their right-hand side terms are linked to our choice

M AN U

of eq. (17) for the concentration profile χ(ξ); they are expected to differ in actual cases, as we will seewill be seen later.

Combining eqs. (19) and (21) gives another interesting relationship: κShnit λ2 φ0 = Shλ2 φ0 = 4Θ2 = 4

(22)

where we have definedappears the Sherwood number associated to the sum of the

also defined:

TE D

reaction ratesconstants Sh = κShnit . We have also defined a A Thiele modulus is Θ=

δa 2

s

(knit + krec ) φ0 λ p φ0 Sh = r0 Dg 2

(23)

EP

which has to be unity in the considered case. The Thiele modulus appears naturally in diffusion/reaction problems [14]; when it is unity, diffusion through the

AC C

chosen depth (here, δa /2) has the same rate as the heterogeneous reaction (here, nitridation and recombination rates are lumped together). It is now possible to retrieve the values of the nitridation and recombination

rate constants from observables and the computed values of the boundary concentrations. For this, it is useful to define the effective nitridation rate : eff knit =

φ0 vr Ω sC b

23

(24)

ACCEPTED MANUSCRIPT

This is the apparenteffective surface reaction rate of the material considered as flat and homogeneous. Equation (19) can be rewritten as : (25)

RI PT

knit

eff 8 vr 1 8 knit r0 = = 3 Ω sCb λ 3 φ0 δa

Very interestingly, we can notice that in our caseit appears that the intrinsic nitri-

dation rate of the fibers, knit , is not very different from the effective nitridation rate

SC

eff eff of the material, knit : indeed, r0 /δa ≈ 10−2 and φ0 ≈ 10−1 , so knit is ≈ 10 · knit .

The intrinsic, fiber-scale recombination rate constant, krec , is retrieved through

M AN U

the evaluation of κ, thanks to eq. (21): κ=

3 Dg 2 δa keff

(26)

nit

Physically, we see an interesting point is that diffusion through the affected depth, in competition with the effective nitridation, defines the recombination/nitridation ratio. Indeed, when recombination is important, the penetration through the porous

TE D

medium, for the same effective recession rate, is lesser. The value of krec is, as a function of observables :

krec

4r0 Dg 2 eff = − k φ0 δa δa 3 nit

! (27)

EP

Reciprocally, knowing the elementaryintrinsic, fiber-scale rate constants, one can retrieve the ablated thicknessaffected depth and the effective nitridation rate

AC C

constant:

 !−1/2   2r0 Ω sC b    δa = p 1+    κφ0 φ0 Sh   r !−1/2    φ 3 Ω C  0 s b eff     knit = 4 knit Sh 1 + κφ

(28a) (28b)

0

The main conclusions of this analytical model are the following: This analytical model calls out some comments. First, it is clear from eq. (28a) that the affected depth is depending on the total reaction/diffusion ratio Sh. 24

ACCEPTED MANUSCRIPT

This ratio is itself dependent on the effective radius of the fibers r0 : this is the most important contribution of the porous medium structure and geometry, in ad-

RI PT

dition to its initial porosity ε0 = 1 − φ0 . Also note that this ratio concerns both reactions, nitridation and recombination, and not only the one that directly alters the structure of the porous medium. Second, the ratio between the effective and

the fiber-related reaction constants is directly linked to the ratio λ between the

SC

affected depth and the fiber radius (see eq. (25)), and to the initial solid fraction.

trogen recombination

M AN U

3.2. Detailed numerical simulations of material recession and radicalatomic ni-

The simulations are carried out using three codes matching three levels of description. First, CFD computations including fully coupled Maxwell and Navier-Stokes equationsNavier-Stokes equations with electromagnetic field equations allow modeling the plasma hydrodynamics produced by the wind tunnel andwithin the plasma

TE D

torch and the test chamber to rebuild the free-stream conditions upstream of the sample; second, using the previously reconstructed free-stream conditions, chemically reacting Navier-Stokes CFD computations coupled to a gas-surface interaction modeling based on macroscopic surface mass and energy balances provide

EP

identification of the effective reaction constants and the active gaseous species concentration at the wall; finally, the Monte-Carlo Random Walks/Simplified March-

AC C

ing Cube (MCRW / SMC) code dedicated to the details ofmodeling inside the porous medium is used to retrieve the local ratereaction constants intrinsic to the

fibers.

3.2.1. CFD computations The Plasmatron free-stream test conditions upstream of the sample have been rebuilt using the VKI procedure. It is composed by the two following steps. First, the subsonic Plasmatron flow field is numerically simulated using a resistive mag25

ACCEPTED MANUSCRIPT

netohydrodynamics solver, characterizing the boundary layer geometryweaklyionized plasma flow around the heat flux probe (considered isothermal at 350 K

RI PT

in this step) under local thermodynamic equilibrium and axisymmetric flow assumption (VKI ICP code [15–17]). Input to this first step are the gas mass flow rate, electrical power, and static pressure ps in the test chamber. From this simula-

tion, hydrodynamic parameters characteristic of the flow are extracted and, along

SC

with the experimentally determined heat flux and Pitotdynamic pressure, serve as

input conditions for the second step of the rebuilding procedure. This step uses the VKI boundary layer code [18] to solve the axisymmetric, chemically-reacting,

M AN U

stagnation-line boundary layer over a catalytic surface (i.e., the copper calorimeter used to measure the heat flux during the test) under chemical nonequilibrium conditions. A catalytic boundary condition is applied to the copper surface (recombination efficiency 0.019) and a Newton method is used to iterate on the boundary layer outer edge temperature T e , until the numerical heat flux matches

TE D

the experimental measurement. Interested readers are referred to [19] for a more thorough explanation. An uncertainty quantification analysis to study the impact of the experimental uncertainties on the rebuilt flow and their propagation in a CFD ablation model was recently presented in [20].

EP

Using the evaluated boundary-layer edgefree-stream conditions, a second threedimensional chemically reacting Navier-Stokes CFD code [21], called AFFS (Ab-

AC C

lation Fortran Flow Solver), works on thecomputes the flow-field around the sam-

ple up to the surface, solving species mass and energy balances at the interface between the material and the outer fluid. It considers through a dedicated ablative boundary condition that the surface is an equivalent flat, impermeable boundary, out of which gases are blown. The species mass balance equations at the surface

26

ACCEPTED MANUSCRIPT

RI PT

are (considering that the normal to the surface points from the gas to the solid): ! N sp X ∂yi α j,rec α j,nit eff eff ρDi,m + ρvb yi = ρ krec (T )ν j,rec y j + knit (T )ν j,nit y j (29) ∂z j=1 ˙ of the Here, the blowing ratevelocity vb is directly related to the mass loss rate m sample: vb = m/ρ ˙

(30)

SC

3.2.2. MCRW/SMC computations: updates to the analytical model results

The detailed computations of the ablation of the material have been performed

M AN U

with the AMA software [3], which has been developed on the principle of MonteCarlo Random Walks (MCRW) to solve a diffusion-heterogeneous reaction problem, in which the interface is explicitly described and can move according to ablation events. AMA, which has been implemented in ANSI C, contains five main parts:. First, a 3-D image containing several phases (fluid/solids) is described by

TE D

discrete cubic voxels (3D pixel). The moving fluid/solid interface is determined by a Simplified Marching Cube (SMC) approach [22]. Mass transfer by diffusion is then simulated by a Brownian motion simulation technique, which is a continuum (grid-free) method to simulate diffusion in which the reactant species is

EP

discretized into random walkers. Heterogeneous first-order reaction on the wall is simulated by a sticking probability adapted to the Brownian motion simulation

AC C

technique. Surface recession is treated by decreasing the grayscale value of the grid node lying closest to the location of the sticking event; when this grayscale value reaches zero, the SMC triangulation is applied locally around the modified node.

A Dirichlet upper boundary condition is simulated using a buffer zone,

where C is maintained constant at C0 . The boundary value Cb is lower than C0 because of gas diffusion between the buffer and the surface and is retrieved by post-processing the results. AMA has been validated by comparison to a 1D analytical model in transient regime and to the 3D mesoscopic scale model in steady 27

ACCEPTED MANUSCRIPT

state. Independent implementations of the same model have been successfully compared to AMA [3, 5, 23].

RI PT

For the purposes of the present work, AMA has been adapted to handle simultaneously nitridation and recombination events. i.e. random walker sticking

events that lead or do not lead to surface recession. The original code was handling only the equivalent of nitridation, i.e. an image grayscale level decrease was

SC

systematically associated to every sticking event. In the improved version, with a probability ratio of Prec / (Pnit + Prec ), the image remains unaltered at the stick-

ing event. The sticking probability for any walker is now the sum of the sticking

M AN U

probabilities for nitridation Pnit and recombination Prec .

The code accepts as input values the concentration of active species near the surface, the diffusion coefficients and kineticreaction rate constants (from which the rate efficiency factors are computed), and the geometry of the image. As an output, in addition to details on the ablated porous medium geometry, the

TE D

thickness of the ablated depth is obtained, using a frame-to-frame image difference detection routine, as well as the effective nitridation rate (or equivalently, the mass loss rate)recession velocity (linked to the mass loss rate through the porous medium density or to the effective nitridation rate through eq. 24). No blowing is

EP

accounted for in AMA; however, using eq. (36)a blowing correction derived from the above mentioned CFD ablation model can be introduced and will be presented

AC C

later on. It allows a good match with the near-surface mass balance modelusing the macroscopic description to correct the effective nitridation rate constant for the blowing.

4. Computationsal results and discussion The results of the CFD computations will be presented first; then, the outcome of the numerical MCRW/SMC method will be discussed. Chaining the two types 28

ACCEPTED MANUSCRIPT

of results requires some corrections, that will be discussed in a third section. Finally, the interpretation of the experiments is given under the form of intrinsic

RI PT

reaction efficiency factors. 4.1. CFD analysisresults

The output of the VKI rebuilding procedure, as computed with the VKI boundary layer software, is the size of the boundary layer and the physical values of free-

SC

stream conditions upstream of the sample in terms of bulk velocity, pressure, temperature, and species mass fractions at the boundary layer outer edge (see Table 3).

M AN U

These freestream conditions are used as inflow constant conditions in order to perform the subsequent subsonic CFD simulations with a reasonably large domain upstream of the sample.

Table 3: Rebuilt boundary-layer edgefree-stream conditions for the nitrogen plasma tests. T ∞ :

TE D

temperature; p∞ : pressure; h∞ : enthalpy; y∞,i : species mass fractions (i = e− , N+ , N, N2 ). T ∞ (K)

p∞ (Pa)

u∞ (m·s−1 )

h∞ (MJ·kg−1 )

y∞,e−

y∞,N+

y∞,N

y∞,N2

P-N-15-T1500

5745

1650

367.7

25.11

4.89·10−9

1.25·10−4

5.15·10−1

4.85·10−1

P-N-15-T2000

10603

1500

1245.4

90.26

1.32·10−5

3.35·10−1

6.65·10−1

2.48·10−5

P-N-15-T2800

12172

1560

2419.1

141.7

2.94·10−5

7.48·10−1

2.52·10−1

5.95·10−7

EP

Test ID

6175

25000

20.8

32.67

2.71·10−9

6.91·10−5

3.02·10−1

6.98·10−1

P-N-250-T2000

6675

25000

41.8

40.52

1.22·10−8

3.10·10−4

5.47·10−1

4.52·10−1

25000

110.5

52.44

4.87·10−7

1.24·10−2

9.72·10−1

1.53·10−2

AC C

P-N-250-T1500

P-N-250-T2800

8717

Figure 6 illustrates the results obtained with AFFS, the CFD code including

ablation, (second step of the procedure described above)using the free-stream conditions provided by the VKI flow rebuilding methodology (first step of the

procedure described above). As can be seen, the grid is sufficiently large to have

29

ACCEPTED MANUSCRIPT

a practically undisturbed flow condition at inflow. The mass fractions (and consequently the concentrations) and fluxes of the chemical species are determined at

RI PT

the reacting sample surface. As shown in Figures 6(c-d), the only gaseous species present in a non-negligible amount at the sample surface are CN, N and N2 . Fur-

thermore, material recession is due solely to nitridation, with carbon sublimation

AC C

EP

TE D

M AN U

SC

(indicated by C3 production at the surface) being entirely negligible.

30

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

Figure 6: Computation of test case P-N-15-T1500 using the second codeAFFS. a) Pressure and temperature; b) Close-up of a) near the sample; c) and d) mass fractions of several species : CN, C3 , N2 , N.

31

ACCEPTED MANUSCRIPT

4.2. MCRW/SMC results We have takenThe images used in these computations were either ideal images

RI PT

made of parallel cylinders aligned along the flow direction and perpendicularly to it, or X-ray µ-CT images extracted from the µ − CT scans presented in section

2.3. Three images with solid volume fractions of 9.1%, 8.1% and 6.6%, respec-

tively, were obtained. The images size was 300 × 300 × 300; a 10-voxel thick

SC

buffer zone has been added above the surface, spaced by a 10-voxel thick pure

gas phase zone. The computations in AMA have been repeated for several values

M AN U

of the kineticreaction rate constants knit and krec and of the diffusion coefficient, and the resulting effective reactivities and ablation thicknesses were compared to the CFD/mass balance solver results and to the experimental measurements, respectively.were extracted and cast in dimensionless form in order to compare them with the analytical model results.

Figure 7 illustrates a snapshot of the detailed results of AMA computations,

TE D

obtained after the establishment of a steady state. The frame-to-frame porosity difference curve helps determining the ablated thickness δa . The concentration profile is not very different from the 3rd -order polynomial taken in the simplified model (eq. (17)), thoughand some irregularities can be encountered. The peak

EP

in the concentration profile around 600 µm can be explained by the presence of a comparatively more porous region as can be seen on the grayscale profile. This

AC C

suggests to further verify the agreement between the numerical and the analytical models. Eqs. (28a-28b) were tested against the obtained numerical dataobtained

on three different images and with varying nitridation and recombination constants values. Figure 8 is a plot of the dimensionlessscaled affected depth λ against (φ0 Sh)−1/2 , as suggested by eq. (28a), taking into account the fact that the last term in parentheses on its right-hand side is very close to 1. The slope is approximately 5 instead of 2: this may arise from (i) the deviation of the concentration profile

32

RI PT

ACCEPTED MANUSCRIPT

Ablation rate (GS difference x frame rate )

SC

4

M AN U

2

0

800

Avg. grayscale

Concentration (walkers/slice)

400

TE D

0

0

400

1200 800 1600 Z Position in image (µm)

2000

Figure 7: Example of AMA output. Left: rendering of the fibrous preform sample ablated on its exposed surface (left) during test P-N-15-T1500; Right, top: porosity (meanaverage grayscale

EP

per slice) profile (gray curve) and numerical ablation rate (frame-to-frame porosity difference × frame rate) (blue curve) along the sample thickness. Right, bottom: snapshot of the numerical N

AC C

concentration profile (in walkers per slice) along the sample thickness. The shaded area represents the affected zone.

33

ACCEPTED MANUSCRIPT

180

Scaled affected depth l = da /r0 (-)

160

RI PT

140 120

y = 5.045 x R² = 0.990

100 80

y = 4.538 x R² = 0.9996

Image 2

SC

y = 4.695 x R² = 0.997

60

Image 1

Image 3

Ideal, parallel

40

Ideal, perpendicular

0 0

5

10

M AN U

20

15

20

25

30

35

40

(f0 Sh)-1/2 (-)

Figure 8: ScaledDimensionless affected depth computed with AMA vs. (φ0 Sh)−1/2 for different porous preform images and various Shnit and κ ratios in different parameter ranges for the different

TE D

images. Images 1, 2 and 3 have been extracted from the interior part of CT scans of samples tested at low pressure, avoiding the ablation affected region. Ideal images are made of a square array of parallel cylinders, oriented either parallel or perpendicular to the flow direction.

EP

from a 3rd -order polynomial, and (ii) deviations of the porous medium from the ideal model ( fiber clustering, dead-end pores, etc . . . ).Instead of the value of 2

AC C

expected from the analytical model, the slope is between 4.5 and 5, depending on the type of porous medium (ideal or real). The most plausible cause for such a discrepancy arises from the existence of a deviation of the concentration profile from a 3rd -order polynomial. Is has to be noted, on the other hand, that the influence of the geometric structure of the medium has little impact on the value of the slope. Another point of the theory that can be checked is the ratio between the local, fiber surface-related nitridation constant and the effective nitridation rate constant, 34

ACCEPTED MANUSCRIPT

given by eq. (28b), also taking into account the fact that the last term between brackets on its right-hand side is very close to 1. As shown in figure 9, instead

RI PT

of finding the expected square root dependency to φ0 /Sh, a 3/4-power dependency is found. Again, this fact probably arises from the ”non-ideality” of the porous medium. The difference of the exponent could again arise from the fact

that the concentration profile is not strictly a polynomial. For the ideal medium,

SC

the slope of the correlation varies from 0.56 in the parallel direction to 0.34 in the perpendicular direction. This is a marked contrast, denoting that the effective

reactivitynitridation rate constant has a stronger sensitivity to the structural details

M AN U

of the porous medium than the affected depth. The real porous medium images present a common trend, indicating the the three samples have similar structures, and the common slope is close to the one of the ideal medium with fibers perpendicular to the flow direction. This indicates that this type of porous medium has

AC C

EP

section 4.3).

TE D

more fibers perpendicular than parallel to the flow, as will be confirmed later (see

35

RI PT

ACCEPTED MANUSCRIPT

1.2

Scaled effective nitridation constant eff knit /knit (-)

y = 0.308 x R2 = 0.983

y = 0.563 x R² = 0.951

0.8

0.6

M AN U

1.0

SC

1.4

Image 1 Image 2

0.4

Image 3

0.2

0.0 0.0

Ideal, parallel

TE D

y = 0.337 x R² = 0.729

1.0

2.0

Ideal, perpendicular 3.0

4.0

5.0

(f0 /Sh)3/4

EP

eff Figure 9: Effective nitridation constant divided by the local fiber nitridation constant knit /knit com-

puted with AMA vs. (φ0 /Sh)3/4 for different porous preform images and various Shnit and κ ratios.

AC C

Images 1, 2 and 3 have been extracted from the interior part of CT scans of samples tested at low pressure, avoiding the ablation affected region. Ideal images are made of a square array of parallel cylinders, oriented either parallel or perpendicular to the flow direction.

36

ACCEPTED MANUSCRIPT

Accordingly, we can it can be concluded that the computations are in fair agreement with the predictions of the analytical model, though with slight devi-

RI PT

ations that can be attributed to the fact that some hypotheses taken for the model are not verified in the computationsthe porous medium is not ideal. A slightly differentnew equation system, that slightly differs from as compared to eqs. (28a-

SC

28b) is obtained:

M AN U

 5.045r0    δa = p     φ0 Sh    φ 3/4    0 eff   knit = 0.308knit Sh

(31a) (31b)

Remembering that Sh = κknit r0 /Dg , this system can be solved in terms of knit and κ, so that it is possible to identify the local nitridation and recombination kinetic rate constants, using as input the experimental recession rate and affected depth:

TE D

 !3/2   r 0 eff,corr    knit = 36.791knit    δa φ0  r    Dcorr krec φ0  g   = 0.692 eff,corr κ =1+    knit r0 δa k

(32a) (32b)

nit

subsection.

EP

where the superscript “corr” refers to data corrected as described in the following

AC C

Let us discuss an interesting example of utilization of eqs. (32a-32b). Exrayon fibers are known to be much more reactive than ex-PAN fibers. This has been proved experimentally with respect to oxidation [7]: the reactivity ratio is 20 between fibers. This comes from the nanoporous, disordered structure of the exrayon fiber, as compared to the denser, more ordered structure of the ex-PAN fiber. The same thing can certainly be expected with nitridation and recombination. As a consequence, for exactly the same other conditions (and assuming a constant recombination/nitridation ratio), replacing ex-PAN fibers by ex-rayon fibers would 37

ACCEPTED MANUSCRIPT

lead to a size reduction of the affected depth down to 22% ((i.e. 20−1/2 ) of the initial value. Then, the effective rate constant would be 201/4 times larger, i.e. 2.15

4.3. Corrections to the porous medium behavior model

RI PT

times the case of ex-PAN fibers.

The superscript “corr” has been introduced in eqs. (32a-32b) to indicate that

some corrections, necessary to account for some physical phenomena neglected

SC

in the presented porous medium model, were used in the derivation of those equations. A first correction has been introduced to account for the “blowing effect”

M AN U

linked to the mass imbalance of the products vs. reactants. Considering only N, N2 and the ablation product CN as the main gaseous species present at the sample

TE D

surface gives, from eq. (29):    ∂yN  eff eff   ρDN,m + ρvb yN = − krec + knit ρyN    ∂z      ∂yN2  eff + ρvb yN2 = krec ρyN ρDN2 ,m    ∂z      MCN eff ∂yCN    + ρvb yCN = k ρyN  ρDCN,m ∂z MN nit

(33a) (33b) (33c)

Summing up these three contributions, and recalling that the diffusive fluxes can-

EP

cel out by definition, we haveone has:

AC C

! MCN eff ρvb = − 1 knit (T )ρyN MN

Substitution in eq. (33a) gives: " ( ! )# MCN ∂yN eff eff = − krec (T ) + knit (T ) 1 + − 1 yN ρyN ρDN,m ∂z MN

(34)

(35)

We It can therefore be considered that N obeys a simple diffusion-reaction boundary condition, in which the corrected nitridation constant is itself a function of composition: ( eff,appcorr knit

=

eff knit (T )

! ) MCN 1+ − 1 yN MN

38

(36)

ACCEPTED MANUSCRIPT

Taking this dependence into account, the analytical model can be compared to what will be presented in the next section.exploited. However, a last correction

RI PT

has to be made. The identification of the rate constants has been carried out for the six tests. The value of the diffusion coefficient inside the porous medium is corrected by the Knudsen effect as

where the Knudsen number is Kn =

` , dp

(37)

SC

Dcorr = Dg (1 + Kn)−1 η(Kn)−1 g

` being the mean free path and d p the pore

M AN U

diameter, taken here as 4 (1 − φ0 ) /σv0 , and η(Kn) is the porous medium tortuosity, computed as a function of Kn with DMC, a validated home-made software [24, 25], as illustrated in figure 10. Note that the tortuosity values are in any case comprised between 1 and 2, which confirms the η = 1 choice made in the analytical model as a correct crude approximation. Moreover, it is seen that the direction z, parallel to the surface normal, has a larger tortuosity, which confirms

direction.

TE D

the fact that the fibrous medium has more fibers aligned perpendicular to the flow

4.4. Determination of the fiber-scale reaction constants

EP

It is customary to characterize surface reactions in terms of reaction efficiency

AC C

γi . In general, the effective rate constants are decomposed as follows: s  E  1 8RT 0 i kieff = γi exp − 4 πMi | RT {z }

(38)

γi

where the term in square root is the mean molecular thermal speed and γi is a wall

collisionreaction efficiency factor.

39

ACCEPTED MANUSCRIPT

RI PT

Sample surface

2.5

SC

2 1.5 1

M AN U

Diffusion tortuosity (-)

3

0.5 0 0.001

0.01

0.1 1 10 Pore Knudsen number (-)

100

1000

TE D

Figure 10: Tortuosity with respect to gas diffusion computed as a function of the Knudsen number.

EP

The direction transverse to the surface (squares) is the most tortuous one.

Table 4: Results for the preform exposed to nitrogen.

Kn

Dcorr g

yN

eff,corr knit

δa

γnit

γrec

eff γnit /γnit

(-)

(m2 ·s−1 )

(-)

(m·s−1 )

(µm)

(-)

(-)

(-)

P-N-15-T1500

0.972

0.0111

0.47

0.449

464

0.00125

0.150

1.46

P-N-15-T2000

2.04

0.0165

0.94

1.210

232

0.00462

0.375

3.07

P-N-15-T2800

1.77

0.0181

0.72

4.12

266

0.0235

0.606

4.67

P-N-250-T1500

0.064

0.00182

0.066

0.0290

430

0.000106

0.0313

1.51

P-N-250-T2000

0.084

0.00264

0.281

0.0259

326

0.000134

0.0729

2.90

P-N-250-T2800

0.108

0.00373

0.685

0.0306

272

0.000284

0.174

7.54

AC C

Test ID

40

ACCEPTED MANUSCRIPT

The results of the analysis performed on the six tested samples are summa-

RI PT

rized in Table 4. It is seen that the wall collision reaction efficiency factors γi h i are in the 10−4 ; 0.6 range. The nitridation constantsefficiency factors increase with temperature and they are markedly lower at high pressure. Figure 11 repre-

sents the nitridation probabilitiesreaction efficiencies, compared to literature data [26–31]. It can fairly be seen that the values found at low pressures are very

SC

consistent with other previous ones and have a similar temperature evolution (in

particular in with those from R refs. [28] and [31]). On the other hand, the values obtained at high pressures are one order of magnitude below literature data

M AN U

values, though displaying approximately the same temperature evolution; moreover, the recombination rate constantsefficiencies are found higher than the nitridation rate constantsefficiencies by 2 orders of magnitude. A source for these strikingpossibly unphysical results may arise from the fact that using a first-order, irreversible kinetics hypothesis is probably a too severe approximation. In partic-

TE D

ular, it seems that an increase of pressure would favor the reversible character of nitridation. Also, applying only a blowing flux correction to the effective nitridation rate instead of simulating the blowing effect with more precision inside the

AC C

EP

porous medium may be another source of error.

41

[28] [30] [29] [29] [27] [26] [26] [31]

RI PT

Nitridation efficiency (-)

ACCEPTED MANUSCRIPT

M AN U

Surface temperature (K)

SC

Helber Carbon 2017

Figure 11: Nitridation wall collision efficiency factor (probability), as compared to literature data.

5. Conclusion & perspectives

TE D

In this work we have developed a methodology aiming at the determination of the local fiber-scale nitridation and nitrogen recombination rate constants from the analysis of experiments carried out in the VKI Plasmatron on a porous carbon preform has been developed. This has involved first to determine the thickness

EP

of the layer that is affected by ablation from tomographic characterizations and image analysis. Then, a simple one-dimensional analytical model has been set

AC C

up and exploited to show how the main geometrical, reaction and transport parameters impact on the observables. An important result is that the thickness of the affected layer varies with the diffusion/reaction rate ratio. Next, an image-

based three-dimensional numerical approach has been set up and validated with respect to the analytical model, bringing a slight update to its results. Inserting experimental values from image analysis and mass loss rate, and data rebuilt by CFD (surface concentration of nitrogen) in the updated analytical and numerical models has led to the identification of local, intrinsic nitridation and recombina42

ACCEPTED MANUSCRIPT

tion rate constants. The low-pressure carbon nitridation rate constants seem in line with past experimental data; at higher pressures, they are markedly lower. On

RI PT

the other hand, rather surprisingly, recombination seems to occur at a higher rate than nitridation.

Several limitations of the present study are: There is a list of limitations in

the present study that should be addressed. First, there have been few measure-

SC

ments of the affected depth, so that the confidence in the values is still moderate. Next, the effect of mechanical erosion has neither been quantified nor accounted

for. The obtained rate constants are therefore overestimates of the actual, purely

M AN U

chemical rate constants. Concerning the chemical phenomena themselves, neither the reversible character of nitridation, including gasification and redeposition steps, nor the possible non-linearity of the rate laws have been considered. Finally, the blowing effect has been incorporated in a rather crude way, as well as the Knudsen diffusion regime. Nevertheless, this work has allowed gaining some

TE D

understanding of the interaction of thean ablative porous medium with ablative atomic gases.Recommendations to deepen and improve the understanding are: This understanding could be pushed further on the basis of a larger number of tests and of samples for each analyzed condition; moreover, the models could be

EP

extended to account for reversible nitridation and for non-linear heterogeneous

AC C

reaction kinetics, for the blowing effect, and for the Knudsen regime. Acknowledgements

The authors wish to thank the European Space Agency for funding the project

AblaCat (AO/1-7664/13/NL/RA), and Dr. Lionel Marraffa (ESA) for helpful discussion. CEA is acknowledged for funding the Ph. D. of C. L. and the Masters internship of X. L.. Dr. Georges Duffa is warmly acknowledged for helpful discussion concerning the analytical model. X-ray µ-CT scans were obtained at the 43

ACCEPTED MANUSCRIPT

Placamat characterization facility of Bordeaux University and CNRS. Dr. Zuheyr Alsalihi is gratefully acknowledged for running the VKI ICP boundary layer sim-

AC C

EP

TE D

M AN U

SC

and Dr. A. Viladegut for his help during the experiments.

RI PT

ulations. P. Collin is acknowledged for his valuable help as Plasmatron operator

44

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

Nomenclature

45

ACCEPTED MANUSCRIPT

Meaning

Unit

A

Dimensionless group (eq.(15a))

-

C

Gas species concentration

mol·m−3

D

Diffusion coefficient

m2 ·s−1

E

Activation energy

J·mol−1

h

Enthalpy

J·kg−1

I1

Prime integral (see eq. (8)

mol·m2 ·s−1

k

Heterogeneous reaction constant

m·s−1

Kn

Knudsen number

-

`

Mean free path

M

Molar mass

m ˙

Mass ablation rate

pd

Dynamic pressure

Pa

ps

Static pressure

Pa

q˙ cw

Cold-wall heat flux

W·m−2

R

Perfect gas constant

8.314 J·mol−1 ·K−1

r

Fiber radius

m

Sh

Sherwood number

-

Ts

MeanAverage surface temperature

K

t

Time

s

vr

Recession rate

m·s−1

z

Position

m

M AN U

SC

RI PT

Symbol

m

kg·mol−1

AC C

EP

TE D

kg·m2 ·s−1

46

Meaning

Unit

α

Reaction order

-

γ

Wall collision efficiency

-

∆

Diffusion/reaction ratio (eq. (13b))

-

δa

Affected depth

m

ε

Pore volume fraction

η

Tortuosity factor

Θ

Thiele modulus

κ

Total/nitridation reaction rate ratio

-

λ

Scaled penetration depth

-

ρ

Mass density

kg·m−3

ρc

Condensation ratio (eq. (13a))

-

ξ

Scaled space coordinate

-

σv

Internal surface area

m−1

τ

Exposure time

s

φ

Solid volume fraction

-

χ

Scaled concentration

-

Ω

Molar volume

m3 ·mol−1

M AN U

SC

Symbol

RI PT

ACCEPTED MANUSCRIPT

-

AC C

EP

TE D

-

47

ACCEPTED MANUSCRIPT

•corr

corrected

•eff

effective

•0

initial value

•b

related to boundary

•i,m

related to ith species in the mixture

•g

related to pure gas

•p

related to porous medium

•s

related to the solid phase

•∞

related to free-stream flow

•nit

related to nitridation

•rec

related to recombination

M AN U

References

Unit

RI PT

Meaning

SC

Symbol

TE D

[1] G. Duffa, Ablative Thermal Protection Systems Modeling, AIAA Education Series, 2013. DOI:10.2514/4.101717. [2] P. Schrooyen, K. Hillewaert, T. E. Magin, P. Chatelain, Fully implicit dis-

EP

continuous galerkin solver to study surface and volume ablation competition in atmospheric entry flows, Int. J. Heat Mass Trans. 103 (2016) 108–124.

AC C

DOI: 10.1016/j.ijheatmasstransfer.2016.07.022.

[3] J. Lachaud, Y. Aspa, G. L. Vignoles, A Brownian motion simulation technique to simulate gasification and its application to ablation, Comput. Mater. Sci. 44 (4) (2009) 1034–1041. DOI:10.1016/j.commatsci.2008.07.015

[4] N. N. Mansour, F. Panerai, A. Martin, D. Y. Parkinson, A. MacDowell, T. Fast, G. Vignoles, J. Lachaud, A New Approach To Light-Weight Ablators Analysis: From Micro-Tomography Measurements to Statistical Anal48

ACCEPTED MANUSCRIPT

ysis and Modeling, no. 2013-2768 in AIAA Papers, American Institute of Aeronautics and Astronautics, 2013, pp. 1–11. DOI:10.2514/6.2013-2768.

RI PT

[5] J. C. Ferguson, F. Panerai, J. Lachaud, A. Martin, S. C. Bailey, N. N. Mansour, Modeling the oxidation of low-density carbon fiber material based on micro-tomography, Carbon 96 (2016) 57–65.

SC

DOI:10.1016/j.carbon.2015.08.113.

[6] J. Lachaud, I. Cozmuta, N. Mansour, Multiscale approach to ablation modeling of phenolic impregnated carbon ablators, J. Spacecraft Rock. 47 (6)

M AN U

(2010) 910–921. DOI:10.2514/1.42681.

[7] J. Lachaud, N. Bertrand, G. L. Vignoles, G. Bourget, F. Rebillat, P. Weisbecker, A theoretical/experimental approach to the intrinsic oxidation reactivities of C/C composites and of their components, Carbon 45 (14) (2007)

TE D

2768–2776. DOI:10.1016/j.carbon.2007.09.034.

[8] B. Bottin, M. Carbonaro, V. V. D. Haegen, S. Paris, Predicted and measured capability of the 1.2 MW Plasmatron regarding re-entry simulation, in: R. A. Harris (Ed.), Proceedings of the 3rd European Symposium on Aerothermo-

EP

dynamics for space vehicles, Vol. SP-426 of ESA Conf. Procs., ESA Publications, Noordwijk, The Netherlands, 1998, pp. 553–560.

AC C

[9] B. Bottin, O. Chazot, M. Carbonaro, V. V. D. Haegen, S. Paris, The VKI Plasmatron characteristics and performance, in: Measurement Techniques for High Enthalpy and Plasma Flows, Vol. RTO-EN-8 of Research and Technology Organization Educational Notes, NATO, Rhode-Saint-Gen`ese (Bel-

gium), 1999, pp. 6–1 – 6–26. [10] B. Helber, O. Chazot, A. Hubin, T. E. Magin, Microstructure and gas-surface interaction studies of a low-density carbon-bonded carbon fiber composite 49

ACCEPTED MANUSCRIPT

in atmospheric entry plasmas, Compos. Part A: Appl. Sci. Manuf. 72 (2015) 96–107. DOI:10.1016/j.compositesa.2015.02.004.

RI PT

[11] B. Helber, A. Turchi, J. B. Scoggins, A. Hubin, T. E. Magin, Experimental investigation of ablation and pyrolysis processes of carbon-phenolic ablators

in atmospheric entry plasmas, Int. J. Heat & Mass Transfer 100 (2016) 810–

SC

824. DOI:10.1016/j.ijheatmasstransfer.2016.04.072.

[12] J. Schindelin, I. Arganda-Carreras, E. Frise, V. Kaynig, M. L. andTobias Pietzsch, S. Preibisch, C. Rueden, S. Saalfeld, B. S. andJean Yves Tinevez,

M AN U

D. J. White, V. Hartenstein, K. Eliceiri, P. Tomancak, A. Cardona, Fiji: an open-source platform for biological-image analysis, Nature Methods 9 (7) (2012) 676—682. DOI:10.1038/nmeth.2019

[13] C. A. Schneider, W. S. Rasband, K. W. Eliceiri, NIH image to ImageJ: 25 years of image analysis, Nature Methods 9 (7) (2012) 671–675.

TE D

DOI:10.1038/nmeth.2089

[14] E. W. Thiele, Relation between catalytic activity and size of particle, Ind. Eng. Chem. 31 (7) (1939) 916–920. DOI:10.1021/ie50355a027.

EP

[15] G. Degrez, D. V. Abeele, P. Barbante, B. Bottin, Numerical simulation of inductively coupled plasma flows under chemical non-equilibrium, In-

AC C

ternational J. Numer. Meth. Heat & Fluid Flow 14 (4) (2004) 538–558. DOI:10.1108/09615530410532286.

[16] D. V. Abeele, G. Degrez, Efficient computational model for inductive plasma flows, AIAA J. 38 (2) (2000) 234–242. DOI:10.2514/2.977.

[17] T. Magin, A model for inductive plasma wind tunnels, Ph.D. thesis, Universit´e libre de Bruxelles, von Karman Institute for Fluid Dynamics, Brussels, Belgium (2004). 50

ACCEPTED MANUSCRIPT

[18] P. Barbante, G. Degrez, G. Sarma, Computation of nonequilibrium high temperature axisymmetric boundary-layer flows, J. Thermophys. Heat Trans.

RI PT

16 (4) (2002) 490–497. DOI:10.2514/2.6723. [19] P. Rini, Analysis of differential diffusion phenomena in high enthalpy flows, with application to thermal protection material testing in ICP facilities, Ph.D.

thesis, Universit´e Libre de Bruxelles, von Karman Institute for Fluid Dynam-

SC

ics, Brussels, Belgium (2006).

[20] A. Turchi, P. M. Congedo, B. Helber, T. E. Magin, Thermochemical

M AN U

ablation modeling forward uncertainty analysis - Part II: Application to plasma wind-tunnel testing, Int. J. Therm. Sci. 118 (2017) 510–517. DOI:10.1016/j.ijthermalsci.2017.04.005.

[21] D. Bianchi, F. Nasuti, R. Paciorri, M. Onofri, Computational Analysis of Hypersonic Flows Including Finite Rate Ablation Thermochemistry, in Procs.

TE D

7th European Workshop on Thermal Protection Systems and Hot Structures, ISSN: 0379-6566, ESA Publications, Noordwijk, The Netherlands (2013). [22] G. L. Vignoles, M. Donias, C. Mulat, C. Germain, J.-F. Delesse, Simplified

EP

marching cubes: An efficient discretization scheme for simulations of deposition/ablation in complex media, Comput. Mater. Sci. 50 (3) (2011) 893 –

AC C

902. DOI:10.1016/j.commatsci.2010.10.027.

[23] G. Vignoles, Y. Aspa, M. Quintard, Modelling of carbon-carbon composite ablation in rocket nozzles, Compos. Sci. Technol. 70 (9) (2010) 1303–1311, DOI:10.1016/j.compscitech.2010.04.002.

[24] G. Vignoles, Modelling binary, Knudsen and transition regime diffusion inside complex porous media, in: G. Battiston (Ed.), Procs. 10th Intl. Conf. on

51

ACCEPTED MANUSCRIPT

Chemical Vapor Deposition, Vol. 5 of J. de Physique IV, EDP Sciences, Les Ulis, France, 1995, pp. 159–166. DOI:10.1051/jphyscol:1995517.

RI PT

[25] G. L. Vignoles, O. Coindreau, A. Ahmadi, D. Bernard, Assessment of geo-

metrical and transport properties of a fibrous C/C composite preform as digitized by X-ray computed micro-tomography. Part II : Heat and gas transport,

SC

J. Mater. Res. 22 (6) (2007) 1537–1550. DOI:10.1557/JMR.2007.0216.

[26] H. W. Goldstein, The reaction of active nitrogen with graphite, J. Phys.

M AN U

Chem. 68 (1) (1964) 39–42. DOI:10.1021/j100783a007.

[27] C. Park, D. Bogdanoff, Shock-tube measurement of nitridation coefficient of solid carbon, J. Thermophys. Heat Trans. 20 (3) (2006) 487–492. DOI:10.2514/1.15743.

[28] T. Suzuki, K. Fujita, K. Ando, T. Sakai, Experimental study of graphite ab-

TE D

lation in nitrogen flow, J. Thermophys. Heat Trans. 22 (3) (2008) 382–389. DOI:10.2514/1.35082.

[29] T. Suzuki, K. Fujita, T. Sakai, Graphite nitridation in lower surface tem-

EP

perature regime, J. Thermophys. Heat Trans. 24 (1) (2010) 212–215. DOI:10.2514/1.43265.

AC C

[30] L. Zhang, D. A. Pejakovic, J. Marschall, M. Dougherty, D. Fletcher, Laboratory investigation of the active nitridation of graphite by atomic nitrogen, J. Thermophys. Heat Trans. 26 (1) (2012) 10–21. DOI:10.2514/1.T3612.

[31] B. Helber, A. Turchi, T. E. Magin, Determination of active nitridation reaction efficiency of graphite in inductively coupled plasma flows, Carbon 125 (2017) 582–594. DOI:10.1016/j.carbon.2017.09.081

52

ACCEPTED MANUSCRIPT

Appendix A. Derivation of eq. (1) Consider a simple array of square cylinders with spacing a and radius r. DurThe porosity ε and solid volume fraction φ are given by : φ=1−ε=π

 r 2 a

2πr a2

M AN U

σv =

(A.1)

SC

The internal surface area σv is given by :

RI PT

ing ablation, the radius will decrease from an initial value r0 to a current value r.

(A.2)

From these two relationships, one can infer: σv =

2φ r

(A.3)

The same relations apply at initial state, i.e.:

 r 2 0

a

√ r0 π a= √ 1 − ε0

EP

therefore :

TE D

φ0 = 1 − ε0 = π

AC C

Inserting eqs. (A.1) and (A.5) in eq. (A.2) gives the desired relationship: √ √ r1 1 − ε 1 − ε0 σv = 2π = 2π √ √ aa π r0 π

53

(A.4)

(A.5)

(A.6)