Radio

oxidation, then alkaline hydrolysis Production of water-soluble molecules Increase of the radionuclides solubility? An example of an Intermediate Level Long-Lived waste package Considered polymer Poly(ester urethane) (PUR) Constituent of gloves for glove boxes Composed of 3 segments synthesized from these precursors: hard segment extender soft segment + 8.9% inorganic fillers + 1,8% cross linking agents + 0,4% pigments H O O O O O H n 25.5% w 63.4% w Objectives Characterizing and quantifying water-soluble molecules created by the alkaline hydrolysis of the non-irradiated and irradiated PUR at different doses Understanding the degradation mechanisms PUR under radiolysis Irradiated PUR under hydrolysis Identifying the products than can complex the radionuclides Being able to model the complexant release kinetics Commissariat à l'énergie atomique et aux énergies alternatives

-Hydrolysis of each chain occurs with equal probability.

-Water diffusion into the polymer does not limit hydrolysis.

-Water-soluble products are immediatly leached out of the polymer.

-Hydrolysis is a complete reaction.

Involved reactions:

Induced rates:

𝑃 2 + 𝑊 → 2𝑃 1 𝑃 3 + 𝑊 → 𝑃 1 + 𝑃 2 … 𝑃 𝑛 + 𝑊 → 𝑃 𝑟 + 𝑃 𝑛-𝑟 (𝑟 = 1, … , 𝑛 -1) 𝑑[𝑃 1 ] 𝑑𝑡 = 2𝑘 𝐻 [𝑊] 𝑖=2 ∞ [𝑃 𝑖 ] 𝑑[𝑃 2 ] 𝑑𝑡 = -𝑘 𝐻 𝑊 𝑃 2 + 2𝑘 𝐻 [𝑊] 𝑖=3 ∞ [𝑃 𝑖 ] … 𝑑[𝑃 𝑛 ] 𝑑𝑡 = -(𝑛 -1)𝑘 𝐻 [𝑊][𝑃 𝑛 ] + 2𝑘 𝐻 [𝑊] 𝑖=𝑛+1 ∞ [𝑃 𝑖 ]

MODELLING THE SOLUBLE FRACTION RELEASE

Yoon's model

After mathematical geniuses work:

𝑃 1 𝜆 1 = 1 -2 1 - 1 𝜇 𝑛0 𝑒 -𝜏 + 1 - 2 𝜇 𝑛0 𝑒 -2𝜏 𝑃 2 𝜆 1 = 1 - 1 𝜇 𝑛0 𝑒 -𝜏 -2 1 - 2 𝜇 𝑛0 𝑒 -2𝜏 + 1 - 3 𝜇 𝑛0 𝑒 -3𝜏 … 𝑃 𝑘 𝜆 1 = 1 - (𝑘 -1) 𝜇 𝑛0 𝑒 -(𝑘-1)𝜏 -2 1 - 𝑘 𝜇 𝑛0 𝑒 -𝑘𝜏 + 1 - (𝑘 + 1) 𝜇 𝑛0 𝑒 -(𝑘+1)𝜏 with 𝜆 𝑘 = 𝑛=1 ∞ 𝑛 𝑘 [𝑃 𝑛 ], 𝜏 = 0 𝑡 𝑘 𝐻 𝑊 𝑑𝑡, 𝜇 𝑛0 = 𝑖=1 ∞ 𝑖 𝑃 𝑖 𝑖=1 ∞ 𝑃 𝑖 = 𝜆 1 𝜆 0
Considering that m n0 is very high and that only P 1 is water-soluble:

𝑷 𝟏 ≈ 𝝀 𝟏 (𝟏 -𝟐𝒆 -𝒌 𝑯 𝒕 + 𝒆 -𝟐𝒌 𝑯 𝒕 ) -Reaction on surface (versus bulk reaction).

-Water concentration inside the PUR depends on temperature.

There are two competing processes (hydrolysis is the predominant process at " high " temperature):

-Plasticizing by water : PUR more porous.

-Remained synthesis reactants and additives leaching at the first steps of the The ratio between mass loss and TOC confirms that adipic acid is the main product.

Whatever the dose, the release of soluble fraction at 40 and 60°C follows Yoon's model which means that the ester groups mainly hydrolyzed with equal probability and diffusion does not limit hydrolysis.

The release of soluble fraction at 23°C does not follow Yoon's model. Several hypotheses to explain this discrepancy: autocatalysis? Plasticizing? Change in water concentration?

Perspectives: improving the Yoon's model to cover the largest range of temperature.

=> Better understanding of the PUR hydrolysis kinetics. Next step: identifying the degradation products that can complex the radionuclides.

INTRODUCTION

  Nuclear waste managementSeveral countries are considering a deep geological disposal. Switzerland, www.nagra.ch Opanilus clay, 300-900 meters depth France, www.cigéo.com Callovo-Oxfordian clay, 500 meters depth Finland, www.posiva.fi Migmatic gneiss + bentonite clay, 400-450 meters depth INTRODUCTION Polymers in nuclear waste packages Polyvinyl chloride, cellulose, ion-exchange resins, polyethylene… What will become polymers inside the geological disposal? water-soluble molecules created by the alkaline hydrolysis of the non-irradiated and irradiated PUR at different doses Understanding the degradation mechanisms PUR under radiolysis Irradiated PUR under hydrolysis Identifying the products than can complex the radionuclides Being able to model the complexant release kinetics INTRODUCTION MATERIALS: A TWO-STEP PREPARATION 1 st step:PUR is irradiated under air using g rays by LABRA ( 60 Co source), dose rate: ~ 0.9 kGy.h -1 , doses: 500 kGy and 1,000 kGy 2 nd step: non-irradiated and irradiated PUR are then hydrolyzed 3 temperature values: room temperature (≈23°C), 40 and 60°C ∆t is variable and depends on the degradation rate : PUR is irradiated under air using g rays by LABRA ( 60 Co source), dose rate: ~ 0.9 kGy.h -1 , doses: 500 kGy and 1,000 kGy 2 nd step: non-irradiated and irradiated PUR are then hydrolyzed 3 temperature values: room temperature (≈23°C), 40 and 60°C ∆t is variable and depends on the degradation rate : PUR is irradiated under air using g rays by LABRA ( 60 Co source), dose rate: ~ 0.9 kGy.h -1 , doses: 500 kGy and 1,000 kGy 2 nd step: non-irradiated and irradiated PUR are then hydrolyzed 3 temperature values: room temperature (≈23°C), 40 and 60°C ∆t is variable and depends on the degradation rate alkaline hydrolysis solution of 1,000 kGy PUR at day 31 at room T: For more details about mass balances: Fromentin et al. (2016) Polymer Degradation and Stability, 128,

Yoon

  𝝀 𝟏 (𝟏 -𝟐𝒆 -𝒌 𝑯 𝒕 + 𝒆 -𝟐𝒌 𝑯 𝒕 ) 𝑘 𝐻 = 𝐴 × 𝑒 -𝐸 𝑎 𝑅𝑇 MODELLING THE SOLUBLE FRACTION RELEASE Yoon's model applied to our experimental data Room temperature → the model does not match… → Temperature affects the way water-soluble products are released. → It is observed at all doses. explained in Bernstein and Gillen (2010) Polymer Degradation and Stability, 95, 1471-1479. degradation at room temperature is faster than expected. Something speeds up the reaction or the process is different at this temperature value.Hypotheses for this discrepancy: Some of the Yoon's model hypotheses are not valid:

  degradation products i.e. adipic acid.* 𝐷 23°𝐶 < 𝐷 40°𝐶 < 𝐷 60°𝐶 * Salazar et al. (2003) Journal of Polymer Science: Part A: Polymer Chemistry, 41, 1136-1151. and 1,000 kGy irradiated PUR hydrolyzed at room temperature (≈23°C), 40 and 60°C: Mass loss and Total Organic Carbon (TOC) show the same evolution.

temperatures, 3 doses Room temperature ≈ 23°C

40°C 60°C

financed by CEA, AREVA NC, EDF. The authors are grateful to Florence Cochin for her scientific collaboration.


Phillip Finck
email: phillip.finck@inl.gov

Tim Valentine
email: valentinete@ornl.gov

James Dyrda
email: james.dyrda@oecd.org

Jim Gulliford
email: jim.gulliford@oecd.org

Jean-Pascal Hudelot
email: jean-pascal.hudelot@cea.fr

Upendra Rohatgi
email: rohatgi@bnl.gov

OECD-NEA EXPERT GROUP ON MULTI-PHYSICS EXPERIMENTAL DATA, BENCHMARKS AND VALIDATION

Keywords: multi-physics, benchmark, validation, experiment, modeling, simulation

Coupled multi-physics computational methods continue to evolve to meet the needs of the R&D community designers, operators and safety regulators, in order to improve predictive accuracy and precision and to evaluate complex operational or accidental scenarios. Novel multi-physics simulation tools are designed to enable rigorous modelling of coupled behaviors between, inter alia, reactor physics, thermal-hydraulics, fuel performance and coolant chemistry. In order to be used to their full potential, these tools will require a more complex array of validation tests due to the multiple length and time scales, as well as the number of physical phenomena being simulated. However, the ability to conduct appropriate validation experiments has either progressed very little or in some cases significantly regressed. Recognition of this divide has led research and industry experts from across the NEA nuclear science community to form a new Expert Group on Multi-physics Experimental Data, Benchmarks and Validation (EGMPEBV)

* . The paper describes the over-arching objectives of the group, in relation to the original motivations and eventual output goals. The specific tasks to be completed under the direction of the two established Task Forces are then presented alongside a depiction their inter-dependence.

INTRODUCTION

Computational analysis methods continue to evolve in many nuclear power countries to meet the needs of the R&D community, designers, operators and safety regulators, to improve predictive accuracy and precision and to evaluate complex operational or accidental scenarios that could have only been addressed by experimental means or simplified bounding calculations in the past. As a result of these developments, computational methods targeted at multi-physics and multi-scale simulations are beginning to be used both for deeper conventional analyses and for new applications. The goal of these codes is to allow modelling of highly complex scenarios at a very high level of spatial, phenomenological and/or temporal resolution and with demonstrated high accuracy. They also seek to enable rigorous modelling of coupled behaviours between, inter alia, reactor physics, thermal-hydraulics, fuel performance and coolant chemistry.

The increasing resolution of such analytical tools does not, however, eliminate the need for suitable validation via comparison to experiments. On the contrary, in order to be used to their full potential, these tools will likely require a more complete array of validation tests as a result of multiple length and time scales, as well as the number of physical phenomena being simulated. They will also require accurate measurements of all terms, including coupling terms, in situations where the appropriate experimental techniques and facilities might not exist today. However, the ability to conduct validation experiments has for some applications progressed at a slower pace than for computational methods, or in some cases it has significantly regressed (through shutdown of facilities, or through retirement of experts); furthermore, as modelling capabilities are reaching deeper levels in single areas and in coupled behaviours, appropriate experimental techniques simply do not yet exist. Recognition of this divide has led research and industry experts from across the NEA nuclear science community to form a new Expert Group on Multi-physics Experimental Data, Benchmarks and Validation (EGMPEBV). The aim of the group is to provide member countries with guidelines and recommendations for validating and improving their novel multi-physics simulations, and access to key experimental data.

MOTIVATION, OBJECTIVES AND GOALS

Validation of multiple physics models requires a wide variety of experimental data, which emphasises the importance of maximising the use of historically accumulated data to avoid the significant cost of performing similar experiments today. The preservation, evaluation and dissemination of such legacy validation data represent a cost-effective path forward to validate modern codes. Identification and prioritisation of key legacy data relevant to modern requirements is therefore one of the group's primary goals. A review and evaluation of the data by experts, both current and contemporary to the experiments in question, will also be undertaken. The target end product will be evaluated and independently reviewed benchmark datasets with quantified uncertainties, which are of significantly greater value to users than the "raw" documentation. This will follow the process established by OECD-NEA in, for example, the ICSBEP and IRPhE projects. The EGMPEBV will also build upon efforts being made by other Expert Groups, such as the Expert Group on Uncertainty Analysis in Modelling (EGUAM), to develop methodologies and recommendations for uncertainty propagation.

One such set of experimental data is from the former Loss-of-Fluid Test (LOFT) facility in the United States where an NEA international project was conducted between 1983 and 1989 † . LOFT data are particularly unique in that they originate from the only large-scale multiple-phenomena test facility that employed a nuclear powered core. Represented members from the United States have therefore made it a priority to collect, analyse and re-model these data so as to re-evaluate the uncertainties and the sensitivity of the measured data for use in multi-physics benchmarks. Support for this activity is being provided by the newly established U.S. Department of Energy's Nuclear Energy Knowledge and Validation Center (NEKVAC), which will work in close collaboration with the Expert Group and other similar national initiatives.

At the same time, while historical data are undoubtedly of great value, they may be limited by various factors; these usually relate to past experiments being targeted at the validation of older codes and the application of those codes. As a result, these historical data can often exhibit limitations including a lack of measurement accuracy or lack of knowledge of the measurement uncertainty in the experimental techniques; an absence of measurements of some critical parameters; and, an obscuring of relationships between individual physics phenomena owing to the output data resolution or simply the inability to allow for a detailed reinterpretation of the experiments based on the documented information.

Part of the role of the EGMPEBV will therefore be to help identify gaps in experimental data, where the scarcity of information is detrimental to the validation efforts of stakeholders. By comparing similar needs, efforts to fill such gaps may be co-ordinated across member countries. Where data does not exist, suitable experiments, facilities and measurement techniques may be proposed and developed to address those specific needs. A significant benefit of this international effort would be the leveraging of experimental capabilities that are likely to go beyond the capacity of any single country to implement to achieve the desired results -a true representation of cost efficiency for all partners involved.

Establishing consensus guidelines for the application of validation data is essential in light of the developing multi-physics code systems. The EGMPEBV will aim to establish standards for evaluating experimental data and determining how these data should be applied to the specific codes and applications in question. The resulting output, along with appropriate phenomena identification and ranking tables (PIRT), could help guide users on the applicability and importance of particular experimental data to reactor phenomena or scenarios of interest.

Finally, guidance for performing robust validation analyses is also important, including for methods to extrapolate uncertainties beyond the validation domain and for estimates of the degree to which stakeholders can rely on the results. This is closely related to the issue of scaling, or when factors must be applied to experimental results because of the extrapolation from the experimental configuration to a full-size application. To compile and make available such recommendations will thus be of significant benefit to the scientific community.

ORGANIZATION AND TASK FORCE ACTIVITES

The organization of the EGMPEBV relies on three Task Forces, focused on 1) Experimental Data, 2) Methods and Standards, and 3) Specific Applications (note however that the third Task Force is not yet fully constituted). The major aspects and expectations of Task Forces one and two are summarized below.

Task Force 1 is focused on experimental data qualification and benchmark evaluation. It aims at providing better and more accurate experimental datasets which support validation of high-fidelity multi-physics modelling and simulation (M&S) tools. At the same time it has to consider the incompleteness of past experimental data sets (in terms of data, documentation or uncertainties) and with the limited number of available multi-physics experimental facilities. Hence, the three main objectives of Task Force 1 are the following:

-reviewing and re-evaluating past experimental datasets for the validation of traditional or novel multi-physics codes, -defining the needs and priorities for new experiments (e.g. pellet-clad interaction, cladding integrity, CRUD-Induced Power Shift and CRUD-Induced Localized Corrosion…), and -establishing and recommending processes for designing new experiments using improved instrumentation, experimental techniques and data treatments with the intent to provide high quality data along with detailed uncertainty information, specifically acquired for validating high fidelity multi-physics M&S tools.

Task Force 2 was designed to develop validation methods and guidelines, and uncertainty qualification for the new family of multi-physics, multi-scale codes in the context of emerging demands such as longer fuel cycles and power uprate.

The capability to validate these codes has made significant progress. Modern validation, verification and uncertainty quantification (VVUQ) techniques enable analysts to extract information from existing experiments in a systematic manner and provide the users with a quantified uncertainty estimate. Currently, there are efforts in the U.S. (e.g. VERA within CASL [START_REF] Turinsky | Update on Capabilities Development at CASL[END_REF], SHARP [START_REF] Siegel | Software design of SHARP[END_REF] and MAM-MOTH [START_REF] Gaston | MOOSE: A parallel computational framework for coupled systems of nonlinear equations[END_REF] within NEAMS) and NURESAFE [START_REF] Chanaron | Advanced multi-physics simulation for reactor safety in the framework of the NURESAFE project[END_REF] in the EU to develop multi-physics/multi-scale analyses and associated validation techniques. These new approaches require evaluation of existing tests for appropriate data or to perform new tests to fill the gaps in required test data.

The principle objectives of Task Force 2 initiative are:

-development of consensus guidelines for validation of multi-physics M&S tools and data, -development of guidelines for performing uncertainty qualification and evaluating ranges of applicability for predicting M&S performance outside of the validation domain, and -identifying needs for specific experiments with the intended purpose of validating multi-physics M&S tools and data.

In the following sub-sections, each of the specific tasks planned under the two Task Forces, in order to meet the stated objectives, are described in detail. The links and inter-dependencies between these tasks are also depicted in Figure 1.

Task Force 1 3.1.1 Task 1: Define Scope of Multi-Physics Applications

The objective of this activity is to define the scope of the MPEBV activity as it relates to other activities within the OECD-NEA as well as to define the terminology that will be used to categorize the types of multi-physics benchmarks and simulation processes. Other Expert Groups under the auspices of the OECD-NEA focus on validation and benchmarking of single-phenomena physics (ICSBEP, IRPhE, etc.) or in some instances dual-phenomena physics (SINBAD, SFCOMPO, etc.). The intent of this Expert Group is on multi-physics M&S. As such, the Expert Group will establish the processes for engaging with the other Expert Groups (such as that on uncertainty analysis in modelling, EGUAM) and Task Forces under the auspices of the NEA. Furthermore, the Expert Group will engage with national validation centers and modeling and simulation development efforts (CASL, NURESAFE, etc.) within the NEA member countries.

Various groups utilize different nomenclature when describing multi-physics M&S as well as characterization of experimental data and types. However, the semantics employed to describe these areas differ among those focused on nuclear safety as compared to those focused on the nuclear science areas. In order to provide some consistency and coordination among these focus areas, the EGMPEBV activity would classify the experimental datasets into three separate areas:

-separate effect tests (SET), -multiple effect tests (MET), -plant measurements and observations (PMO).

The choice of these designations is to provide consistency with existing nomenclature employed by much of the nuclear community as well as to minimize confusion when differences exists between those in the nuclear safety fields and those in the nuclear science fields. This task will also provide concise and clear definition on the differences between novel (N) and traditional (T) methods, what is meant by "coupling", how uncertainty methods have been used for sensitivity and uncertainty analyses, etc., a description and definition of the scope of traditional (T) multi-physics applications will be developed along with a description and definition of the scope of novel (N) multi-physics applications.

Task 2: Summary Reports on the Current Status and Expected Needs for Validation of Multi-Physics Modeling and Simulation Tools

One of the objectives of this task is to develop a report that summarizes the use of best estimate plus uncertainty (BEPU) evaluations in industry and regulation based on traditional (T) multi-physics calculations supplemented by sensitivity analysis (uncertainty quantification (UQ)). To reduce excessive conservatism associated with the design of some operational and safety systems, industry representatives and regulators have started using best-estimate analyses coupled with probabilistic risk assessment to evaluate plant operation under normal and off-normal conditions in lieu of the use of worst case scenario approaches. The move toward best-estimate analyses (and thus toward enhanced uncertainties analysis) in lieu of bounding analyses has resulted in the need to better model the physical phenomena in nuclear reactors. This report will summarize the approaches for benchmarking and validation of traditional multi-physics modeling and simulation tools.

A second objective of this task is to develop a report that will summarize the current approaches to migrate from traditional (best-estimate plus uncertainties -T/BEPU) to novel (high-fidelity first principles with embedded uncertainty quantification -N/UQ) multi-physics modeling and simulations. This report will include a summary of the activities in countries or multi-national approaches such as CASL, SHARP, MAMMOTH, NURESAFE etc. for validation of specific problems.

The final objective of this task is to develop a summary report on the availability of experimental data that will be needed to validate both traditional and novel modeling and simulation (M&S) codes. This report should describe the following: a. Available reactor physics data and benchmarks; b. Available fuel modeling data and benchmarks; c. Available core thermal-hydraulics data and benchmarks; d. Available system thermal-hydraulics data and benchmarks; e. Available traditional (T) multi-physics data and benchmarks; and f. Anticipated needs for validation of traditional (T) and novel (N) multi-physics tools.

Task 3: Summary Report on the Major Challenges and Priorities for Validation of Multi-Physics

Modeling and Simulation Tools

Some of the key issues facing the nuclear industry include the aging of the nuclear fleets throughout the world, the need to extract more from the nuclear fuel with higher burn-up, the desire to increase the power output of existing reactors, the need to develop alternate accident-tolerant fuels, and the need for better understanding of normal and off-normal operating conditions, to name just a few. Several operational challenges for nuclear power plants such as crud-induced power shifts, crud-induced localized corrosion, pellet-cladding interactions, and grid-to-rod fretting cannot be adequately addressed using conventional M&S tools and analyses; hence, conservative estimates of the impacts of these phenomena on reactor operations are often assumed to ensure the safe operation of nuclear power plants. In addition, significant safety issues such as departure from nucleate boiling, loss of coolant accidents and reactivity initiated accidents are often analyzed using rather conservative estimates of the likelihood of such events as an added precaution. Such excess conservatism may lead to less than optimal operation of the plants.

This task will produce a report that addresses the uncertainty treatments in the validation process for these challenge problems and the potential benefits to reducing the excess conservatism. What is the impact of best estimates as compared to high fidelity simulations as compared to other as yet to be defined approaches? The report should address the views of both the vendor/operator that desires to make the optimal use of a plant/facility as compared to the desire of a regulator to ensure the safety of the plant/facility (performance vs. safety). This report should also address the differences in the challenges/problems that might arise when considering different reactor design (Gen. II&III, Gen. IV, etc.). The EGMPEBV focus is on currently operated Light Water Reactors (LWRs) from Gener-ation II, next to be built Generation III and Generation III+ LWRs as well as Small Modular Reactors (SMRs) and potentially representatives of Generation IV -High Temperature Reactors (HTRs) and Fast Reactors (FRs).

The part of the report discussing the challenges and priorities for validation of traditional (T) multi-physics modeling and simulation tools will address: a. Using the existing data; b. Creating new data -priorities and mechanisms; c. Including uncertainty quantification in traditional (T) multi-physics simulations (BEPU methodologies) and validation of uncertainty quantification methods: status, needs and priorities.

The experience and expertise in uncertainty propagation in traditional (T) multi-physics calculations and applications of BEPU methodologies will be summarized including experiences by countries as well as large international projects such as the OECD-NEA LWR Uncertainty Analysis in Modeling (UAM) and International Atomic Energy Agency (IAEA) HTR UAM Coordinated Research Program (CRP).

Task 4: Recommendations and Implementation of Processes for Evaluating Existing Experimental Data for Multi-Physics Modeling and Simulation

The objective of this activity is to define and implement the processes for evaluating past experimental data sets and validation efforts including the quantification of uncertainty in experimental data, the ranking of important physical phenomena, the identification of important measured parameters, the processes for converting raw experimental data to measured parameters, the processes for converting the measured parameters to modeled parameters if they differ, etc. Furthermore, guidance will need to be provided on the level of acceptability of uncertainty data for those experiments that are to be considered benchmark quality and those experiments that will only be considered for general purpose testing. Special attention has also to be given to the completeness and the representativeness of existing experimental data sets.

Consideration will be given to the conversion of "measured" parameters to parameters that can be directly simulated using multi-physics methods. The analysis should include the conversion of the "raw" data to the "measured" data. This should also include the treatment of uncertainties from the conversion process and also any uncertainties from the approximation methods employed in the simulations. Focus will be given to the existing data with uncertainties from nuclear power plants operation and tests i.e. plant measurements and observations (PMO). The implementation of this activity will result in evaluations of experimental data provided by participants form the various member countries. This activity is closely related to task 7 that involves the development of the benchmark models from the evaluated experiments. This will entail establishing a peer review process and database systems for collection and retention of the evaluated experiments. The resulting evaluations and databases will be shared among all member countries of the Expert Group. The availability of new facilities for conducting multi-physics experiments is of a concern. No country among the NEA membership maintains or operates a full complement of facilities that would be needed to study and evaluate all of the possible multi-physics phenomena and the coupling of these phenomena. The cost for such facilities often exceeds that which could be borne by a single country. As such, collaborative arrangements need to be considered for optimal use of limited resources among the NEA member countries.

The need to have access to prototypical experimental conditions with a complete and good quality experimental dataset (hence a complete and adapted set of instrumentations and experimental techniques) should also be addressed. A key outcome of this task will be to review and update the NEA report on "The Research and Test Facilities Required in Nuclear Science and Technology". This updated review should include a comprehensive review of the experimental facilities in Russia be-cause of the breadth of facilities and capabilities that have not before been catalogued by the NEA. This report should address the challenges with proprietary data sharing and how this could limit the validation of multi-physics M&S tools and data.

The report will address also what data needs to be collected from nuclear power plants, especially from the new plants during their commissioning tests and subsequent operation. The new plants have better measurement capabilities, as demonstrated by the OECD VVER-1000 Kalinin-3 Coupled Code Benchmark, where high-quality multi-physics data was collected during a commission test at the start-up of Kalinin-3 unit, and subsequently used to develop an international benchmark. Similar opportunities exist with the commissioning of AP1000, EPR-1600 and other new plants worldwide. A report will also be generated to identify the high-priority needs for the broader nuclear industry. A critical activity of this group will be to identify mechanisms to foster international collaboration in fulfilling the needs for those high-priority experiments identified in the summary report, including a recommendation for an initial demonstration of such a collaborative effort.

Task 6: Summary Report on Traditional and Novel Measurement Methods for Validation of Multi-Physics Modeling and Simulation Tools

The objective of this activity is to develop a report on traditional measurement methods and the limitations of those methods that could be improved to be more representative of the parameters that can be simulated. Identify "ideal" improvements for some measurements that would provide more fundamental parameters for direct comparison with simulations. There is also a need to address the evolution of measurement methods using modern tools as compared to those used in past experimental facilities (novel vs. traditional measurement methods).

Task 7: Summary Report and Implementation of the Guidance for Developing Multi-Physics

Benchmarks

The objective of this activity is to define the process and methodologies for evaluating, documenting, and utilizing multi-physics benchmarks. This effort will build from the experience of the NEA in the development of the ICSBEP and IRPhE benchmarks as well as the traditional (T) multi-physics benchmarks using experimental data such as OECD Ringhals BWR Stability benchmark, OECD/NRC BWR TT benchmarks, OECD/DOE/CEA V1000 CT benchmark, OECD Kalinin-3 Coupled Code benchmark, OECD/NRC Oskarshamn-2 BWR Stability benchmark and the OECD LWR UAM benchmarks. The established procedures within OECD-NEA for developing benchmarks for traditional (T) multi-physics tools will be summarized. The guidance, developed by Task Force 2, will include the requirements and processes for describing and documenting the physical characteristics of multi-physics benchmarks, the requirements for quantifying uncertainties in the physical characteristics and the measured parameters, the processes for equating the actual measured parameters and their associated uncertainties to modeled parameters, etc. The implementation of the guidance will result in benchmark evaluations that will be archived and distributed to all member country participants in the Expert Group.

Task Force 2 3.1.1 Task 3: Summary Report on the Existing Practices for Multi-physics Validation

Code developers, research organizations, and professional societies have developed various practices for validation of M&S tools and data. The approaches may vary depending on the limitations of the M&S tools and data as well as the intended application of such M&S tools. A report will be developed that summarizes some of the primary practices, including recommendations and guidelines that have been developed within the member countries of the NEA. Some of the key factors that will be included in the summary report are as follows:

1 The summary report will identify several commonly used assessment frameworks and summarize the general approach, the strengths, the weaknesses and ease of use of these frameworks. The report will also include an overview of the PIRT process along with the key issues that must be considered when identifying relevant physical phenomena and the coupling among the phenomena such as thermal-hydraulics and neutronics. The report may include a summary of the principle issues to consider when coupling physical phenomena that have varying degrees of fidelity in the temporal, spatial and/or energy domains. Likewise, the report will summarize current approaches for estimating total uncertainty at various levels of the validation hierarchy and the cautions and challenges in such estimates. The report will also provide recommendations as to the treatment of experimental uncertainty and the inference of measured data with respect to simulated responses of interest. In some instances, the underlying physical phenomena are approximated in the simulations and engineering approximations may be used to model observed behavior using rudimentary models that rely on "calibrations." The report will include recommendations for the use of "calibration" and the propagation of uncertainty from this use. A critical topic in estimating the total uncertainty at the application conditions of interest is how to estimate the model form uncertainty using validation metrics. Finally, the report will include a summary of the primary recommendations for the use of multi-physics modeling and simulation tools extrapolated beyond the validation domain.

Task 4: Examples of Implementing Validation Processes for Novel Methods

The objective of this task is to develop a summary report on the current progress of individual member countries or multi-national organizations in the application of validation approaches for specific problems. Member countries will be encouraged to annually submit reports to the Task Force that describe the validation of advanced M&S tools and the approaches employed therein. The CASL program in the U.S. and the NURESAFE program in the European Union have already made significant progress in validating advanced M&S tools for specific problems, which will be compiled by the Task Force in an annually updated report.

Task 5: Development of Validation Matrices for Specific Problems

The nuclear industry is required to provide sufficient evidence regarding the predictability of reactor behavior under certain situations and scenarios. The analysis of these off-normal events often entails the use and validation of M&S tools to ascertain the degree of certainty regarding the predictability of a reactor under such off-normal conditions. Furthermore, the nuclear industry has identified a few specific "challenge" problems in which the use of advanced M&S tools and data could be used to reduce excess conservatism or improve predictability of safety concerns. These challenge problems are likely to be dependent on the reactor type and operational constraints. As such, the Task Force will undertake to develop validation matrices for "standardized" scenarios and challenge problems such as, for example, departure from nucleate boiling, crud induced power shifts, fuel-cladding and cladding-coolant interactions. Member countries will be asked to identify specific standardized scenarios or challenge problems and then to develop a validation matrix for the specific problem by identifying the validation data for the various physical phenomena and at different levels of the validation hierarchy. The Task Force will define a standardized format for the validation matrix and thereafter annually compile the member country contributions into a summary report to be submitted to the NEA.

The development of the validation matrices will serve as the starting point for identifying gaps in the validation data sets and or identifying the availability of data sets that could be shared among member countries. The Task Force will review the member country contributions provided under task 5 and identify gaps in validation data for the same problem. Additionally, the comparison among the member country submissions for the same problem could serve to identify data sets that might be shared among the member countries. Following each update, the Task Force will develop a summary report that identifies gaps in the validation matrices for specific problems. If appropriate, the Task Force may provide recommendations for new experiments or for the sharing of existing experimental data.

CONCLUSIONS

It is clear that the planned objectives of the EGMPEBV are extremely challenging and will require significant effort to achieve. However, the participants have already committed to a demanding schedule of work and have achieved good progress on a number of the initial tasks. Furthermore, additional contributors to the group are still being sought in order to provide the diverse expertise necessary to achieve the scope of objectives described. Despite the breadth of the task at-hand, the ultimate aim is clearly understood by the members, which is for the nuclear industry to be able to realize the potential benefits of novel methods over traditional ones.

Finally, it is noteworthy that a third Task Force is currently under formation to address specific applications of multi-physics methods, through application of the outcomes of the first two Task Forces. Its intention will be to concentrate on novel high accuracy start-up measurements from Russian VVERs, and potentially also the challenges presented by new, or newly restarted, fuel transient testing reactors.

Figure 1 .

 1 Figure 1. Organization and dependencies between activities of the two Task Forces.

3. 1 . 5

 15 Task 5: Needs, Options, Recommendations and Mechanisms for Conducting Experiments Specifically for Validation of Multi-Physics Modeling and Simulation Tools

Figure 2 .

 2 Figure 2. Example validation pyramid for the important physical phenomena in pellet-clad interaction (PCI) simulation.

† For more information on the LOFT Programme, see www.oecd-nea.org/jointproj/loft/


C Hurley

F Martin

L Marchetti

C Blanc

E Andrieu

J Chêne

Role of grain boundaries in the diffusion of deuterium in nickel base alloy 600 studied by thermal desorption mass spectroscopy

d) e) a) b) c)

The system could be defined one of three ways: (i.) system where GBs are 2 H diffusion short circuits, (ii.) system where GBs act as potential trap sites or (iii.) Fick's second law. Systems (i.) and (ii.) require the determination of many temperature dependent parameters that may not readily be available in literaturetherefore for simplicity (iii.) was initial used.

To use Fick's second law to describe the system, the following assumptions were made:

Context and objectives

The role grain boundaries (GBs) play in the diffusion of hydrogen in polycrystalline alloys has long been debated, with some researchers reporting a "short-circuit" [1][2][START_REF] Harris | Metall Trans A[END_REF][START_REF] Chanfreau | [END_REF][START_REF] Brass | [END_REF][6][START_REF] Dieudonne | Defect Diffusion Forum[END_REF], slowing [START_REF] Yao | [END_REF]9] or even mixed effect depending on grain size and orientation [6,10,11] of GBs on the transport of hydrogen across a metal membrane. This study aims to: 1. address the role of GBs in the diffusion of deuterium ( 2 H) in a model alloy 600 with a grain size of tens of micrometers and 2. derive the diffusion coefficient of 2 H in this alloy through coupled experimental analysis and simulation. (ii) A600-pc polycrystalline elimination of residual trap sites (i) A600-sc single crystal

Thermal desorption mass spectroscopy (TDS)

(1) Cathodic charging of sample with deuterium ( 2 H) at one face (opposite face and sides are protected with a resistant varnish) (2) Transfer from electorchemical cell to TDS system (3) 2 H desorption flux (J ²H ) monitoring in function of t or T a) "traditional" analysislinear temperature ramp b) "hybrid" analysisextended (RT or HT) isothermal period followed by linear temperature ramp Table 2. Experimental conditions for "traditional" and "hybrid" TDS analysis.

RT aging aims to enhance short circuit grain boundary diffusion [START_REF] Brass | [END_REF] while at HT grain boundary dominated diffusion should be minimized [START_REF] Dieudonne | Defect Diffusion Forum[END_REF] 𝜕𝑐 𝜕𝑡 = 𝐷 𝜕 2 𝑐 𝜕𝑥 2

Choice of model and initial hypotheses

• clear similarities in peak temperature, peak intensity and peak form between normalized experimental polycrystalline and single crystal spectra were observed, • diffusion coefficient for 2 H in polycrystalline A600 was derived and validated by fitting an experimental spectra [14,15] and • single crystal A600 TDS spectra were successfully simulated and well reproduced using the derived D from A600-pc.

f) g) h)

1. Acquisition of experimental TDS spectra for A600-pc and A600-sc a) traditional b) hybrid -RT c) hybrid -HT similar intensities, peak temperatures and spectral forms 4. Simulate A600-sc TDS spectra using the 2 H diffusion coefficient derived from A600-pc f) traditional g) hybrid -RT h) hybrid -HT very good reproduction of normalized A600-sc spectra using A600-pc 2 H coefficient

• continue this work using other industrial alloys (i.e. stainless steels) and grain sizes • evaluate using other charging (T and t), aging (T and t) and TDS (φ) conditions • improving numerical model and fitting routine As in this alloy grain boundary interactions can be neglected more focus can be placed on studying the interactions between hydrogen and trap sites such as chromium carbides  this work will be presented in the oral presentation on Tuesday at 10:00 pm.

  a) the diffusion coefficient (D) follows an Arrhenius relationship, b) 2 H diffusion along GBs can be neglected compared to the 2 H flux from interstitial diffusion in the grain volume, c) the sample is free of 2 H before cathodic charging d) constant imposed surface concentration during cathodic charging Role of grain boundaries in the diffusion of deuterium in nickel bas alloy 600 studied by thermal desorption mass spectroscopy

Perspectives 2 . 3 .

 23 Derive 2 H diffusion coefficient in polycrystalline A600 (d) Validate derived diffusion coefficient by simulating aged spectra (e) superior reproduction of normalized aged spectrum using derived diffusion coefficient as compared to literature diffusion coefficients [12,13] 𝑫 𝐜𝐦 𝟐 • 𝐬 -𝟏 = 𝑫 𝟎 • 𝐞𝐱𝐩 -𝑬 𝑫 𝑹𝑻 𝑫 𝟎 = 1.0 ± 0.5 • 10 -2 cm 2 s -1 𝑬 𝑫 = 45 ± 4 kJ. mol -1 Direction de l'Energie Nucléaire Direction déléguée aux Activités Nucléaires de Saclay Département de Physico-Chimie grain boundaries do not play a significant role in the transport of 2 H in alloy 600 for the grain size and experimental conditions studied.

  1,2, † , Frantz Martin 1 , Loïc Marchetti 1,3 , Christine Blanc 2 , Eric Andrieu 2 , Jacques Chêne 4 1 DEN-Service de la Corrosion et du Comportement des Matériaux dans leur Environnement (SCCME), CEA, Université Paris-Saclay, F-91191, Gif-sur-Yvette, France 2 Université de Toulouse, CIRIMAT, UPS/INPT/CNRS,ENSIACET, 4 allée Emile Monso, 31030, Toulouse Cedex 4 3 DEN, MAR, DTCD, SECM, LCLT, CEA Marcoule, 30207 Bagnols sur Cèze 4 CNRS/CEA UMR 8587, CEA Saclay, 91191 Gif-sur-Yvette, France † VTT Technical Research Centre of Finland Ltd. Nuclear Reactor Materials, P.O. Box 1000, 02044 VTT, Finland

	Materials and method
		C	Mn	Si	S	P	Ni	Cr	Cu	Co	Fe	Ti	Al
	A600-pc	0.06	0.82	0.31 <0.001 0.008	base	15.8	0.01	0.01	9.6	0.196 0.164
	A600-sc	n.s.	0.01	0.08	0.01	0.01	base	16.99	0.02	0.01	5.57	<0.01	-
	Conclusions

Table 1 .

 1 Chemical composition of the polycrystalline A600 and A600-like single crystal used in this study

	"Model material" fabrication
	polycrystalline
	industrial A600

(banded grain structure: 55±7 µm & 120±5 µm)

Therefore it can be concluded that... The authors would like to thank EDF for their financial support and M.C. Lafont (CIRIMAT, ENSIACET) her help with material observations.


J Rebufa

E Le Guyadec

D Mazuyer

F Thouverez

Effects of textured journal bearings on vibrations of lightweight rotors

Keywords: Textured bearings, Multi-scale homogenization, Non-linear rotor dynamics 1

A dynamic model of a rotating shaft on two textured hydrodynamic bearings is presented. The hydrodynamic mean pressure is computed using multi-scale periodic homogenization and is projected on a flexible shaft with internal damping. The steady state solutions of an unbalance response are analysed with two sinusoidal texturing patterns. The stability zone and the amplitude of limit cycle are presented. The harmonic balance method applied to the monolithic resolution of the fluid structure coupling is efficient and allows parametric studies of stability.

INTRODUCTION

Advanced rotating machineries, using small rotors or working at high rotating speeds, require a complete control of the vibration components. Lubrication systems, such as hydrodynamic journal bearings play a crucial role in the dynamical behaviour of the whole machinery and have a significant influence on structures vibrations. Predictions of the shaft vibrations are complex due to the nonlinear behaviour of the hydrodynamic bearing. Dynamic coefficients of bearings ( [START_REF] Lund | Review of the Concept of Dynamic Coefficients for Fluid Film Journal Bearings[END_REF]) are widely used for computing the stability of equilibrium points. Nevertheless, they do not take into account non-linear effects of the bearings, that may lead to unstable motion ( [START_REF] Myers | Bifurcation Theory Applied to Oil Whirl in Plain Cylindrical Journal Bearings[END_REF], [START_REF] Vorst | Steady-state behaviour of flexible rotordynamic systems with oil journal bearings[END_REF], [START_REF] Zheng | Nonlinear Dynamic Behaviors of a Complex Rotor-Bearing System[END_REF]). Furthermore, cavitation in the lubricant appears and changes drastically the pressure distribution. Several cavitation models could be retained leading in most cases to different behaviours. The recent adaptation of the Reynolds cavitation model ( [START_REF] Zheng | Nonlinear Dynamic Behaviors of a Complex Rotor-Bearing System[END_REF]) and the Jacobson-Floberg-Olsson ( [START_REF] Giacopini | A Mass-Conserving Complementarity Formulation to Study Lubricant Films in the Presence of Cavitation[END_REF]) to a linear complementary problem is straightforward and efficient. Besides, recent improvements in surface texturing have shown several enhancements in hydrodynamic lubrication such as friction reduction or improvement of drag force ( [START_REF] Etsion | State of the Art in Laser Surface Texturing[END_REF]). One of the first study suggesting that micro irregularities could improve the lubrication has been made almost 50 years ago by Hamilton ([7]). Several studies have been following in order to model the surface state influence with its statistical properties ( [START_REF] Tzeng | Surface Roughness Effect on Slider Bearing Lubrication[END_REF], [START_REF] Patir | An Average Flow Model for Determining Effects of Three-Dimensional Roughness on Partial Hydrodynamic Lubrication[END_REF]). Among them Patir and Cheng introduced the idea of averaging the Reynolds equation by introducing the flow factors method. Elrod has been one of the first to apply the perturbation method with multiple-scale analysis to Reynolds equation ( [START_REF] Elrod | A General Theory for Laminar Lubrication With Reynolds Roughness[END_REF]). The homogenization of Reynolds equation has been mathematically formalized by Bayada et al. with mass conservation model of cavitation ( [START_REF] Bayada | A Double Scale Analysis Approach of the Reynolds Roughness Comments and Application to the Journal Bearing[END_REF], [START_REF] Bayada | An Average Flow Model of the Reynolds Roughness Including a Mass-Flow Preserving Cavitation Model[END_REF]) and has also been developed recently in a transient form and a compressible form by Almqvist et al. ([13], [START_REF] Almqvist | Homogenization of a Reynolds equation describing compressible flow[END_REF]). In these models the assumptions of Reynolds equation need to be verified at the two different scales. Many works have been devoted to analyse the validity of the Reynolds equation in presence of surface texturing as reported De Kraker et al. in [START_REF] De Kraker | Development of a texture averaged Reynolds equation[END_REF]. In their article they adapted the flow factor method with a local CFD computation. The improvement of the computing power also allowed several deterministic static study of the bearing ( [START_REF] Tala-Ighil | Effect of textured area on the performances of a hydrodynamic journal bearing[END_REF], [START_REF] Brizmer | A Laser Surface Textured Journal Bearing[END_REF]). However, to consider the surface effects in a numerical analysis a very fine mesh is needed leading to prohibitive computational time in a dynamic analysis. Homogenization techniques provide a rigorous way of averaging the influence of a periodic texturing pattern. One can study the influence of a micro scale texturing pattern on the macro scale with a relatively coarse mesh, and take into account the average pressure distribution and cavitation zone. Eventually, the mean pressure profile can be projected on the structure at each time step to study the dynamics of the shaft lubricated by two textured bearings. Effects of the texturing pattern on the dynamic behaviour of the system are not well known and need to be properly analysed. So far, few works have been devoted on the stability and transient analysis of the journal bearing with roughness effect ( [START_REF] Ramesh | Stability of Rough Journal Bearings Using Nonlinear Transient Method[END_REF], [START_REF] Turaga | Non-Linear Transient Stability Analysis of a Rigid Rotor Supported on Hydrodynamic Journal Bearings with Rough Surfaces[END_REF], [START_REF] Lin | Application of the Hopf bifurcation theory to limit cycle prediction of short journal bearings with isotropic roughness effects[END_REF]). The resulting dynamic system is complex and efficient motion studies must be used to allow parametric analysis of surface texturing influence. The Harmonic Balance Method (HBM) is an efficient tool for periodic solution research. To convert in the frequency domain complex nonlinearities such as hydrodynamic bearing resultant force the Alternate Frequency/Time domain method (AFT) is efficient and widely used in rotor dynamic systems ( [START_REF] Cameron | An Alternating Frequency/Time Domain Method for Calculating the Steady-State Response of Nonlinear Dynamic Systems[END_REF], [START_REF] Sarrouy | Global search of non-linear systems periodic solutions: A rotordynamics application[END_REF], [START_REF] Grolet | Computing multiple periodic solutions of nonlinear vibration problems using the harmonic balance method and Groebner bases[END_REF]). Moreover, the algebraic expression of the motion equation in the frequency domain allows substantial reduction of the problem to nonlinear degrees of freedom ( [START_REF] Nacivet | A dynamic Lagrangian frequencytime method for the vibration of dry-friction-damped systems[END_REF]). In this study we present a homogenized model taking into account a periodic texturing pattern on a rotating flexible shaft. A monolithic system is presented where mass conserving cavitation and surface texturing are taken into account in the finite hydrodynamic bearing model. The HBM/AFT method with interface reduction is used to study the vibration amplitude of the unbalance response and the stability of limit cycles. The particular cases of longitudinal and transverse sinusoidal texturing will be studied.

FLUID FILM FORCE COMPUTATION WITH HOMOGENIZATION TECHNIQUE

The Reynolds equation with mass conservation algorithm is often used in hydrodynamic lubrication in its isoviscous and unstationary form ( [START_REF] Almqvist | Homogenization of the unstationary incompressible Reynolds equation[END_REF], [START_REF] Almqvist | Homogenization of a Reynolds equation describing compressible flow[END_REF]):

𝜕 𝜕𝑥 ( ℎ 3 12 𝜇 𝜕𝑝 𝜕𝑥 ) + 𝜕 𝜕𝑦 ( ℎ 3 12 𝜇 𝜕𝑝 𝜕𝑦 ) = 𝑅𝛺 2 𝜕𝜃ℎ 𝜕𝑥 + 𝜕𝜃ℎ 𝜕𝑡 𝑝r = 0 , 𝑝 ≥ 0 , r ≥ 0
In this equation 𝑝 is the pressure relative to the cavitation pressure, usually different from the ambient pressure, ℎ describes the film height, 𝜇 the lubricant viscosity, 𝑅 the shaft radius, 𝛺 the shaft rotating velocity. Like in the Elrod algorithm ( [START_REF] Elrod | A Cavitation Algorithm[END_REF]) the film is viewed as a mixture and 𝜃 represents the relative mixture density. The relative motion of the shaft is in the 𝑥 direction, and 𝑦 represents the transverse direction to the flow.

In order to take into account a small variation of amplitude, one can decompose the film height with the plain bearing height ℎ 0 and the periodic local height due to the texturing pattern ℎ 𝜉𝜁 ℎ(𝑥, 𝑦, 𝜉, 𝜁) = ℎ 0 (𝑥, 𝑦) + ℎ 𝜉𝜁 (𝜉, 𝜁)

Coming from the idea of separating the two different scales, the film height can be expressed

ℎ 𝜖 (𝑥, 𝑦) = ℎ (𝑥, 𝑦, 𝑥 𝜖 , 𝑦 𝜖 ) (2)
The main idea of multiple scales method relies on the asymptotic expansion of the solution

𝑝 = 𝑝 0 + 𝜖𝑝 1 + 𝜖 2 𝑝 2 + ⋯ (3)
With these assumptions the homogenized problem is to find the pressure 𝑝 0 solution to the homogenized problem. The expressions of the coefficients 𝐴 𝑖𝑗 and 𝐵 𝑖 are available in the general case in several articles ( [START_REF] Bayada | An Average Flow Model of the Reynolds Roughness Including a Mass-Flow Preserving Cavitation Model[END_REF], [START_REF] Almqvist | Homogenization of the unstationary incompressible Reynolds equation[END_REF]). In this paper we make the assumptions that the film height can be decomposed in the following form:

ℎ(𝑥, 𝑦 , 𝑥𝜖, 𝑦𝜖) = ℎ 1 (𝜃, 𝑦, 𝜃𝜖) ℎ 2 (𝜃, 𝑦, 𝑦𝜖) (4)
With these assumptions the coefficient expressions are available in [START_REF] Bayada | An Average Flow Model of the Reynolds Roughness Including a Mass-Flow Preserving Cavitation Model[END_REF] :

𝐴 11 = ℎ 2 3 h1 -3 ̃ , 𝐴 22 = ℎ 1 3 h2 -3 ̃ , 𝐴 12 = 𝐴 21 = 0 𝐵 1 = ℎ 1 -2 h1 -3 ̃ ℎ 2 ̃ , 𝐵 2 = 0 , Θ = 1 ℎ 2 ̃ℎ1 -2 ̃ ( 𝜃 0 ℎ 2 h1 2 ̃)
Where ⋅ ̃ represents the average operator. The homogenized Reynolds equation in that case is In the case when only one surface is textured it is possible to simplify the problem to a stationary problem ( [START_REF] Almqvist | Homogenization of the unstationary incompressible Reynolds equation[END_REF]). The film height is also independent of which side is textured due to Reynolds hypothesis. In order to compute the mean pressure distribution the finite element method can be used ( [START_REF] Giacopini | A Mass-Conserving Complementarity Formulation to Study Lubricant Films in the Presence of Cavitation[END_REF]). Moreover the problem can be written as a linear complementary problem ([4], [START_REF] Giacopini | A Mass-Conserving Complementarity Formulation to Study Lubricant Films in the Presence of Cavitation[END_REF]). The discretized problem can be written then

𝜕 𝜕𝑥 ( 𝐴 11 12𝜇 𝜕𝑝 0 𝜕𝑥 ) + 𝜕 𝜕𝑦 ( 𝐴 22 12𝜇 𝜕𝑝 0 𝜕𝑦 ) = 𝑅Ω 2 𝜕Θ𝐵 1 𝜕𝑥 + 𝜕Θℎ 0 𝜕𝑡 (5)
{ ℍ(𝐱) 𝐩 ̅ + ℚ 𝐫̅ = 𝐟 f (𝐱) 𝐩 ̅ T 𝐫̅ = 0 , 𝐩 ̅ ≥ 0 , 𝐫̅ ≥ 0 (6)
where 𝐩 ̅ represents the dimensionless discrete pressure vector and 𝐱 the shaft displacements vector.

The dimensionless pressure is obtained with the following formula:

𝑝̅ = 𝑝 12𝜇Ω ( 𝑅 𝐶 ) 2
The boundary conditions in this study are only the atmospheric value of pressure at the bearing edge. As there is no supply groove, the rotating frame is even easier to program for the bearing pressure computation because the fluid grid is moving. Depending of the choice of the cavitation model (Reynolds or JFO), the matrix ℚ can be either the identity matrix ( [START_REF] Zheng | Nonlinear Dynamic Behaviors of a Complex Rotor-Bearing System[END_REF]) or the mass conservation matrix ( [START_REF] Giacopini | A Mass-Conserving Complementarity Formulation to Study Lubricant Films in the Presence of Cavitation[END_REF]).

The matrix ℍ and the vector f f are expressed:

ℍ 𝑖𝑗 (𝐱) = ∫ [𝐴 11 (ℎ) 𝜕𝑊 𝑖 𝜕𝑥 𝜕𝑁 𝑓 𝑗 𝜕𝑥 + 𝑅 2 𝐿 2 𝐴 22 (ℎ) 𝜕𝑊 𝑖 𝜕𝑦 𝜕𝑁 𝑓 𝑗 𝜕𝑦 ] Γ 𝑑Γ 𝐟 f 𝑖 (𝐱) = - 1 2 ∫ 𝑊 𝑖 Γ 𝜕𝐵 1 (ℎ) 𝜕𝑥 𝑑Γ
𝑊 𝑖 and 𝑁 𝑓 𝑗 state for the test function and the interpolation function of the finite element method, respectively. However, the ℍ matrix and 𝐟 f vector are non-linearly dependent of the shaft position. For each new position of the journal the matrix coefficient must be computed with the integration of the flow factors. In order to compute efficiently the matrix ℍ(𝐱) and the vector f f (𝐱) for each time step of the solver method, a polynomial development of the homogenized factor is used. The algorithm allows also the tilt of the shaft in the bearing. For example, one can write the average film height with a misaligned shaft in the pure squeeze rotating coordinates ( [START_REF] Childs | Journal Bearing Impedance Descriptions for Rotordynamic Applications[END_REF]):

ℎ = 1 + 𝑒 1 𝑦 cos(2𝜋𝑥) + 𝑒 2 (1 -𝑦)cos (2𝜋𝑥)
where 𝑒 1 and 𝑒 2 are the eccentricity of the two bearing edges. . From the trinomial expansion we obtain:

(1 + 𝜖 1 𝑦 cos(𝑥) + 𝜖 2 (1 -𝑦)cos (𝑥)) n = ∑ [ n! (n -k)! (k -j)! j! ] 𝑒 1 𝑘-𝑗 𝑦 𝑘-𝑗 𝑒

𝑑Γ

We have then

ℍ 𝑖𝑗 = ∑ [ n! (n -k)! (k -j)! j! ] (𝛼 𝑛 𝑎 (𝑘-𝑗)𝑗 + 𝑅 2 𝐿 2 𝛽 𝑛 𝑏 (𝑘-𝑗)𝑗 ) 0≤𝑗≤𝑘≤𝑛≤3 𝑒 1 𝑘-𝑗 𝑒 2 𝑗
As the 10 matrices 𝑏 𝑝𝑞 and 𝑑 𝑝𝑞 are computed at the program initialization the assembly of the matrix ℍ 𝑖𝑗 is efficient. The same idea is easily applied to assemble 𝐟 f .

Where 𝜎 and 𝜖 are the main component of the stress tensor and the strain tensor in the beam neutral axis direction, respectively. Only the viscous damping is modelled, the hysteretic damping is taken equal to zero. The dynamic equation of a rotor with internal damping at constant rotating velocity is:

𝕄𝐱̈+ [𝕂 𝑏 + Ω𝔾]𝐱̇+ [𝕂 + Ω𝕂 𝑐 ]𝐱 = 𝐟 d (Ω) + 𝐟 s + 𝐟 b (8)
where 𝕄, 𝔾, 𝕂, 𝕂 𝑏 , 𝕂 𝑐 , 𝐟 s , and 𝐟 d , are the mass matrix, gyroscopic matrix, stiffness matrix, structural damping matrix, circulatory matrix, static load force vector, and unbalance force vector, respectively. The gyroscopic matrix 𝔾 arises from the rotating movement of the Timoshenko beam. The circulatory matrix comes from the elastic deformation of the rotating beam with the viscoelastic law ( [START_REF] Zorzi | Finite Element Simulation of Rotor-Bearing Systems With Internal Damping[END_REF]). With the constitutive law (7) the damping matrix is directly related to the stiffness matrix: 𝕂 𝑏 = 𝜂𝕂 . The discretised displacement vector 𝐱 is made of 4 degrees of freedom per node ( [START_REF] Zorzi | Finite Element Simulation of Rotor-Bearing Systems With Internal Damping[END_REF]). In order to project the fluid pressure on the structure one can express the virtual work of the bearings pressure on the shaft:

𝛿𝑊 = ∫ ∫ (cos(𝜃) 𝛿𝑢 + sin(𝜃) 𝛿𝑤) 2 𝜋 0 𝑝 𝑅 𝑑𝜃 𝐿 0 𝑑𝑦 (9)
As the inertia forces are neglected in Reynolds equation the choice of the pressure grid can be made in the rotating frame ( [START_REF] Childs | Journal Bearing Impedance Descriptions for Rotordynamic Applications[END_REF]).

After discretization:

𝛿𝑊 = 𝛿𝐱 T ℙ𝐩 (10)
The matrix ℙ is constructed directly from the fluid interpolation functions and the structure test functions.

Application of the Harmonic Balance Method

The HBM consists in a search of solution of the form of a truncated Fourier series: 𝑥 = 𝑋 0 + ∑ X nc cos(𝑛𝜔𝑡) + 𝑋 𝑛𝑠 sin (𝑛𝜔𝑡) 𝑁 ℎ 𝑛=1 [START_REF] Bayada | An Average Flow Model of the Reynolds Roughness Including a Mass-Flow Preserving Cavitation Model[END_REF] Introducing the following notations

𝑋 = [𝑋 0 X nc X nc … X nc X N h c ] 𝑇 𝑇 = [𝐼 cos(𝜔𝑡) 𝐼 sin(𝜔𝑡) 𝐼 … cos(𝑛𝜔𝑡) 𝐼 sin(𝑛𝜔𝑡) 𝐼 ] 𝑇 ∇ = 𝑑𝑖𝑎𝑔(0 𝑛×𝑛 ∇ 1 … ∇ N h ) with ∇ k = 𝑘 [ 0 𝐼 -𝐼 0 ] (13)
leads to the expression of the displacement, velocity and acceleration with the frequency domain unknowns vector:

𝑥 = 𝑇𝑋 𝑥̇= 𝜔 𝑇𝛻𝑋 𝑥̈= 𝜔𝑇𝛻 2 𝑋 ( 14
)
By replacing these expressions in [START_REF] Bayada | A Double Scale Analysis Approach of the Reynolds Roughness Comments and Application to the Journal Bearing[END_REF] and after a projection on the trigonometric functions basis one obtain a new nonlinear system without any time dependence:

(𝜔 2 𝑁 𝑀 + 𝜔𝑁 𝐶 + 𝑁 𝑘 )𝑋 = 𝐹 + 𝐹 𝑁𝐿 (𝑋) (15)
with 𝑁 𝑀 = 𝑑𝑖𝑎𝑔(𝑀, … , 𝑀) a block diagonal matrix of 2𝑁 ℎ + 1 blocks, 𝑁 𝐶 and 𝑁 𝐾 are constructed identically with the matrices [𝕂 𝑏 + Ω𝔾] and [𝕂 + Ω𝕂 𝑐 ], respectively.

The force F is composed by the static and dynamics load on their respective trigonometric projection.

The 𝐹 𝑁𝐿 represents the non-linear projection in the frequency domain. To compute this vector, the Alternate Frequency/Time domain has been used ( [START_REF] Cameron | An Alternating Frequency/Time Domain Method for Calculating the Steady-State Response of Nonlinear Dynamic Systems[END_REF]). From the frequency unknown vector 𝑋 , the position and velocity are constructed from equation [START_REF] Almqvist | Homogenization of a Reynolds equation describing compressible flow[END_REF]. Then the fluid model allows us to compute the resulting force in the time domain from solving equation ( 8) and projecting the pressure with the projection matrix in [START_REF] Elrod | A General Theory for Laminar Lubrication With Reynolds Roughness[END_REF]. A Discrete Fourier Transform (DFT) of the resulting discrete time vector leads to the contribution of the non-linear force on each harmonic. This method can be resumed by the following scheme:

𝑋 → (𝑥, 𝑥̇) → 𝑓 𝑛𝑙 (Ω, 𝑥, 𝑥̇) → 𝐹 𝑛𝑙 (𝑋)

A basic Newton-Raphson algorithm is used to solve the system [START_REF] De Kraker | Development of a texture averaged Reynolds equation[END_REF]. Once the solution 𝑋 the system ( 15) is found the velocity and displacement in the time domain are reconstructed with equations ( 14). Then, a stability analysis is performed. There are many ways to perform such an analysis with the HBM ( [START_REF] Peletan | A comparison of stability computational methods for periodic solution of nonlinear problems with application to rotordynamics[END_REF]).

In this paper we used the exponentials method to compute the monodromy matrix. The Floquet multipliers, i.e. the eigenvalue of the monodromy matrix are then used to analyse stability of the synchronous motion. An unstable limit cycle is characterized by at one Floquet multiplier with a complex modulus greater than 1 ([29]).

NON-LINEAR RESPONSE TO MASS UNBALANCE

For a given rotating speed one can consider the stability zone of the unbalance response. The parameters of the study are given in table 2. The only static load applied to the horizontal rotor is its weight. In that case a study of a perfectly balanced shaft leads to highly unstable motions. Van de Vorst et al. showed that for lightweight rotors the mass unbalance can achieve stability of the motion ( [START_REF] Vorst | Steady-state behaviour of flexible rotordynamic systems with oil journal bearings[END_REF]). Besides in the case where the cavitation pressure is sub-ambient, the full film lubrication can happen when the pressure is not low enough for the fluid to vaporize in the bearing. The Sommerfeld pressure profile is thus obtained and leads to highly unstable motions. Adding mass unbalance can solve the problem by forcing cavitation. Although the mass unbalance distribution stabilizes the system, it also increases the amplitude of stable vibration. The aim of this paper is to analyse the impact of surface texturing on the minimum mass unbalance to achieve stability. We also focus on the minimum stable vibration amplitude that can be obtained with a certain texturing pattern. Numerical integration shows that the Neimark-Sacker and the flip bifurcation point seem to be subcritical. The minimum mass unbalance to achieve stability for a certain shaft rotating velocity is plotted on the figure 5. For computational time purpose, the Reynolds cavitation hypothesis has been used to solve Reynolds equation in that example. A straightforward dichotomy method has been used to find the Neimark-Sacker bifurcation point. The computation of Floquet multipliers shows that varying the mass unbalance parameter give rise to a generalized Hopf bifurcation or Neimark-Sacker bifurcation (figure 4) in the case of sinusoidal transverse texturing. In the case of sinusoidal longitudinal texturing as well as for plain bearings the parameters of the study give rise to a flip bifurcation. For both systems the minimum unbalance decreases while the rotating speed increases. This means that the system can be unstable during the start-up and stabilise himself after a certain threshold speed. On the diagram we can notice that the minimum unbalance mass is slightly lower for transverse texturing. For example at 2500 rad/s the mass unbalance needed to stabilize the system with sinusoidal transverse texturing (hr = 0.5) is 10% lower than with the sinusoidal longitudinal texturing (hr = 0.5). However, if we pay attention to the vibration amplitude obtained with the boundary stable system we notice that the amplitude value with transverse texturing is 30% higher than with longitudinal texturing. This study is focused on the impact of textured bearing on small rotating shafts vibration. We have presented a powerful method to study the stability of the unbalance response on the system based on multi-scale homogenization of Reynolds equation. The monolithic resolution of the coupled system of equation is performed. The harmonic balance method gives fast and accurate results to study unbalance response and stability of limit cycles.

We noticed that transverse texturing patterns lower the minimum mass unbalance needed for stability although it increases the vibration amplitude. On the opposite the longitudinal texturing pattern allows substantial reduction of the vibration amplitude. However it increases the minimum unbalance needed for stable motion. This study has been performed for one dimensional sinusoidal texturing of finite bearing. However, homogenized flow factor for arbitrary textured bearing could be computed and approximated by polynomials separately and added to the algorithm. Furthermore, the algorithm can be adapted for the averaged flow factors of De Kraker ( [START_REF] De Kraker | Development of a texture averaged Reynolds equation[END_REF]) to take into account recirculation of lubricant in the asperities.

BIBLIOGRAPHY

Figure 1 -

 1 Figure 1 -Homogenized pressure profile with sinusoidal longitudinal pattern with 𝒉 𝒓 = 𝟎. 𝟒 𝑪 𝒓 with the bearing parameter of table 2 with tilted shaft the relative eccentricity 𝒆 𝟏 = 𝟎. 𝟒, 𝒆 𝟐 = 𝟎. 𝟏 A polynomial approximation of 𝐴 11 (ℎ) at the order 3 is taken: 𝐴 11 (ℎ) = ∑ 𝛼 𝑛 ℎ 𝑛 3 𝑛=0 and 𝐴 22 (ℎ) = ∑ 𝛽 𝑛 ℎ 𝑛 3 𝑛=0

Figure 2 - 1 A

 21 Figure 2 -Polynomial interpolation of flow factors for transverse texturing pattern

Figure 3 -

 3 Figure 3 -Rotor -bearing system configurationThis allows taking into account the tilt of the shaft in the bearing. Hence, the dynamic system of equation of the system becomes

Figure 4 -

 4 Figure 4 -Neimark-Sacker and period-doubling bifurcation of limit cycles with shaft rotating speed variation for transverse texturing and period doubling bifurcation for longitudinal texturing pattern 𝒉 𝒓 = 𝟎. 𝟓 𝑪 𝒓

Figure 5 -

 5 Figure 5 -Minimum mass unbalance for stable motion and vibration amplitude of the corresponding boundary stable system 𝒉 𝒓 = 𝟎. 𝟓 𝑪 𝒓 for each rotating speed

Figure 6

 6 shows the vibration amplitude with the limit unbalance for stable motion, for different amplitude of sinusoidal texturing patterns. The shaft rotating speed is 2500 rad/s for all the simulation points. The diagram clearly shows the decrease of longitudinal texturing patterns, while with transverse texturing patterns it significantly increases for ℎ 𝑟 ≥ 0.1 𝐶 𝑟 . Below this value we notice a significant variation of the minimum unbalance. For transverse sinusoidal texturing pattern with ℎ 𝑟 ≈ 0.1 𝐶 𝑟 , both the decrease of the minimum unbalance and the vibration amplitude are achieved.

Figure 6 -

 6 Figure 6 -Vibration amplitude as a function of texturing pattern amplitude for transverse and longitudinal sinusoidal texturing at 2500 rad/s

Table 1 -Homogenized factors for transverse and longitudinal sinusoidal texturing

 1

		Transverse	Longitudinal
	ℎ_𝜖(𝑥, 𝑦, 𝑥𝜖, 𝑦𝜖) 𝐴 11 (𝑥, 𝑦) 𝐴 22 (𝑥, 𝑦) 𝐵 1 (𝑥, 𝑦)	ℎ(𝑥, 𝑦) + ℎ 𝑟 sin(𝑥𝜖) 2 (ℎ(𝑥, 𝑦) 2 -ℎ 𝑟 2 ) 5 2 2 2 ℎ(𝑥, 𝑦) 2 + ℎ 𝑟 ℎ(𝑥, 𝑦) 3 + 3 2 ℎ(𝑥, 𝑦)ℎ 𝑟 2 2ℎ(𝑥, 𝑦)(ℎ(𝑥, 𝑦) 2 -ℎ 𝑟 2 ) 2ℎ(𝑥, 𝑦) 2 + ℎ 𝑟 2	ℎ(𝑥, 𝑦) + ℎ 𝑟 sin(𝑦𝜖) ℎ(𝑥, 𝑦) 3 + 3 2 2 ℎ(𝑥, 𝑦)ℎ 𝑟 (ℎ(𝑥, 𝑦) 2 -ℎ 𝑟 5 2 ) 2 2 2 2 ℎ(𝑥, 𝑦) 2 + ℎ 𝑟 ℎ(𝑥, 𝑦)

Table 2 -Rotor-bearing parameters

 2

	Shaft parameters		Hydrodynamic bearings parameters
	Radius	𝑅 = 2.5 𝑚𝑚	Bearing length	𝐿 𝑝 = 5 𝑚𝑚
	Young Modulus	𝐸 = 211 𝐺𝑃𝑎	Radial clearance	𝐶 𝑟 = 10 𝜇𝑚
	Poisson coefficient	𝜈 = 0.3	Oil viscosity	𝜇 = 0.002 𝑃𝑎. 𝑠
	Density	𝜌 𝑠 = 7800 𝑘𝑔/𝑚 3	Cavitation pressure	𝑝 𝑐𝑎𝑣 = -1 𝑏𝑎𝑟
	Shaft length Internal damping coef.	𝐿 𝑎 = 130 𝑚𝑚 𝜂 = 2.10 -7	Disks parameters Density	𝜌 𝑑 = 2700 𝑘𝑔/𝑚 3
	Coupling radial stiffness	8400 𝑁/𝑚	Radius	𝑅 = 7.5 𝑚𝑚
	Coupling angular stiffness	0.014 𝑁/°	Length	𝐿 𝐷 = 5 𝑚𝑚
	Rotating velocity range	2.500 -3.000 𝑟𝑎𝑑/𝑠 ~23.900 -28.600 𝑟𝑝𝑚


Julia Agullo
email: julia.agullo@cea.fr

Christian Bataillon
email: christian.bataillon@cea.fr

Laurent Trenty
email: laurent.trenty@andra.fr

Corrosion modelling of iron for long term prediction in nuclear waste repository

Keywords: DPCM, nuclear waste, passive iron, CALIPSO, long term prediction

The long-term behavior of nuclear waste canisters in geological repository depends on metallic materials and geochemical features of the surrounding. In clay, H 2 and metallic cations induce geochemical transformations which change pH and redox potential. So the 1000-years behavior of canisters could be predicted by linking geochemical and corrosion models.

The Diffusion Poisson Coupled Model (DPCM) has been proposed as corrosion model for passive iron. To ensure feedback from geochemistry on corrosion rate, backward reactions were introduced. The DPCM is a mechanistic model involving 40 parameters to solve the Poisson and Fick equations in moving boundaries situation. The purpose of this paper is to present the strategy adopted to set up parameters of the CALIPSO code which solves the equation of the DPCM.

It has been shown that one half of these parameters could be evaluated from equilibrium considerations. Then specific experiments could be simulated to set up some parameters. For instance, dissolution rate of magnetite layer monitored by XANES, layer thickness growth monitored by ellipsometry, potentiostatic and Mott-Schottky experiments have been simulated to set up parameters of: dissolution and release processes, inner oxide growth and electrostatic description of outer interface. This strategy was efficient because the DPCM can simulate a large panel of stationary and non-stationary experiments.

Introduction

The long-term behavior of nuclear waste canisters in geological repository depends on metallic materials and geochemical features of the surrounding. In clay, H 2 and metallic cations induce geochemical transformations which change pH and redox potential. So the 1000-years behavior of canisters could be predicted by linking geochemical and corrosion models.

The Diffusion Poisson Coupled Model (DPCM) has been proposed as corrosion model for passive iron. To ensure feedback from geochemistry on corrosion rate, backward reactions were introduced as shown on figure 1. The DPCM is a mechanistic model involving 40 parameters to solve the Poisson and Fick coupled equations in moving boundaries situation. This model is explained in detail in reference [START_REF] Bataillon | Corrosion modelling of iron based alloy in nuclear waste repository[END_REF]. In summary, the DPCM model describes an oxide layer which contains mobile charges (electrons and ionic species). The boundary conditions of the oxide layer are described by kinetic equations in order to solve the Fick equations and by equations describing the state of electric charge of adjacent phases of the oxide layer to the resolution of the Poisson equation. Two particular kinetic steps monitor independently the movement of interfaces of the oxide layer (oxide-metal ≡ inner and oxideenvironment ≡ outer). Thereby, the locations of interfaces are unknown in the DPCM model. The CALIPSO code solves the equations of the Diffusion Poisson Coupled Model. The numerical methods implemented in the CALIPSO code are described in detail in reference [START_REF] Bataillon | Numerical methods for the simulation of a corrosion model with moving oxide layer[END_REF]. In summary, these methods are based on a fully implicit resolution of coupling between Fick, Poisson and moving interfaces equations. This implementation solves the DPCM model in either steady state or dynamic situation. A finite volume scheme with a discretization of Sharfetter & Gummel for convective-diffusive flux is used.

Figure 1: Reaction pathways of the DPCM (Diffusion Poisson Coupled Model) [1]

The purpose of this paper is to present the strategy adopted to set up parameters of the CALIPSO code associated to the DPCM. Some parameters like diffusion coefficients [START_REF] Backhaus-Ricoult | Defects and Cation Diffusion in Magnetite (VII): Diffusion Controlled Formation of Magnetite During Reactions in the Iron-Oxygen System[END_REF] are available in literature. For estimation of others parameters, several ways were considered: thermodynamic equilibrium and some specific experiments dedicated to elementary fluxes like reductive dissolution of magnetite, steady state passive current measurements or oxide layer growth ellipsometry monitoring.

Thermodynamic equilibrium

Description of the system

Considering the equilibrium reaction (1):

Fe 3+ + 3 2 H 2 → ← Fe + 3 H + ( 1
)
The difference of potential, E, between the metal and the solution corresponds to an equilibrium electrochemical potential. This equilibrium could be divided in 2 parallel reaction paths as depicted in figure 2: one for the cations and another for the electrons. In the framework of kinetic any equilibrium can be considered as a steady state where backward and forward paths run as the same rate.

Each path could be considered as an electrochemical equilibrium which are referenced in the literature as standard potential. In the anodic part:

Fe → ← Fe 3+ + 3 e (2)
In the cathodic part:

3 ( H + + e → ← 1 2 H 2 ) (3)
It must be outlined that for electrochemical equilibrium the potential in all phases are homogeneous. Moreover interfaces are always assumed to be uncharged. Hence the voltage drop between the adjacent phases are the so-called potential of zero charge (PZC). Finally the global potential between the solution and metal is the electrochemical potential from the Pourbaix diagram [START_REF] Pourbaix | Atlas of Electrochemical Equilibria in Aqueous Solutions[END_REF]. The potential profile in the system solution/magnetite/metal of figure 2 at the chemical equilibrium ( 1) is presented in figure 3.

Figure 2: Reaction pathways corresponding to the chemical equilibrium (1).

Electronic equilibrium of solution/magnetite/iron system

In the DPCM model, the kinetic of the electronic exchange between the iron and the oxide is described by the following equation:

J e = m e 1 c e 1 -k e 1 c e metal ( 4
)
where c e 1 and c e metal are the electron concentrations in the oxide at the inner interface and in metal, respectively. m e 1 and k e 1 are the kinetic constants of Richardson.

m e 1 = √ k B T 2 π m * and k e 1 = √ k B T 2 π m ** (5)
where m * et m ** are effective masses of electron in the oxide and in metal respectively and k B Boltzmann constant.

The equilibrium corresponding at J e = 0. The solution of this equation gives the concentration of electrons in equilibrium in the oxide at the inner interface. After dimensionless this concentration by N 1 eq = c e 1(eq) Ω ox where Ω ox is molar volume of the oxide and from the integration developed in [START_REF] Bataillon | Corrosion modelling of iron based alloy in nuclear waste repository[END_REF], the following relation is obtained:

N 1 eq = N metal ln [ 1+exp [ -γ ( E -φ 1) ]] (6)
where

N metal = Ω ox k B T n DOS Fe √ m ** m *
. In this expression, n DOS Fe is the average density of state for iron. Knowing (E -φ 1 ) = Δ φ 1 pzc (cf. figure 3), the equation ( 6), it is possible to calculate Δ φ 1 pzc : As (E -φ 1 ) = Δ φ 1 pzc and Δ φ 0 pzc = φ 0 = φ 1 (cf. figure 3), the following relation is obtained:

Δ φ 1 pzc = - 1 γ ln [ exp ( N 1 eq N metal ) -1 ] (7)
[E(pH) -Δ φ 0 pzc ] = Δ φ 1 pzc (8)
The system described on figure 2 also corresponds to an electrochemical equilibrium between iron and magnetite:

3 Fe + 4 H 2 O → ← Fe 3 O 4 + 8 H + + 8 ē (9)
The potential of this equilibrium is:

E = E Fe/ Fe 3 O 4 0 - ln 10 γ pH (10) At 25°C, E Fe / Fe 3 O 4 0 = -0.085 V/NHE (Normal Hydrogen Electrode).
In equation ( 8), the potential E is the electrochemical equilibrium between iron and magnetite in contact of solution (cf. Relation (10)) and depends of pH solution. But Δ φ 1 pzc is a characteristic of ironmagnetite junction which is independent of solution pH. So Δ φ 0 pzc should be a pH function and its pH dependence should be the same as potential E(pH).

In the CALIPSO code, the value of Δ φ 0 pzc is calculated by the following relation:

Δ φ 0 pzc ( pH ) = Δ φ pzc pH =0 - ln 10 γ . n pH pzc0 . pH (11)
The pH dependence of solution is monitored by the value of n pH pzc0 . This dependence should be identical to the one of the Fe/ Fe 3 O 4 / H + couple equilibrium electrochemical potential (cf. relation (10)). Consequently, n pH pzc0 = 1. Introducing (11) and (10) in (8), the following reaction is obtained:

Δ φ pzc pH =0 = E Fe/ Fe 3 O 4 0 -Δ φ 1 pzc with n pH pzc0 = 1 (12)
It should be noted that the description of the local electronic balance at the inner interface and the fact that the potential profile in the magnetite layer is homogeneous at thermodynamic equilibrium leads to the identification of potential Δ φ 1 pzc for inner interface and Δ φ 0 pzc for outer interface. Moreover, the potential Δ φ 0 pzc depends on pH solution and is expressed as an electrochemical potential with reference to the NHE potential. Until now, no data on the kinetic constants ratio values have yet been obtained because Butler-Volmer type reactions have not yet been introduced. The electronic exchange reaction between magnetite and iron does not follow a type of Butler-Volmer law but a law of Richardson. In the following paragraphs, Butler-Volmer type law will be introduced.

The kinetics of the electronic exchange between the oxide and the solution is the kinetic redox of H + /H 2 couple:

J e = m e 0 a H 2 n e H 2 exp ( b e 0 γ φ 0 ) -k e 0 10
-n e pH . pH c e 0 exp ( -a e 0 γ φ 0 )

At the equilibrium, J e = J 0 . After rearranging the following relation is obtained:

k e 0 m e 0 = Ω ox 10 [ n e pH -(a e 0 + b e 0 )] . pH exp [ -(a e 0 + b e 0 ) γ Δ φ 1 pzc ] N 1 eq (Δ φ 1 pzc ) . exp [[ (a e 0 + b e 0 ) -2 n e H 2 ] γ E Fe/ Fe 3 O 4 0 ] (14)
In this relation, the left part is a kinetic constants ratio and is pH independent. So it implies that:

n e pH = (a e 0 +b e 0 ) (15)
The following relation, which is a constant, is finally, obtained:

k e 0 m e 0 = Ω ox exp [-n e pH γ Δ φ 1 pzc ] N metal ln [ 1 + exp ( Δ φ 1 pzc )] exp [[ n e pH -2 n e H 2 ] γ E Fe/ Fe 3 O 4 0 ] (16)
The relation ( 16) does not define a unique solution since it depends on the choice of the reactions of orders n e pH and n e H 2 . In the case of a local electrochemical balance it is useful to take as orders reactions the values of the stoichiometric coefficients of the reaction (3). A convenient solution would be:

n e pH = 1
and n e H 2 = 1 2 then a e 0 + b e 0 = 1 and

k e 0 m e 0 = Ω ox exp [-γΔ φ 1 pzc ] N metal ln [ 1 + exp (Δφ 1 pzc )] (17)
The ratio of kinetic constants depends only on the potential of zero charge of the inner interface and the concentration dimension. The sum of the Butler -Volmer coefficients is equal to 1. This is an expected result for the redox process of the H + /H 2 couple since reactant and product are species in solution.

The same processing could be applied for the Fe 3+ /Fe couple. And finally the description of the thermodynamic equilibrium with the DPCM leads to the following relation listed in the Table 1.

Table 1: Kinetic constants ratio for interfacial reactions for solution/magnetite/iron system described in figure 2.

Mobile charge Outer interface Inner interface

Cations [START_REF] Schmuki | In Situ X-Ray Absorption Near-Edge Spectroscopic Study of the Cathodic Reduction of Artificial Iron Oxide Passive Films[END_REF] that is the thickness and the potential evolution with the charge density. As the current density was constant (-5µA/cm 2 ), the charge density could be converted in time. In the framework of the DPCM, this experimental set up could be described by the figure 4.

m Fe 0 k Fe 0 = 2 P m -2 exp [ 3γ ( E Fe/ Fe 3+ 0 -Δ φ 1 pzc )] a 1 0 +b 1 0 = 1 m Fe 1 k Fe 1 = ( P m 2 -1 ) exp [ 3 ( a 1 L + b 1 L ) γ Δ φ 1 pzc ] Electrons k e 0 m e 0 = Ω ox exp [ γ ( E H + /H 2 0 -Δ φ 1 pzc )]
In this figure, the ferric release reduction can be discarded since the potential is very cathodic.

As the film has been deposited by sputtering, oxygen vacancies are likely very low and as there is no growth between magnetite and tantalum, therefore the oxygen vacancies flux can be discarded. As consequences the experience is monitoring by 2 independent steps. The Oxide Host Lattice Dissolution which is linked with the thickness decrease and the constant k e 0 which is linked to the potential drop during the reduction. It must be highlighted that the kinetic constant for electronic exchange must be those of the junction between magnetite and tantalum. The n DOS Ta should take into account in the paragraph on the Electronic equilibrium of solution, magnetite iron system. Both interfacial reactions depend on the pH being that of the pH buffer (8.4).

Oxide Host Lattice Growth reaction

The rate of the oxide growth at the inner interface (cf. figure 1) determines the motion rate of this interface. The motion rate of the outer interface was determined previously. When both motion rate are equal steady state is achieved. In the literature [START_REF] Lu | Transient growth and thinning of the barrier oxide layer on iron measured by real-time spectroscopic ellipsometry[END_REF], the variation of the steady state thickness as the applied potential is available (figure 9 herein [START_REF] Lu | Transient growth and thinning of the barrier oxide layer on iron measured by real-time spectroscopic ellipsometry[END_REF]). It has been assumed that the process is irreversible as in the Point Defect Model [START_REF] Lu | Transient growth and thinning of the barrier oxide layer on iron measured by real-time spectroscopic ellipsometry[END_REF]. As the oxide growth rate depends on the inner interfacial potential this is obvious that this rate depends on the applied potential through the resolution of the Poisson equation which defines the potential profile in the system. This could be simulated by the CALIPSO code in searching for values of k ox 1 (cf. figure 1) and the associated Butler-Volmer coefficient ( a ox 1 ) in order to fit the experimental curves available in [START_REF] Lu | Transient growth and thinning of the barrier oxide layer on iron measured by real-time spectroscopic ellipsometry[END_REF]. This is important to notice that this set up is likely the more difficult to achieve since in practice the thickness range could be easily got but the right slope not due to the potential dependence.

Conclusions

In this paper a strategy of the CALIPSO set up organised step by step allowing to estimate parameters values without using multifitting processing which is often tricky to implement as soon as the number parameters increased over 3 or 4. The first step was based on equilibrium states which allowed to reduce the number of independent parameters since kinetic constant ratio are constant. In the second step, the dissolution rate of the oxide at the outer interface is available from the literature. It is the major parameter linked to the pH of depassivation or Flade potential. Knowing the potential and pH ranges were passive film exist in the third step, the variation of the steady state thickness with applied potential has allowed to estimate the parameter value of the oxide growth. From this point, the last step is to set up the value for k Fe 0 (figure 1) which controls the level of the anodic passive current. Until now the parameters set up was based on data and experiments which did not deal with corrosion situation. Then it is possible to check of the parameters values set in simulating a corrosion experiment using the evolution of the free corrosion potential and the final damage.

Figure 3 :

 3 Figure 3: Potential profile in the solution/magnetite/metal system of figure 2 at the chemical equilibrium (1)

XANES study of the cathodic reduction of artificial iron oxide passive films [ 5 ]

 5 Schmuki et al. have studied the galvanostatic reduction of Fe 2 O 3 and Fe 3 O 4 sputter deposited thin films by XANES (X-ray Absorption Near Edge Spectroscopy) [5]. Two informations are available in Fig 8 herein

Figure 4 :

 4 Figure 4: Reaction pathways of the DPCM for tantalum covered magnetite layer.

Acknowledgement

The authors would like to thank Andra for financial support.