Damien Campello

Nicolas Tardif

Marwa Moula

Marie-Christine C Baietto

Michel Coret

Jean Desquines

Identification of the steady-state creep behavior of Zircaloy-4 claddings under simulated Loss-Of-Coolant Accident conditions based on a coupled experimental/numerical approach

Keywords: Zircaloy-4, LOCA, Ballooning tests, FEMU

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

After Fukushima accident in March 2011, increasing attention has been given to nuclear reactor core coolability. The full analysis of a LOCA requires complex coupled physical modeling including thermal-hydraulics and mechanical behaviors of fuel rods. During its operating life, the water in the primary loop of a Pressurized Water Reactor (PWR) is pressurized at 155 bar and heated by fuel pellet stacks up to 320 • C. The fuel claddings are also internally pressurized by filling gas and fission gas release. When a breach appears in the loop, the water is depressurized limiting the coolant efficiency. Both the water and the fuel rods are then heated and the claddings are loaded with internal pressure. It can lead to rod ballooning and, potentially, to rod bursting. Re-flooding the nuclear core reactor interrupts this accidental sequence. The fuel rod ballooning can impede the core cooling capacity by flow blockage.

Until the cladding bursts and reaches the nearby fuel rods, a relatively constant plane stress loading is representative of hypothetical LOCA conditions. In the present study, secondary creep data are quantified from Zircaloy-4 (Zr-4) cladding internally pressurized at high temperature (up to 850 • C).

The PERFROI program [START_REF] Repetto | The R&D perfroi project on thermal mechanical and thermal hydraulics behaviors of a fuel rod assembly during a loss of coolant accident[END_REF] aims at improving the thermal-hydraulicmechanical modeling of the entire nuclear reactor core during a postulated LOCA. Several aspects involved in the burst of the LOCA are independently studied in order to further understand the interactions between several fuel rods. The present study focuses on the secondary creep behavior of Stress Relieved Annealed (SRA) Zr-4.

Cladding temperature during a LOCA can reach up to 1200 • C. The steady-state creep behavior of Zr-4 claddings is here investigated for temperatures around the phase transformation from α to (α + β)-phase (ranging from 750 to 850 • C) [START_REF] Forgeron | Experiment and modeling of advanced fuel rod cladding behavior under loca conditions: Alpha-beta phase transformation kinetics and edgar methodology[END_REF]. Cladding micro structure is modified by irradiation under normal operations, however, irradiation damages are annealed early in the ballooning phase of a LOCA [START_REF] Torimaru | Changes in mechanical properties of irradiated zircaloy-2 fuel cladding due to short term annealing[END_REF]. The specific influence of some other parameters such as corrosion and hydriding will not be discussed in this paper.

The secondary creep behavior of as-received Zr-4 is usually [START_REF] Kassner | Fundamental of creep in metals and alloys[END_REF][START_REF] Hayes | Creep of zirconium and zirconium alloys[END_REF] modeled using a power-law detailed in equation [START_REF] Repetto | The R&D perfroi project on thermal mechanical and thermal hydraulics behaviors of a fuel rod assembly during a loss of coolant accident[END_REF]. The temperature T and stress σ influence on the steady-state creep-rates εp ss are respectively defined by the activation energy Q and the Norton exponent n: Material behavior identification is increasingly assisted with FEM [START_REF] Réthoré | Robust identification of elasto-plastic constitutive law parameters from digital images using 3D kinematics[END_REF][START_REF] Avril | Overview of identification methods of mechanical parameters based on full-field measurements[END_REF] for solving inverse problem. FEMU, originally performed by Kavanagh and Clough in [START_REF] Kavanagh | Finite element applications in the characterization of elastic solids[END_REF], relies on an error minimization approach in order to match as fairly as possible observed and the calculated data. The cost function choice is particular to the experimental setup and processing data method [START_REF] Passieux | Multiscale displacement field measurement using digital image correlation: Application to the identification of elastic properties[END_REF][START_REF] Ienny | Identification from measurements of mechanical fields by finite element model updating strategies: a review[END_REF] as much as the constitutive behavior law to be identified.

εp ss (σ, T ) = A • σ n • exp � - Q R.T � , (1)

Material and methods

The specimens wee 90 mm long SRA Zr-4 (see the chemical composition in Table 1) claddings tube samples cut by an electro discharge machining.

Link et al. also reported such measurements in [START_REF] Link | Failure of zircaloy cladding under transverse plane-strain deformation[END_REF].

Test procedure and setup

The test sample was connected onto a 10 kN servo-hydraulic tensile machine using custom grips [START_REF] Tardif | Experimental study of the fracture kinetics of a tubular 16MnNiMo5 steel specimen under biaxial loading at 900 and 1000 • C. application to the rupture of a vessel bottom head during a core meltdown accident in a pressurized water reactor[END_REF]. An enclosure was fitted to the machine to perform the tests under inert environment (see Figure 1). Argon flushing was controlled during the entire test. Tube heating rate was set to Ṫ = 3 • C.s -1 using an induction system controlled by a K-type thermocouple measurement regulation loop (Fig. 2a). The specimen is internally pressurized with argon. A proportional pressure regulator controlled the internal pressurization level. Three pressure levels were stepwise established (Fig. 2b). The pressure-induced bottom-end effect was mechanically compensated with the add of an axial compressive force F , calculated using the initial inner radius R i and the internal pressure P i (see equation ( 3)). The setup components close to the tube are illustrated in Fig. 3. The design of the coil geometry, the induction heating device and the cooling of the grips led to a temperature gradient along with the specimen height. Four cameras were set up over the enclosure circumference as depicted in Fig. 4. The compressive force was measured by a 10 kN load cell, and the pressure by 5 MPa gauge sensors plugged into the upper and lower jaws. Several thermocouples were spot welded in the apparent ROI over the specimen's outer surface.

F = -P i • π • R 2 i ( 3
)

Experiment processing

The data measured by the camera set 1 (Fig. 4) are detailed below. The two key digital image analysis methods are reminded.

2-Dimensional Digital Image Correlation

2D-DIC is performed using the software Ufreckles 4 [START_REF] Réthoré | A fully integrated noise robust strategy for the identification of constitutive laws from digital images[END_REF]. Axial and hoop displacement fields were calculated solving the conservation of the optical flow 4 Developed by J. Réthoré at LaMCoS, INSA de Lyon equation [START_REF] Sutton | Determination of displacements using an improved digital correlation method[END_REF][START_REF] Chu | Applications of digital-imagecorrelation techniques to experimental mechanics[END_REF]. A black undercoating (ULFALUX thermo-coating 1200 • C) was sprayed on the outer tube surface. A white speckle was then applied using BND painting. The speckle size ranged from 10 to 120 μm. Ufreckles uses a global FE-DIC approach. A typical Q4 element size was 30x30 px 2

(1 9 μm). From the gradient tensors, logarithmic strains were computed.

True equivalent creep strain was calculated over the sample surface assuming material incompressibility (equation ( 4)) during high temperature creep.

vp zz + vp θθ + vp rr = 0 ε vp eq = � 2 3 (ε 2 zz + ε 2 θθ + ε 2 rr ) � 1 2 (4)
Typical processed areas are depicted in Fig. 4. They were small quadrant in front of the camera assumed as a plan during the processing. In the case of 2D-DIC, the displacements were extracted on the generator in front of the cameras using the FE interpolation. Systematic errors due to initial and ballooned geometries were estimated lower than 1% of the true equivalent strain up to 0.12 hencky strain [START_REF] Campello | Design and validation of a multi-modal setup for the study of zirconium alloys fuel rods behavior under thermal mechancal conditions simulating a severe accident[END_REF]. Creep curves measured on the surface of a single specimen at five different axial positions are plotted in Fig. 5a.

The steady-state creep-rate axial distributions at three different time steps are plotted in Fig. 5b versus normalized axial position. The different slopes correspond to the 3 internal pressure loadings (Fig. 2b).

Near Infra-Red Thermography

Near Infra-Red Thermography (NIRT) was used as a radiometer [START_REF] Meriaudeau | Temperature imaging and image processing in the steel industry[END_REF] in order to measure the thermal axial distribution. Light intensity I, captured by the silicon sensor of the cameras dedicated to NIRT, is directly linked to the specimen surface radiation. The thermocouples (referenced in Fig. 4) Tc1, Tc2 and Tc3 allowed to determine the radiometric model detailed in equation ( 5), as detailed in [START_REF] Campello | Secondary creep behavior of zr-4 claddings under loca conditions[END_REF]. Due to the speckle, an intensity threshold was performed to process the pixels whose emissivity was the same.

T = K 1 ln( K 2 I + 1) . (5)
About ten profiles were extracted at different test times. The temperature distributions at t 1 , t 2 and t 3 are plotted versus the normalized axial position in Fig. 5c.

Test phenomenology

The specimen heating induces a thermal expansion of the tube. The axial temperature distribution leads to a stiffer cladding material close to the plugs. The internal pressure loading induces (from point 1 to 2 in Figs. 2a and6) a hoop stress σ θθ and a variation of axial stress σ zz across the wall thickness due to bending. Note that the strains are measured over the outer surface of the sample where the stress state is not uniaxial. The combined effect of wall bending and temperature gradient induces non uniform strain rate along the specimen height. The initial tubular geometry tends to balloon, leading to a complex heterogeneous stress state.

Indeed, non-linear geometric effects respectively add compressive (a) and tensile (b) axial stresses on the cladding outer surface near the ROI borders and at its maximum formation location. In addition to this structural effect, the direction of pressure application follows the normal to the inner surface (c). The average axial stress across the wall thickness increases with tube ballooning so as the temperature axial heterogeneity. The hottest zones are located closest to the induction coils.

Finally, thermal and structural effects must be taken into account, relying on a FEMU process, in order to obtain an accurate evaluation of the secondary creep behavior of Zr-4 claddings.

FEMU algorithm 2.4.1. Test modeling

An overview of the different finite element model updatings can be found in Réthoré [START_REF] Réthoré | Robust identification of elasto-plastic constitutive law parameters from digital images using 3D kinematics[END_REF]. The modeling of the boundary conditions affects the whole identification process. Experimental imperfections such as misalignment can be directly measured by DIC. It is thus often better to use the measured displacements as boundary conditions. This is of great importance when the borders are far (with respect to the characteristic dimensions) from the processed zone. Nevertheless, some kinematics assumptions have to be formulated in order to determine the boundary conditions through the thickness in a 3D or axi-symmetric calculation.

In the present study, since the bending of the wall induces non-uniform displacements through the thickness, the whole specimen, up to the clamped ends, was modeled. Please recall that specific hydraulic clamping bushes guarantee a very small misalignment as demonstrated in [START_REF] Tardif | Determination of anisotropy and material hardening for aluminum sheet metal[END_REF] with the very reproducible measurements performed during the tension-torsion test.

Tests were modeled as axisymmetric using the commercial software Abaqus [START_REF] Hibbett | ABAQUS/standard: User's Manual[END_REF]. The part of the specimen located between plugs (58 mm long) was meshed using 6 regular CAX4I elements across its section (see Figure 7c).

Element sizes were respectively 100 and 190 μm in the radial and axial directions.

Radial displacement of the two outer corner nodes and axial displacement of the outer middle node were blocked. The internal pressure was imposed using experimental data onto the inner tube surface.

The temperature axial profile was measured in the ROI using NIRT. Additional thermocouple measurements reported temperatures close to 500 • C at the plug locations. A sine sum model was fitted to these data. At least, three temperature profiles are calculated per internal pressure step. The sine sum models were implemented using the user subroutine UTEMP. A linear time interpolation was performed between these test times. An example of experimental temperature data and the fitted ones are plotted in Fig. 7a versus the axial position. The discrepancy is below 1 • C, as plotted in Fig. 7b.

Elastic behavior was modeled by a temperature-dependent Young's Modulus E(T ) [START_REF] Hagrman | Matpro-version 11: a handbook of materials properties for use in the analysis of light water reactor fuel rod behavior[END_REF] (equation ( 6) with T the temperature in Kelvin) and a Pois- son' ratio ν of 0.342. Thermal expansion was described by the coefficient α th equal to 6•10 -6 •C -1 .

E (T ) = 116.1 • 10 3 -59 • T [M P a] (6)
The creep constitutive law of equation ( 1), still to be determined, was implemented using the CREEP subroutine of the ABAQUS program.

The model relied on two kinds of calculation steps, a static one corresponding to the specimen's heating and a time-dependent viscous one.

The latter one was activated at the end of heating. The non-linear geometric effects were taken into account. Calculations were performed using Abaqus/Standard solver.

Algorithm

The main goal of this study is the determination of the constitutive creep law detailed in equation ( 1), and thus the quantification of A, n and Q in order to be able to predict the creep behavior of Zr-4 under LOCA conditions. FE calculations were performed and the error er (see Fig. 8) between the experimental measured data and the numerical one were computed. The cost function was chosen to only accounts for the 3 steady-state creep-rate axial distributions experimentally determined. The cost function was weighted using the experimental measurements in order to ensure the same influence of the three creep-rate profiles (that do not have the same magnitude). The problem was linearized using the steady-state creep-rate logarithms. The parameters optimization scheme is presented in Fig. 8. The method is based on a Gauss-Newton algorithm. A virtual testing, without noise, was performed for validation purposes.

The error on the parameters was lower than 0.2%.

Results & discussion

The proposed methodology was applied to the test conditions described in Figs. 2 and4. Experimental data resulting from the kinematics and thermal processing are plotted in Figs. 5b and5c.

3.1. Application to a test

Constitutive behavior law parameters

Camera set 1 was first considered.The FEMU strategy led to the determination of the parameters A = e 32,8 M P a -n .s -1 , n = 3.41 and Q = 466 kJ.mol -1 .The creep-rates calculated at timestamps t 1 , t 2 and t 3 of the last FEMU were compared to the DIC measurements in Fig. 9.

The numerically calculated steady-state creep-rates are in good agreement with the 2D-DIC measurements. As reported by Hayes et al.in [START_REF] Hayes | Creep of zirconium and zirconium alloys[END_REF] Zr-4 is well described by Norton-type equations.

Equivalent strains and stresses

The constitutive behavior law was determined using the three steadystate creep-rates axial distributions.

The equivalent creep strains are in good agreement with experiments as illustrated in Fig. 10. The consistency between experiment and FE calculations demonstrates that the cladding behavior is governed by secondary creep only up to 12% strain. The values of the von Mises stress σ vM , observed from FE calculation, at the outer tube surface are plotted in Fig. 11 versus normalized axial position at the three times of interest. As the three steady-state creep-rate calculated and measured curves are consistent, the magnitudes of first and last equivalent stress distributions delineate the validity range in terms of stress of the determined constitutive behavior law. Hence, the set of parameters accounts for the creep behavior from 16 to 30 MPa and from 830 to 850 • C. Note that the duration of the specimen exposure to high temperature is representative of the LOCA conditions, since no thermal constant level was necessary before the application of the mechanical loading.

Axisymmetric assumption

The whole methodology relies on the axisymmetric assumption (modeling, both thermal and kinematics field measurements). We applied a thickness measurement in a r-θ plane at the end of a test as depicted in Fig. 12 in order to validate this assumption. Note that the hencky strain at the end of the third loading reached 0.3 and was thus far above the strain range used for the identification of the constitutive law (at most 0.12).

Analysis of nine tests

The proposed methodology was applied to several testings in order to determine the effect of both stress and temperatures levels on the creep behavior of Zr-4 claddings.

Specimens

Nine experiments were performed under the thermal-mechanical conditions reported in Table 2.

Results

Table 3 summarizes the parameters derived from each test using one set of cameras. Note that the complex specimen preparation and setup potentially induced the loss of some of the thermocouples. The set of cameras was chosen in order to ensure the validity of the NIR thermography measurements and to lower the angular distance between the selected cameras.

Note that the creep behavior can not be modeled systematically using the same constitutive law all over a single experiment gage surface. As illustrated by the parameters obtained for AR-5 specimen, the Norton exponent varied significantly even for a limited stress variation. For an accurate determination of the parameters, the test data were split into two data sets. A first powerlaw was determined using the first two loadings, and then a second finite element model was updated using the last two loading steps. The change of geometry induced by the first internal pressure was taken into account to identify the constitutive parameters of the second step of the test.

The root mean square deviation < S n > detailed in equation ( 7) is also given in Table 3 in order to characterize the accuracy of the modeling. Note that there are about hundred measurement points, distributed along the generator, per internal pressure loading.

< S n >= � � � � 1 n n � i=1 � εfem eq -εexp eq εexp eq � 2 (7)
The steady-state creep-rates are plotted within their domain of validity in Fig. 13. One box corresponds to one identification step. The ranges in terms of stress and temperature were respectively determined using the FE calculation and the NIRT.

The observed consistency between the experiments is noteworthy and speaks not only to the repeatability of the performance, but also to the processing and the modeling of the experiments.

Comparison with literature data

In the literature [START_REF] Kaddour | Experimental determination of creep properties of zirconium alloys together with phase transformation[END_REF][START_REF] Rosinger | The steady-state creep of Zircaloy-4 fuel cladding from 940 to 1873 K[END_REF][START_REF] Busby | High temperature, time-dependent deformation in internally pressurized Zircaloy-4 tubing (LWBR Development Program)[END_REF], a great number of tests is required to identify creep parameters for given (or not) stress values. Furthermore, each test induces its own uncertainty.

The present experimental set up, combined with the measurement techniques, allowed us to identify creep parameters for stress and temperature ranges, and to reduce the number of tests.

The influences of both stress σ vM and the temperature T on the steadystate creep rate εss eq were directly extracted from the present study. The Norton exponents, representing the effect of the stress on the creep-rate and associated to each test, are plotted in Fig. 14.

The temperature T α+β , at which the phase transformation α → α + β begins, is represented in Figs. 14 and 15 using a purple line.

The highest Norton exponent determined above T α+β is 3.41. It was identified using the three internal pressure loadings of the AR-2 experiment so the stress's influence is constant from 16 to 30 MPa in this domain. Garde et al. [START_REF] Garde | Micrograin superplasticity in zircaloy at 850 c[END_REF] performed uniaxial tensile creep tests combined with metallography analysis. They evidenced a grain boundary sliding creep mechanism for [START_REF] Kaddour | Experimental determination of creep properties of zirconium alloys together with phase transformation[END_REF][START_REF] Rosinger | The steady-state creep of Zircaloy-4 fuel cladding from 940 to 1873 K[END_REF] in the α-phase domain.

The activation energies, representing the effect of temperature on steadystate creep-rates, are plotted in Fig. 15. were used to assess the lower stress α-phase constitutive behavior law. The data addressing to σ vM (M P a) ∈ [START_REF] Hibbett | ABAQUS/standard: User's Manual[END_REF]38], and εss eq ≤ 3.10 -4 derived from tests on specimens AR-5,7 and 8 were used to calculate the higher stress α-phase constitutive behavior law. The data of AR-6 and 9 specimens were finally used for modeling the highest creep-rates for temperatures from 780 to 809 • C. The parameters, their validity range, and the root mean square deviation are reported in Table 4.

The deviation < S n > of each of the identification remains lower than 20%. Finally, the activation energies determined at low and high stresses are consistent for low creep-rates. The activation energy determined for εss eq > 3.10 -4 s -1 is higher than the energies calculated independently with each test data. Q and n magnitudes are consistent with those found in the -The setup enables to capture a 20 • C thermal heterogeneity. The thermal distribution perfectly correlates with the strain rate distribution. In addition three different internal pressures were applied during a test. A secondary creep constitutive power-law can be fully identified from a single test.

Hence, the methodology enables significant reduction of the number of tests required to determine the parameters.

-The methodology was applied to nine experiments performed for temperatures ranging from 747 to 848 • C, and internal pressures varying from 1 to 5 MPa. The steady-state creep-rate measurements of the nine tests were in a good agreement between each others. The constitutive creep-law parameters extracted from the FEMU were consistent with those in Ref.

  with A being a third material parameter. Creep mechanisms directly influence n and Q creep parameters. In the α-phase domain, Norton exponents close to 1 are associated with a diffusion creep mechanism whereas for n ∈ [4 -6] a dislocation mechanism prevails [6, 7, 4]. Concerning the activation energy Q, possible creep mechanisms diffusion and dislocation are identified for values ranging respectively from 50 to 150 kJ/mol and 20 to 300 kJ/mol. Few data address the creep behavior in the (α + β)-phase domain, beginning at 810 • C for as-received Zr-4. The experimental setup described in Campello et al. [8, 9, 10, 11] allows to perform uniaxial tests on Zr-4 alloys at a hoop stress ranging from 5 to 55 MPa, for temperature ranging from 750 to 850 • C. The results provided steady-state creep rate and temperature axial distribution along with the specimen. The rich database was then fully analyzed using Finite Element Model Updating (FEMU) in order to identify the creep parameters.

  By coupling high resolution kinematics and thermal fields addressing to heterogeneous thermal-mechanical conditions, the methodology proposed is able to accurately determine secondary creep parameters in a local stress and temperature domain. It also significantly reduces the number of tests needed to map the constitutive parameters according to the thermal-mechanical conditions of interest, and the associated random errors. This paper first details for each test the experimental setup and the data processing providing creep-rate and temperature profiles along the Region Of Interest (ROI) [8]. The test's modeling and the identification algorithm are then detailed, and validated. An experiment performed at 850 • C is itemized, and the modeling options are discussed. Finally, the results of nine experiments are detailed and the influence of temperature and stress are demonstrated in the tested thermal-mechanical conditions. Data from several experiments are correlated in order to identify the parameters of three constitutive creep-laws addressing possible diffusion and dislocation creep mechanisms governing the deformation in the α-phase domain.

Figure 1 :

 1 Figure 1: Schema and upper view of the experimental setup.

Figure 2 :

 2 Figure 2: Thermal and pressure histories during a test performed at 850 • C.

Figure 3 :

 3 Figure 3: The geometries of the sample and the ROI are respectively detailed in a) and b). The test analysis datum is represented in c).

Figure 4 :

 4 Figure 4: Distribution of the cameras over the enclosure.

  Temperature axial profiles.

Figure 5 :

 5 Figure 5: The creep curves at five different axial locations (z is detailed in the legend) are plotted in (a) and creep-rate axial distributions calculated within the test periods t 1 , t 2 and t 3 are plotted in (b) versus the normalized axial position. Temperature axial profiles at the time periods associated to the creep-rate distribution calculation are plotted in (c).

Figure 6 :

 6 Figure 6: Evolution of the specimen during the deformation process (see Fig.2a).

Figure 7 :

 7 Figure 7: (a) Measured and fitted temperature axial distributions, and (b) their difference are plotted versus the normalized axial location. (c) The boundary conditions and the mesh of the modeling are represented.

Figure 8 :

 8 Figure 8: FEMU algorithm

Figure 9 :

 9 Figure 9: Comparison between the experimental and calculation creep-rate profiles based on camera set 1.

Figure 10 :

 10 Figure 10: Equivalent creep strains over time.

Figure 11 :

 11 Figure 11: FE calculated equivalent stress distributions.

Figure 12 :

 12 Figure 12: True hoop strains calculated from the measurements of internal and outer radius.

Figure 13 :

 13 Figure 13: Steady-state creep-rate of all the full testing performed are plotted versus temperature and stress.

Figure 14 :

 14 Figure 14: Norton exponents

Figure 15 :

 15 Figure 15: Activation energies.

3. 3 .

 3 Toward a correlation of the tests in the α-phase domain The maps of the parameters n and Q, based on the identification of constitutive laws using the data of each individual test, illustrated the two creep mechanisms expected in the α-phase domain. For the purpose of severe accident detailed modeling, accounting for multi-physical phenomena, the data of several tests were combined to determine creep laws with validity domain covering a large range of temperature and stress. The data addressing to σ vM (M P a) ∈ [12, 26] obtained from tests on specimens AR-4,5,7 and 8

  Conclusions and perspectivesA ballooning test was designed to investigate as-received SRA Zr-4 cladding secondary-creep behavior under simulated LOCA thermal-mechanical conditions. The specimens exposures duration to high temperature is representative of the LOCA conditions.An axially heterogeneous thermal-mechanical loading was applied on a tubular specimen. Thermal and kinematic fields were measured along a generator. The large amount of data generated by the full field measurements was then used to determine the Zr-4 cladding steady-state creep behavior based on the von Mises criterion by FEMU:

[ 7 , 6 ]

 76 . Norton exponents usually associated to diffusion and dislocation creep mechanisms were determined in the α-phase domain for stresses respectively beyond and above a critical stress of 26 MPa. Constitutive laws addressing each of the mechanisms are proposed for the Zr-4 α-phase creep behavior modeling. The influence of the stress on the creep-rate is modified in the early stage of the two-phases domain. It could possibly correspond to the superplasticity exhibited for Zirconium and its alloys by some authors[START_REF] Garde | Micrograin superplasticity in zircaloy at 850 c[END_REF][START_REF] Zwigl | Transformation superplasticity of zirconium[END_REF] in the (α + β)-domain. This transition in n is supported by the increase of the activation energy observed around 810 • C. The magnitude of Q in the α-phase domain was fully consistent with literature data. Future extensions of this research will focus on the influence of the working environment conditions (hydriding, oxydation) on the creep behavior of Zr-4 cladding under LOCA conditions. The actual running conditions have been selected in order to investigate the influence of combined microstructural, and thermal-mechanical conditions on the creep mechanisms. Results on pre-hydrided specimens will be compared with their as-received twin tests to highlight the influence of hydrogen presence. The proposed identifi-

Table 1 :

 1 Stress Relieved Annealed Zircaloy-4 nominal chemical composition (in wt.%)

	Sn	Fe	Cr	O	C	Zr
	1.29 0.21 0.11 0.13 0.016 bal.
	The cladding outer radius is 4.75 mm, and wall thickness w is 0.57 mm.
	Eccentricities Ξ were evaluated between 1 and 3% using equation (2). This
	rather small eccentricity prevents the specimen from early non-axisymmetric

Table 2 :

 2 Settings of the experiments.

	Specimen Temperature ( • C) Internal pressures (MPa)	FEMU
		[min, max]	P 1	P 2	P 3	identification
	AR-1	[831, 848]	1.04 1.65 2.24	2
	AR-2	[831, 848]	2.17 2.89 3.51	1
	AR-3	[810, 828]	2.75 3.33 3.86	2
	AR-4	[781, 799]	1.48 2.66		1
	AR-5	[778, 798]	2.75 3.33 3.86	2
	AR-6	[780, 800]	4.03 4.52 5.04	2
	AR-7	[747, 767]	3.05 3.70 4.89	2
	AR-8	[755, 775]	2.60 3.35 4.25	2
	AR-9	[785, 810]	3.90 4.30		1

Table 3 :

 3 Parameters resulting from the FEMU method.

	Specimen Loadings	ln(A)	n	Q	< S n >
			A in MPa -n .s -1		kJ/mol
	AR-1	1 & 2	46.6	1.08	543	0.079
	-	2 & 3	47.5	2.15	576	0.051
	AR-2	1 & 3	32.8	3.41	466	0.048
	AR-3	1 & 2	14.2	2.77	283	0.062
	-	2 & 3	17.8	2.89	317	0.062
	AR-4	1 & 2	6.60	1.87	192	0.058
	AR-5	1 & 2	9.48	1.16	198	0.024
	-	2 & 3	-1.33	4.69	206	0.024
	AR-6	1 & 2	-1.19	4.42	192	0.028
	-	2 & 3	-5.18	5.08	178	0.040
	AR-7	1 & 2	3.01	2.52	178	0.0341
	-	2 & 3	-3.45	4.71	186	0.027
	AR-8	1 & 2	0.47	1.27	123	0.032
	-	2 & 3	-7.68	4.39	142	0.062
	AR-9	1 & 2	6.16	3.98	223	0.056

Table 4 :

 4 Parameters identified in the α-phase domain.

	T	σ vM	εss eq	A	n	Q	< S n >
	• C	MPa	s -1	MPa -n .s -1		kJ/mol
	748-798 12-26		e 2.46	1.66	150	0.095
	748-798 26-38 < 3.10 -4	e -6.95	4.95	163	0.094
	780-809 30-47 > 3.10 -4	e 12.2	3.92	293	0.19
	literature describing diffusion and dislocation creep mechanisms of Zr-4 in
	these temperature and stress ranges.

cation method should also be appropriated to identify steady-state creep laws of hydrided specimens. The oxygen influence on the creep behavior of Zr-4 claddings will be further investigated using pre-hydrided and pre-oxidized specimens.

Acknowledgment

The study was performed in the framework of the PERFROI ANR project


Yiran Zhao

Julie Descamps

Brieg Le Corre

Yoan Léger

Alexander Kuhn

Neso Sojic
email: neso.sojic@enscbp.fr

Gabriel Loget
email: gabriel.loget@cnrs.fr

Wireless Anti-Stokes Photo-induced Electrochemiluminescence at Closed Semiconducting Bipolar Electrodes

Keywords: . Photoelectrochemistry, Electrochemiluminescence, Bipolar Electrochemistry, Semiconductors, Silicon

Wireless electrochemical systems constitute a rapidly developing field. Herein, photo-induced electrochemiluminescence (PECL) is studied at Si-based closed bipolar electrodes (BPEs) for designing anti-Stokes systems that can convert IR into visible photons, without direct electrical contact. We show that protection of the anodic emitting pole of the BPE allows the triggering of bright and longstanding emission under the synergetic actions of an external bias and IR illumination. Photoactive nand p-type Si BPEs are studied with front-side and back-side illumination, respectively, and non-photoactive n + -Si BPEs are studied in the dark. Two electrochemiluminescent (ECL) systems ([Ru(bpy) 3 ] 2+ /TPrA and L-012) are tested and we show that the onset bias and the anti-Stokes shift can be controlled by the ECL system that is employed. These advances, rationalized by simulations, will be useful for the design of original PECL systems for chemical sensing or photodetection.

via a mechanism referred to as photo-induced ECL (PECL). [START_REF] Zhao | Electrochemiluminescence with Semiconductor (Nano)Materials[END_REF] Taking advantage of the absorption properties of the different SC materials and emission properties of ECL luminophores, PECL can be employed to convert light in different ways. [START_REF] Laser | Semiconductor Electrodes. Photo-Induced Electrogenerated Chemiluminescence and up-Conversion at Semiconductor Electrodes[END_REF][START_REF] Luttmer | Electrogenerated Chemiluminescence: 34. Photo-Induced Electrogenerated Chemiluminescence and Up-Conversion at Semiconductor Electrodes[END_REF][START_REF] Yu | Luminescence Amplification at BiVO 4 Photoanodes by Photoinduced Electrochemiluminescence[END_REF][START_REF] Yu | Photoinduced Electrochemiluminescence at Nanostructured Hematite Electrodes[END_REF][START_REF] Vogel | Spatiotemporal Control of Electrochemiluminescence Guided by a Visible Light Stimulus[END_REF][START_REF] Zhao | Photoinduced Electrochemiluminescence at Silicon Electrodes in Water[END_REF][START_REF] Zhao | Metal-Insulator-Semiconductor Anodes for Ultrastable and Site-Selective Upconversion Photoinduced Electrochemiluminescence[END_REF] For instance, it can be used to perform anti-Stokes conversion of IR into visible radiation at low applied potentials by using a narrow bandgap photoelectrode (such as Si) and an ECL system emitting visible light. [START_REF] Vogel | Spatiotemporal Control of Electrochemiluminescence Guided by a Visible Light Stimulus[END_REF][START_REF] Zhao | Photoinduced Electrochemiluminescence at Silicon Electrodes in Water[END_REF][START_REF] Zhao | Metal-Insulator-Semiconductor Anodes for Ultrastable and Site-Selective Upconversion Photoinduced Electrochemiluminescence[END_REF] Bipolar electrochemistry is a fundamentally different wireless electrochemical technique, compared to classic electrochemistry, that induces electron-transfer reactions with spatial control on the surface of a conducting object. [START_REF] Fosdick | Bipolar Electrochemistry[END_REF][START_REF] Bouffier | Bipolar Electrochemistry[END_REF] This technique contrasts with the usual 2-or 3electrode setups employed in conventional electrochemical experiments because, in bipolar electrochemistry, the electrode of interest is immersed in the electrolyte without any wire connection. Instead, two external feeder electrodes connected to a power supply are located on either side of the conductive object to generate a potential drop in the electrolyte. At a sufficiently high potential drop, oxidation and reduction take place simultaneously at the poles of the conducting object, referred to as bipolar electrode (BPE). [START_REF] Crooks | Principles of Bipolar Electrochemistry[END_REF] Over the last decade, bipolar electrochemistry has become a unique approach for micro-and nanostructuration, [START_REF] Loget | Bipolar Electrochemistry: From Materials Science to Motion and Beyond[END_REF][START_REF] Bradley | Creating Electrical Contacts between Metal Particles Using Directed Electrochemical Growth[END_REF][START_REF] Loget | Bipolar Anodization Enables the Fabrication of Controlled Arrays of TiO 2 Nanotube Gradients[END_REF] motion generation, [START_REF] Gupta | Designing Tubular Conducting Polymer Actuators for Wireless Electropumping[END_REF] and electroanalysis. [START_REF] Bouffier | Bipolar (Bio)Electroanalysis[END_REF][START_REF] Rahn | Recent Advancements in Bipolar Electrochemical Methods of Analysis[END_REF] Bipolar electrochemistry can be achieved in two main configurations, referred to as open and closed bipolar cells. In open bipolar cells, the BPE is positioned between the feeder electrodes which are immersed in the same electrolyte. In closed bipolar cells, the BPE separates the cell into two compartments, each containing a feeder electrode, allowing for employing different electrolytes. [START_REF] Gamero-Quijano | Closed Bipolar Electrochemistry in a Four-Electrode Configuration[END_REF] Bipolar electrochemistry has been coupled with photoelectrochemistry at SCs for the design of Boolean logic gates [START_REF] Loget | Logic Gates Operated by Bipolar Photoelectrochemical Water Splitting[END_REF][START_REF] Loget | Light-Driven Bipolar Electrochemical Logic Gates with Electrical or Optical Outputs[END_REF][START_REF] Melvin | Wireless Dual Stimuli Actuation of Dye Sensitized Conducting Polymer Hybrids[END_REF] and selective electrodeposition. [START_REF] Ongaro | Asymmetric Modification of TiO 2 Nanofibers with Gold by Electric-Field-Assisted Photochemistry[END_REF][START_REF] Tiewcharoen | Anisotropic Metal Deposition on TiO 2 Particles by Electric-Field-Induced Charge Separation[END_REF] It has been coupled to ECL for sensing, [START_REF] Rahn | Recent Advancements in Bipolar Electrochemical Methods of Analysis[END_REF][START_REF] Zhan | Electrochemical Sensing in Microfluidic Systems Using Electrogenerated Chemiluminescence as a Photonic Reporter of Redox Reactions[END_REF][START_REF] Wu | Electrochemiluminescence on Bipolar Electrodes for Visual Bioanalysis[END_REF][START_REF] Li | Spatially-Resolved Multicolor Bipolar Electrochemiluminescence[END_REF][START_REF] De Poulpiquet | Dual Enzymatic Detection by Bulk Electrogenerated Chemiluminescence[END_REF][START_REF] Wang | Intracellular Wireless Analysis of Single Cells by Bipolar Electrochemiluminescence Confined in a Nanopipette[END_REF] motion tracking, [START_REF] Sentic | Light-Emitting Electrochemical "Swimmers[END_REF][START_REF] Eßmann | Wireless Light-Emitting Electrochemical Rotors[END_REF] and imaging. [START_REF] Villani | Mapping the Distribution of Potential Gradient in Bipolar Electrochemical Systems through Luminol Electrochemiluminescence Imaging[END_REF][START_REF] Anderson | Electrochemiluminescence (ECL)-Based Electrochemical Imaging Using a Massive Array of Bipolar Ultramicroelectrodes[END_REF][START_REF] Ismail | Enhanced Bipolar Electrochemistry at Solid-State Micropores: Demonstration by Wireless Electrochemiluminescence Imaging[END_REF][START_REF] Qin | Label-Free Electrochemiluminescence Imaging of Single-Cell Adhesions by Using Bipolar Nanoelectrode Array[END_REF] In this letter, Si-based photoactive BPEs are investigated for the generation of upconversion PECL in a closed bipolar cell. We show that anodic anti-Stokes PECL can be triggered wirelessly on bipolar electrodes composed of nor p-type Si, leading to the generation of visible photons under IR illumination and the application of an external bias. We used Si as a BPE material, first of all, because of its narrow bandgap (1.1 eV) allowing for anti-Stokes PECL. [START_REF] Zhao | Electrochemiluminescence with Semiconductor (Nano)Materials[END_REF][START_REF] Vogel | Spatiotemporal Control of Electrochemiluminescence Guided by a Visible Light Stimulus[END_REF][START_REF] Zhao | Photoinduced Electrochemiluminescence at Silicon Electrodes in Water[END_REF][START_REF] Zhao | Metal-Insulator-Semiconductor Anodes for Ultrastable and Site-Selective Upconversion Photoinduced Electrochemiluminescence[END_REF] Second, because, depending on its doping, Si can be photoinactive (e.g., with heavy n + -doping) or photoactive, with an electrochemical behavior that can be tuned from photocathodic (with moderate p-doping) to photoanodic (with moderate n-doping). As shown in Scheme 1, two photoactive systems were studied, which differ in the type of doping of the Si BPE: p-Si or n-Si (n + -Si was also employed as a non-photoactive reference surface) and in different illumination configurations (back-side or front-side). In this work, we illuminate the photoactive pole of the BPE (the cathodic/anodic pole for the p-Si/n-Si BPE) so photogenerated minority carriers can readily react with the dissolved species. Back-side illumination refers to the situation where absorption and PECL emission occur at different sides of the BPE, front-side illumination refers to the situation where absorption and PECL emission occur on the same side.

In all cases, the Si BPE (thickness of ~500 µm) separated the two independent electrolyte compartments, each containing a carbon rod feeder electrode connected to an external power supply. The feeder anode compartment was filled with a solution of phosphate buffer saline (PBS) containing 0.2 M [Fe(CN) 6 ] 3-and 0.2 M [Fe(CN) 6 ] 4-(electrolyte 1 in Scheme 1a). The feeder cathode compartment was filled with the ECL electrolyte, composed of 5 mM tris(bipyridine)ruthenium(II) ([Ru(bpy) 3 ] 2+ ) and 0.1 M tri-n-propylamine (TPrA) in PBS (pH 7.4, electrolyte 2 in Scheme 1a), in the first part of this work. Before the experiments, each side of the Si BPE was modified to ensure efficient interfacial charge transfer. The side of the Si surface immersed in the ECL electrolyte (that will be the anodic pole ( + ) under external polarization) was modified with a protective SiO x /Ir coating (see SI for more detail) to avoid anodic decomposition during PECL, as recently reported by our group. [START_REF] Zhao | Metal-Insulator-Semiconductor Anodes for Ultrastable and Site-Selective Upconversion Photoinduced Electrochemiluminescence[END_REF] The side of the Si surface immersed in the [Fe(CN) 6 ] 3-/[Fe(CN) 6 ] 4-electrolyte (that will act as the cathodic pole ( -)) was freshly hydrogenated with dilute HF to generate Si-H and to remove the native insulating SiO x layer, which can otherwise impede charge transfer.

Before studying the bipolar cells of Scheme 1, each interface was investigated in the dark (dashed lines in Figure 1) and under IR illumination (850 nm, solid lines in Figure 1) by voltammetry in a conventional 3-electrodes cell connected to a potentiostat. The left part of Next, the behavior of moderately doped electrodes was studied. As shown in Figure 1b, the response of the p-Si/SiO x /Ir anode is not significantly affected by IR illumination and very similar to that obtained for n + -Si/SiO x /Ir (Figure 1a). This indicates, that, in this regime, the p-Si/SiO x /Ir is forward-biased and charge transfer is controlled by majority charge carriers (h + ). It, therefore, behaves as a conventional anode. In contrast, the left part of Figure 1b shows that the cathodic behavior for p-Si-H is strongly affected by IR illumination. Indeed, this electrode produced a negligible dark current, but a high cathodic photocurrent under illumination. The onset potential for the reduction of [Fe(CN) 6 ] 3-at the illuminated p-Si-H, is more positive than that recorded at the non-photoactive n + -Si-H cathode (Figure 1a), due to the generation of a photovoltage at the interface. This behavior is in good agreement with that of a reverse-biased junction exhibiting rectification in the dark and a charge transfer controlled by photogenerated minority carriers (e -), in other words, a photocathode. n-Si-based electrodes were also studied.

As shown in Figure 1c, in this case, the opposite behavior was observed. Indeed, the n-Si-H cathode sustained cathodic current in the dark and under illumination, and the n-Si/SiO x /Ir photoanode exhibited a strong rectification in the dark and a photoactivity under illumination. In the latter case, the ECL onset potential shifted towards a more negative value compared to the non-photoactive n + -Si/SiO x /Ir (Figure 1a), indicative of conventional photoanode behavior in which the junction is reverse-biased and charge transfer controlled by photogenerated minority carriers (h + ). These experiments showed that p-Si and n-Si-based surfaces, appropriately treated or coated, behave as photocathodes in the [Fe(CN) 6 ] 3-/[Fe(CN) 6 ] 4-electrolyte and photoanodes in the ECL electrolyte, respectively. This is in good agreement with the flatband potential (V fb ) values (Table S1, Figure S1) obtained by capacitance measurements for the four interfaces based on photoactive Si, which confirm that p-Si is in depletion (accumulation) when used as a cathode (anode) and n-Si is in depletion (accumulation) when used as an anode (cathode). An important feature that will be used in the following, is that IR illumination can be used to control the current flow when Si is depleted.

Next, the two interfaces studied in each panel of Figure 1 were combined to make the BPE assemblies depicted in the inset of Figure 2a. In these experiments, the external bias applied between the feeder electrodes was varied and spectra were simultaneously recorded from the   pole of the BPE using a spectrometer. First, the BPE prepared with non-photoactive n + -Si was analyzed in the dark. As shown in Figure 2a, this BPE generated detectable ECL (red emission with a maximum intensity localized at 630 nm) when the applied bias was higher than 1.7 V.

Beyond this, ECL intensity increased with bias, attesting a higher current through the BPE. This value is higher than the threshold required to trigger the redox reactions that can be deduced from the LSV of Figure 1a (~1 V, based on the above-mentioned onset potentials). This can be explained by thermodynamic limitations at the feeder cathode, the requirement of reaching a high current density to generate detectable luminescence, and potential losses at the interfaces.

The two other configurations, based on moderately doped photoactive Si BPEs, are more interesting as they exhibit a strong dependency on illumination conditions. In these experiments, the effect of IR illumination of the photoactive pole of the BPE (  for the n-Si BPE,   for the p-Si BPE, as based on the experiments of Figure 1b,c) was studied.

Figure 2b

presents the results obtained for a BPE prepared from a p-Si substrate. It can be seen here that the BPE is inactive in the dark (empty red circles), due to the rectifying behavior of the reverse-biased p-Si-H/electrolyte junction. However, under back-side IR illumination, the BPE generates PECL at biases higher than 1.4 V (full red circles). This shows that, from that onset bias, electrons photogenerated at the BPE cathodic pole can reduce [Fe(CN) 6 ] 3-and holes can simultaneously induce PECL at the anodic pole. The PECL emission was bright and could be easily visualized with naked eyes, as shown in the inset of Figure 2b. Figure 2c presents the results obtained for a BPE prepared from an n-Si substrate, which has many similarities with the previously described photoactive BPE. Indeed, it is inactive in the dark (empty blue circles), and active under front-side illumination with IR light. Again, PECL is measured at biases higher than 1.4 V (full blue circles), meaning that, above that threshold, photogenerated holes can trigger PECL and electrons can simultaneously reduce [Fe(CN) 6 ] 3-. The PECL (produced on photoactive BPEs) starts at a lower applied bias than ECL (produced on the non-photoactive BPE, Figure 2a), because of photovoltage generation at the photoactive pole of the BPE. Thus, these BPEs are wireless upconverting systems that can be controlled via two stimuli, i.e., the external bias and IR illumination. This is also clearly seen in Figure 2d, which presents the PECL intensity emitted from a p-Si BPE as a function of the IR illumination state (on or off) and the external bias (0 or 1.8 V). Here, we can observe that the BPE emits PECL only when both stimuli are applied and that PECL instantaneously responds to bias and illumination. A general concern regarding such systems is their stability because they are based on Si, an unstable electrode material that tends to deactivate, especially when employed for oxidation reactions. [START_REF] Bae | Strategies for Stable Water Splitting via Protected Photoelectrodes[END_REF] As seen in Figure 2e, n-Si BPEs constantly illuminated and subjected to an external potential of 1.6 V can generate PECL for 22 h, due to the protection provided by the SiO x /Ir coating. [START_REF] Zhao | Metal-Insulator-Semiconductor Anodes for Ultrastable and Site-Selective Upconversion Photoinduced Electrochemiluminescence[END_REF] To further elucidate this behavior, simulations were performed for the n-Si-based BPE, using the parameters obtained with the three-electrodes setup (Table S1, Figure 1) and treating the two phases in contact with Si (electrolyte at the  -pole and Ir at the   pole) as Schottky contacts having a difference in Fermi levels controlled by the external bias. The results, discussed in the supplementary information (Section 2), allowed determining the band diagram of the BPE (Figure S2) and confirmed that the BPE presents rectification, photoactivity, and decreased onset bias due to photovoltage generation. Finally, our BPE systems were assessed with another ECL electrolyte, that is, a 0.1 M NaOH solution containing 3 mM 8-amino-5-chloro-7-phenyl-pyrido [3,4-d]pyridazine-1,4(2H,3H)dione) (L-012). This luminophore is a luminol analog that emits at 490 nm and that has been previously used for Stokes-type PECL. [START_REF] Yu | Luminescence Amplification at BiVO 4 Photoanodes by Photoinduced Electrochemiluminescence[END_REF][START_REF] Yu | Photoinduced Electrochemiluminescence at Nanostructured Hematite Electrodes[END_REF] In addition to its bright ECL emission, a particularity of L-012 is that its oxidation (leading to the ECL emission) occurs at a much lower potential than that of [Ru(bpy) 3 ] 2+ . Thus, we anticipated that L-012 could be used to trigger wireless PECL with BPEs at a lower applied external bias. As shown in Figure 3a, this anodic system was first studied by conventional voltammetry with a 3-electrode setup using pand n-Si/SiO x /Ir anodes. As it was the case previously (Figure 1b), p-Si/SiO x /Ir was not photoactive in this potential range and triggered L-012 oxidation with an onset potential of 250 mV. Conversely, n-Si/SiO x /Ir was photoactive and allowed inducing L-012 oxidation at an onset potential as low as -0.03 V under IR illumination. Again, the considerable negative shift of onset potential (-280 mV) is caused by the photovoltage generated at the n-Si/SiO x /Ir junction. Then, BPE systems were designed as previously, except that the [Ru(bpy) 3 ] 2+ /TPrA ECL electrolyte (electrolyte 2 in Scheme 1) was replaced by the L-012 containing electrolyte. Figure 3b shows the ECL intensity as a function of the external bias for the non-photoactive n + -Si-based BPE studied in the dark. Here, ECL can be detected at a bias >0.4 V, a value which is considerably lower compared to that measured (1.7 V) for the non-photoactive BPE operating with [Ru(bpy) 3 ] 2+ /TPrA (Figure 2a). This is caused by the lower thermodynamic threshold required to drive current through the bipolar cell. Figure 3c and 3d present the results obtained with the photoactive p-Si and n-Si-based BPEs illuminated from the back-side and front-side, respectively. Here, we can see a behavior very similar to that observed previously in Figure 2b,c except that the PECL onset is shifted to a lower onset bias. Indeed, the electrodes exhibited a poor activity in the dark, but a strong PECL emission under IR illumination at biases >0.2 V. As shown in Figure 3e, in this case, the emission spectrum was shifted to lower wavelengths (maximum emission of 490 nm vs 630 nm for Ru(bpy) 3 ] 2+ ). The anti-Stokes shift can be measured as the difference between the wavelength of maximum emission (490 nm, blue spectrum in Figure 3e) and the excitation wavelength (850 nm, brown spectrum in Figure 3e). To the best of our knowledge, this anti-Stokes shift of -360 nm is the largest reported for PECL systems so far. [START_REF] Zhao | Electrochemiluminescence with Semiconductor (Nano)Materials[END_REF] To conclude, we have reported that closed BPEs based on photoactive por n-Si can be used for the wireless generation of anti-Stokes PECL under IR irradiation with bright emissions in the visible spectrum. The protection of the anodic pole of the BPE with a SiO x /Ir coating grants stability to the BPE, which can operate for a long time (> 20 h) under illumination. Both illumination modes, front-side and back-side were demonstrated, and, depending on the ECL electrolyte, the onset bias for PECL emission and the anti-Stokes shift can be both tuned. The back-side illumination mode can be particularly interesting for reducing the noise background as, in this case, all the incident light is absorbed by the BPE. In particular, we have shown that using L-012 as a luminophore allows a shift of 360 nm and an onset bias of only 0.2 V. These results present advances for the manufacturing of PECL systems that can operate only with a light stimulus and can be important for the design of chemical sensors or photodetectors.

ASSOCIATED CONTENT

Supporting Information.

The experimental details, the simulation description and results, supplementary figures, and a supplementary table are provided in supporting information.

Scheme 1 .

 1 Scheme 1. Up-converting bipolar photoelectrodes. a) p-type and b) n-type bipolar photoelectrodes operating with back-side and front-side illuminations, respectively.

Figure 1a -

 1a Figure 1a-c (colored in grey) shows the linear sweep voltammograms (LSVs) of 3 types of freshly prepared Si-H cathodes immersed in the [Fe(CN) 6 ] 3-/[Fe(CN) 6 ] 4-electrolyte. The right

Figure 1 .

 1 Figure 1. Voltammetric study of the junctions with a three-electrode setup. a) LSVs recorded on n + -Si-based electrodes. Inset: schemes of the interfaces studied at the electrodes (the Si surface is the working electrode). b) LSVs recorded on p-Si-based electrodes. c) LSVs recorded on n-Si-based electrodes. LSVs shown in the left (grey square) were obtained at Sibased cathodes in 0.2 M [Fe(CN) 6 ] 3-/0.2 M [Fe(CN) 6 ] 4-in PBS. Right LSVs (white square) were obtained at Si/SiO x /Ir-based anodes in 5 mM [Ru(bpy) 3 ] 2+ and 100 mM TPrA in PBS. Full line-LSVs and dashed line-LSVs were recorded under illumination (850 nm) and in the dark, respectively. Scan rate = 50 mV s -1 .

Figure 2 .

 2 Figure 2. Upconversion PECL at BPEs. a-c) ECL intensity as a function of applied bias for a) an n + -Si-based BPE in the dark, b) a p-Si-based BPE in the dark (empty circles) and under backside illumination at 850 nm (full circle), and c) an n-Si-based BPE in the dark (empty circles) and under front-side illumination at 850 nm (full circle). Left insets in a-c: schemes representing the BPEs and the illumination modes. Right inset in b): photograph showing the PECL emission at the p-Si-based BPE. d) PECL intensity recorded at a p-Si-based BPE when varying the potential and the light inputs. e) PECL intensity as a function of time for an n-Si-based BPE with an applied bias of 1.6 V and a front-side illumination. The electrolytes consist of 0.2 M [Fe(CN) 6 ] 3-/0.2 M [Fe(CN) 6 ] 4-in PBS and 5 mM [Ru(bpy) 3 ] 2+ and 100 mM TPrA in PBS.

Figure 3 .

 3 Figure 3. Low-bias upconversion PECL at BPEs. a) LSVs recorded at n-Si/SiO x /Ir (blue) and p-Si/SiO x /Ir (red) the dark (dashed curve) and under illumination (850 nm). The electrolyte was composed of 3 mM L-012 in 0.1 M NaOH, scan rate = 50 mV s -1 . Inset represents the interface. b-d) ECL intensity as a function of applied bias for b) an n + -Si-based BPE in the dark, c) a p-Si-BPE in the dark (empty circles) and under back-side illumination at 850 nm (full circle), and d) an n-Si-based bipolar photoelectrode in the dark (empty circles) and under back-side illumination at 850 nm (full circle). Inset represents the BPE. e) Normalized spectra of the emission of the LED used in this work (brown spectrum), PECL emission of [Ru(bpy) 3 ] 2+ (red spectrum), and L-012 (blue spectrum). In b,d) The electrolytes consists of 0.2 M [Fe(CN) 6 ] 3-/0.2 M [Fe(CN) 6 ] 4-in PBS and 3 mM L-012 in 0.1 M NaOH.

AUTHOR INFORMATION

Notes

The authors declare no competing financial interests.

ACKNOWLEDGMENT

This work was funded by ANR (LiCORN, ANR-20-CE29-0006).


Michael Döhler
email: michael.doehler@inria.fr

Qinghua Zhang
email: qinghua.zhang@inria.fr

Laurent Mevel
email: laurent.mevel@inria.fr

Kalman predictor subspace residual for mechanical system damage detection

Keywords: Structural health monitoring, damage detection, fault diagnosis, residual design, subspace system identification, vibration analysis

For mechanical system structural health monitoring, a new residual generation method is proposed in this paper, inspired by a recent result on subspace system identification. It improves statistical properties of the existing subspace residual, which has been naturally derived from the standard subspace system identification method. Replacing the monitored system state-space model by the Kalman filter one-step ahead predictor is the key element of the improvement in statistical properties, as originally proposed by Verhaegen and Hansson in the design of a new subspace system identification method.

INTRODUCTION

The detection of damages based on vibration analysis has been widely applied to structural health monitoring (SHM) [START_REF] Farrar | An introduction to structural health monitoring[END_REF]. In this approach, early sign of damages are modeled as changes in the parameters of the underlying mechanical system, then damage detection amounts to parameter change detection. A particular difficulty for such vibration analysis-based SHM is due to the absence of known system inputs, since the structural excitations are usually only caused by natural ambient disturbances. In this sense the considered vibration analysisbased SHM is an output-only monitoring problem.

Among model-based methods for damage detection [START_REF] Carden | Vibration based condition monitoring: a review[END_REF], it is quite natural to perform parameter change detection by comparing model parameters estimated from data collected at different dates. Nevertheless, due to the difficulty for distinguishing the intrinsic characteristics of the monitored structure and the characteristics of unknown natural excitations, both being present in estimated model parameters, it is difficult to automatize such methods. Alternatively, other methods based on direct model-data matching are more suitable for automated damage detection. In this approach, current measurement data are directly confronted to a reference model. For instance, such methods include nonparametric change detection based on novelty detection [START_REF] Worden | Damage detection using outlier analysis[END_REF], whiteness tests on Kalman filter innovations [START_REF] Bernal | Kalman filter damage detection in the presence of changing process and measurement noise[END_REF], and the local asymptotic approach to change detection [START_REF] Benveniste | The asymptotic local approach to change detection and model validation[END_REF].

The two approaches, based on repeated model parameter estimation and on direct model-data matching, are closely related. In particular, for each parameter estimation method based on the minimization of some error criterion, it is quite straightforward to design its counterpart for direct model-data matching following the lo-cal asymptotic approach to change detection [START_REF] Benveniste | The asymptotic local approach to change detection and model validation[END_REF][START_REF] Basseville | Detection of Abrupt Changes -Theory and Application[END_REF]. More recently, a model-data matching method has been inspired by the subspace system identification method, which is a statespace model estimation method with no direct connection to an error criterion minimization, originally developed in the automatic control community [START_REF] Van Overschee | Subspace Identification for Linear Systems[END_REF][START_REF] Ljung | System Identification -Theory for the User[END_REF]. The essential element of this model-data matching method is the design of a residual derived from the subspace system identification method [START_REF] Basseville | Subspace-based fault detection algorithms for vibration monitoring[END_REF]. Based on this residual and its associated statistical evaluation as in [START_REF] Basseville | Subspace-based fault detection algorithms for vibration monitoring[END_REF], extensions have been proposed to consider changing process noise properties [START_REF] Döhler | Subspace-based fault detection robust to changes in the noise covariances[END_REF][START_REF] Döhler | Subspacebased damage detection under changes in the ambient excitation statistics[END_REF] or uncertainties related to the reference model [START_REF] Viefhues | Statistical subspace-based damage detection with estimated reference[END_REF], as well as damage localization and quantification (corresponding to fault isolation and estimation) with the proposed residual [START_REF] Döhler | Fault detection, isolation and quantification from Gaussian residuals with application to structural damage diagnosis[END_REF][START_REF] Allahdadian | Towards robust statistical damage localization Fig. 3. Estimated non-centrality parameter of the test for the 7.5% damage case in spring 2 for different sample lengths N . Fig. 4. Probability of detection of the 7.5% damage case in spring 2 for different sample lengths N . via model-based sensitivity clustering[END_REF]. Successful applications to field data analysis have been reported e.g. in (Döhler et al., 2014a;[START_REF] Gres | A comparison of damage detection methods applied to civil engineering structures[END_REF].

The purpose of the present paper is to propose a new residual related to the subspace system identification method, with an important improvement in its statistical properties, compared to the residual proposed in [START_REF] Basseville | Subspace-based fault detection algorithms for vibration monitoring[END_REF]. This improvement will simplify residual evaluation for damage detection.

In most methods for fault diagnosis, residuals are computed from sensor signals such that their behaviors reflect the occurrence of the monitored faults [START_REF] Blanke | Diagnosis and Fault-Tolerant Control[END_REF]Volosencu, 2018;Escobet et al., 2019). The subspace residual proposed in [START_REF] Basseville | Subspace-based fault detection algorithms for vibration monitoring[END_REF] has a zero mean when the monitored system is in its nominal state, and a non zero mean when some damage appears. These properties are useful for damage detection. However, the computation of this residual mixes noises of different time instants, leading to complex statistical properties. Even under the assumptions that the noises involved in the system are white, the statistical properties of this residual remain complex. The same complexity is also known in the error terms of the subspace system identification method [START_REF] Van Overschee | Subspace Identification for Linear Systems[END_REF].

The main idea in the new residual design of this paper is to replace the state-space model of the monitored system by its Kalman one-step ahead predictor, while the other steps of this design remain similar to those of [START_REF] Basseville | Subspace-based fault detection algorithms for vibration monitoring[END_REF]. This new design has been inspired by a recent subspace system identification method [START_REF] Verhaegen | N2SID: Nuclear norm subspace identification of innovation models[END_REF] which avoids mixing noises of different time instants in its error terms. Though this subspace system identification method relies on a non convex optimization, for the damage detection problem considered in the present paper, the related subspace residual computation remains as simple as in [START_REF] Basseville | Subspace-based fault detection algorithms for vibration monitoring[END_REF]. This paper is organized as follows. The considered damage detection problem is stated in Section 2. The standard subspace residual is recalled in Section 3. The new subspace residual is presented in Section 4. The evaluation of the new residual is considered in Section 5. A simulation example is presented in Section 6. The paper is then concluded by Section 7.

THE CONSIDERED DAMAGE DETECTION PROBLEM

The vibration behavior of mechanical structures subject to unknown ambient excitations can be described by the differential equation

M Ẍ (t) + C Ẋ (t) + KX (t) = f (t),
(1) where t ∈ R denotes the time; M, C, K ∈ R m×m are mass, damping, and stiffness matrices, respectively; X (t) ∈ R m is the displacement vector of the m degrees of freedom of the structure; and f (t) ∈ R m×m is the vector containing the external unmeasured forces, treated as random disturbances. It is assumed that the vibration amplitude is small so that the linear differential equation ( 1) is accurate enough.

Observed at r sensor positions by displacement, velocity or acceleration sensors at discrete time instants t = kτ (with sampling period τ ), system (1) is then described by a discrete-time state space system model [START_REF] Juang | Applied system identification[END_REF]

x k+1 = Ax k + w k y k = Cx k + v k , (2)
where the state vector x k = X (kτ ) T Ẋ (kτ ) T T ∈ R n with n = 2m, the measured output vector y k ∈ R r and the system matrices

A = exp(A c τ ) ∈ R n×n , A c = 0 I -M -1 K -M -1 C , (3)
C = L d -L a M -1 K L v -L a M -1 C ∈ R r×n , (4)
with selection matrices L d , L v , L a ∈ {0, 1} r×m indicating the positions of displacement, velocity or acceleration sensors, respectively. The state noise w k ∈ R n and output noise v k ∈ R r are unmeasured and assumed to be of zeromean, white and Gaussian.

In this paper, the considered damages are changes in the structural stiffness properties of system (1), corresponding to changes related to the parameters of structural elements. A new subspace residual will be designed with an associated residual evaluation method for damage detection.

THE STANDARD SUBSPACE RESIDUAL

Based on the well known subspace method for system identification (Van Overschee and De Moor, 1996), a subspace residual has been proposed in [START_REF] Basseville | Subspace-based fault detection algorithms for vibration monitoring[END_REF] and succesfully applied to mechanical system damage detection [START_REF] Döhler | Subspace-based fault detection robust to changes in the noise covariances[END_REF][START_REF] Döhler | Subspacebased damage detection under changes in the ambient excitation statistics[END_REF][START_REF] Döhler | Fault detection, isolation and quantification from Gaussian residuals with application to structural damage diagnosis[END_REF][START_REF] Viefhues | Statistical subspace-based damage detection with estimated reference[END_REF]. Let us first shortly recall this "natural" subspace residual, in order to motivate the new residual proposed in the next section.

To design a residual based on the state-space model (2), the unknown state vector x k has to be handled somehow. The subspace approach eliminates x k . For this purpose, it follows from (2) that

       y k = Cx k + v k y k+1 = CAx k + Cw k + v k+1
. . .

y k+l = CA l x k + CA l-1 w k + • • • + Cw k+l-1 + v k+l (5)
for any integer l > 0. Then, for some integer s > n,

    y k y k+1
. . .

y k+s-1     =     C CA . . . CA s-1     x k +         v k Cw k + v k+1 . . . C s-1 j=1 A s-1-j w k+j-1 + v k+s-1         . ( 6
)
Given a data set containing collected sensor measurements y 1 , y 2 , . . . , y N , define the Hankel matrices

Y =      y 1 y 2 • • • y N -s+1 y 2 y 3 . . . . . . . . . . . . y s y s+1 • • • y N      , (7)
W =      w 1 w 2 • • • w N -s+1 w 2 w 3 . . . . . . . . . . . . w s w s+1 • • • w N      , (8)
V =      v 1 v 2 • • • v N -s+1 v 2 v 3 . . . . . . . . . . . . v s v s+1 • • • v N      , (9)
the Toeplitz matrix

T =     0 0 • • • 0 C 0 0 . . . . . . CA s-2 • • • C 0     , (10)
the extended observability matrix

O =     C CA . . . CA s-1     , (11)
and the state sequence stack matrix

X = [x 1 x 2 • • • x N -s+1 ] . (12)
Apply ( 6) to each column of Y , then, the results are appropriately arranged as

Y = OX + T W + V.
(13) This equation shows that, if the noises terms T W +V were neglected, the matrix Y ∈ R (sr)×(N -s+1) would have a rank not larger than those of O and X. Moreover, without the noises terms, each column of Y would belong to the subspace spanned by the columns of O. This is the key element behind the subspace system identification method.

Assume that the rank of the extended observability matrix O ∈ R sr×n is n (otherwise the size n of the state vector x k would have to be reduced). Let S ∈ R sr×(sr-n) be the left kernel matrix of O, i.e., a full column rank matrix such that S T O = 0. Given O, such an S is typically computed through the singular value decomposition (SVD) of O. Then, left multiply both sides of ( 13) by S T cancels out the term OX, yielding

S T Y = S T T W + S T V. (14)
In [START_REF] Basseville | Subspace-based fault detection algorithms for vibration monitoring[END_REF], a residual for damage detection has been proposed as follows.

A reference data set is first collected, from which a nominal state-space model is estimated by a subspace identification method [START_REF] Van Overschee | Subspace Identification for Linear Systems[END_REF], yielding an extended observability matrix O 0 and the associated left kernel matrix S 0 , with the subscript "0" indicating the nominal case.

Given a new data set, including a Hankel matrix Y as defined in ( 7) and another Hankel matrix Y -similar to Y but filled with past y k , so that Y -is uncorrelated with the noise matrices W and V defined in ( 8) and ( 9), which contain respectively noises entries w k and v k with k > 0.

Then the subspace residual is defined as

Z 1 √ N S T 0 Y Y - T , ( 15
)
where N is the data sample size, Y -serves as an "instrumental variable" [START_REF] Soderstrom | System Identification[END_REF][START_REF] Benveniste | Nonstationary consistency of subspace methods[END_REF] matrix attenuating the noise terms W and V in Y as expressed in (13), due to the fact that Y - is uncorrelated with W and V .

According to [START_REF] Basseville | Subspace-based fault detection algorithms for vibration monitoring[END_REF], when N → +∞, the limiting distribution of Z is Gaussian, of zero mean if the system remains in the nominal state when the data matrices Y and Y -are collected, or of non-zero mean otherwise. This result is based on a central limit theorem, for which the denominator √ N in (15) is necessary [START_REF] Basseville | Subspace-based fault detection algorithms for vibration monitoring[END_REF].

Now let us have a close view of the term T W in (13).

Each column of this matrix product is expressed in the last term of (6), mixing the noises terms w k of different time index k. Though it is assumed that w k and v k are white noises, the mixing effect of the Toeplitz matrix T "colorizes" the noise terms in (13).

The left hand side of ( 14) may also serve as a residual, because the term OX containing the unknown states OX has been canceled out by S T . However, the remaining noise terms in ( 14) are complex due to mixing effect of the Toeplitz matrix T . Instead of analyzing the complex noises terms, the "instrumental variable" matrix Y -has been introduced, and the limiting distribution of the residual Z is used for hypothesis testing, based on a central limit theorem [START_REF] Basseville | Subspace-based fault detection algorithms for vibration monitoring[END_REF].

Another issue of this subspace residual is that the Toeplitz matrix T and the extended observability matrix O involve powers of the matrix A up to A s-1 . If the considered system has unstable modes, these powers may involve large numerical values.

In order to simplify residual evaluation, the next section will propose a new subspace residual avoiding mixed noise terms. Simpler statistical analysis will become possible, without resorting to a central limit theorem.

THE KALMAN PREDICTOR SUBSPACE RESIDUAL

Now let us recall the stationary Kalman filter applied to the state-space system (2):

xk+1|k = Ax k|k , (16a) e k+1 = y k+1 -C xk+1|k , (16b) xk+1|k+1 = xk+1|k + Ke k+1 , (16c)
where xk+1|k and xk+1|k are respectively the one-step ahead state prediction and the updated state estimation, e k is the Kalman innovation, K is the stationary Kalman gain, usually computed by solving an algebraic Riccati equation.

Combine (16a) and (16c) to yield

xk+1|k = Ax k|k-1 + AK(y k -C xk|k-1 ) (17) = A(I n -KC)x k|k-1 + AKy k , ( 18
)
and rewrite (16b) as

y k = C xk|k-1 + e k . ( 19
)
The following results will be based on the "state-space model"

xk+1|k = A(I n -KC)x k|k-1 + AKy k (20a) y k = C xk|k-1 + e k (20b)
by treating AKy k as an "input" term.

Define Ā A(I n -KC), (21)
Ō     C C Ā . . . C Ās-1     , (22) X x1|0 x2|1 • • • xN-s+1|N-s , (23) T
    0 0 • • • 0 CAK 0 0 . . . . . . C Ās-2 AK • • • CAK 0     , (24)
E      e 1 e 2 • • • e N -s+1
e 2 e 3 . . . . . . . . . . . .

e s e s+1 • • • e N      . ( 25
)
Following the same procedure deriving ( 13) from ( 2), now the following equation is derived from (20):

Y = Ō X + T Y + E, (26)
where Y is still as defined in (7).

Remarkably, the noise term E in (26) as defined in (25) contains simple entries e k , unlike the mixed noise term T W in (13). It is well known that the Kalman innovation sequence e k , as defined in (16b), is a white noise under standard assumptions, as those made in this paper. This simple noise term E in ( 26) is clearly an advantage compared to the terms T W + V in (13), with the Toeplitz matrix T mixing w k as detailed in the last term of (6). The simple noise term E will be helpful for residual evaluation.

For residual generation, like in the case of the previous section, a reference data set is first collected, from which a nominal state-space model is estimated by a subspace identification method [START_REF] Van Overschee | Subspace Identification for Linear Systems[END_REF], and the associated stationary Kalman gain K is computed. The corresponding matrix Ō is then denoted by Ō0 to indicate the nominal case. Accordingly, the left kernel matrix of Ō0 , denoted by S0 , is computed (typically through the SVD of Ō0 ), so that ST 0 Ō0 = 0. A new residual is then defined as

Ξ ST 0 (I sr -T )Y. ( 27
)
If the data matrix Y has been collected when the monitored system remains in its nominal conditions, then (26) holds with Ō = Ō0 . In this case, the analysis form of the new residual writes

Ξ = ST 0 (Y -T Y ) (28) = ST 0 ( Ō0 X + E) (29) = ST 0 E. (30)

EVALUATION OF THE NEW RESIDUAL

The full statistical analysis of the residual Ξ should take into account the Hankel matrix structure of E as detailed in (25). For simplicity, the following residual analysis will be adopted in this work.

Notice in ( 25) that each column of E is correlated with its s -1 neighboring columns at each side, but not with the other columns. The left multiplication of E with ST 0 does not mix the columns of E, hence each column of Ξ = ST 0 E is not correlated to the columns of Ξ other than the s -1 neighboring columns. This property of Ξ implies that "downsampling" Ξ by keeping one column out of every s columns will result in a residual matrix with uncorrelated columns.

Let Ξ s,1 denote the matrix containing the 1st, (s + 1)-th, (2s + 1)-th,. . . , (n c s + 1)-th columns of Ξ, and Ξ s,2 denote the matrix containing the 2nd, (s + 2)-th, (2s + 2)-th,. . . , (n c s + 2)-th columns of Ξ, and so on, with n c = N/s . Similarly, the "downsampled" notations E s,1 , E s,2 , . . . are also defined.

Then vec(E s,1 ) = [e 1 e 2 . . . e snc ] T is the vector containing the innovation sequence until sample s • n c , and vec(Ξ s,1 ) (

I nc ⊗ ST 0 )vec(E s,1 ). ( 31
)
To compute the χ 2 test T s,1 = vec(Ξ s,1 ) T cov(vec(Ξ s,1 )) -1 vec(Ξ s,1 ), (32) we need cov(vec(E s,1 )) = I ncs ⊗ Σ = I nc ⊗ (I s ⊗ Σ), ( 33) where Σ = cov(e k ) is the innovation covariance (computed as part of the Kalman filter), which is invertible. Then

cov(vec(Ξ s,1 )) = I nc ⊗ ST 0 (I s ⊗ Σ) S0 . (34)
Define Σ ST 0 (I s ⊗ Σ) S0 . Since S0 is of full column rank and Σ is invertible, it follows that Σ is invertible. Hence, cov(vec(Ξ s,1 ))

-1 = I nc ⊗ Σ-1 , (35)
and the χ 2 test is computed as T s,1 = vec(Ξ s,1 ) T (I nc ⊗ Σ-1 )vec(Ξ s,1 ). ( 36) Similarly, T s,2 , . . . , T s,s are computed from Ξ s,2 , . . . , Ξ s,s , then

T s 1 s (Ξ s,1 + Ξ s,2 + • • • + Ξ s,s ) (37)
is the test over all the columns of Ξ.

SIMULATION EXAMPLE

In this example, a simulated mass-spring chain with eight elements (Fig. 1) is considered for damage detection. The mass, damping and stiffness matrices of the nominal structural model are defined based on the masses m The test values of the nominal state are used to set up an empirical threshold to decide between hypotheses H 0 and H 1 . For a 1% type I error, the probability of detection for the 7.5% damage is 94.8%, and for the 15% damage the probability of detection is 100%.

1 = m 3 = m 5 = m 7 = 1, m 2 = m 4 = m 6 = m 8 = 2, stiffnesses k 1 = k 3 = k 5 = k 7 = 1000, k 2 = k 4 = k 6 = k 8 =

CONCLUSIONS

In this paper, based on a recent result on subspace system identification, a new residual generation method is proposed in order to improve the statistical properties of the existing subspace residual naturally derived from the standard subspace system identification method. By replacing the monitored system state-space model with the Kalman filter one-step ahead predictor, the new subspace residual avoids mixing noise terms.

Fig. 1 .

 1 Fig. 1. Mass-spring chain with four sensors. stiffness decrease in spring 2, 1000 datasets of length N = 10,000 are simulated and the test statistic T s defined in (37) is computed for s = 10. The test values are shown in the respective histograms in Fig. 2. With N = 10,000, s = 10, model order n = 16 and m = 4 sensors, matrix S0 has m•s-n = 24 columns, and matrix Ξ s has 999 columns. The histograms in the faulty state are well separated with increasing fault size.

Fig. 2 .

 2 Fig. 2. Histogram of χ 2 statistic values for the mass-spring chain in the nominal state, and 7.5% and 15% damage in spring 2. Next, the influence of different sample lengths N on the test are examined. Based on 1000 simulated datasets in the nominal and in the 7.5% damage state for each considered N , the non-centrality parameter of the test is evaluated and shown in Fig. 3. As expected, the noncentrality parameter appears to be linear in N . The resulting probability of detection assuming a 1% type I error for each case is shown in Fig. 4.