Ahmed Zouari
email: ahmed.zouari@mines-paristech.fr

During a Reactivity Initiated Accident (RIA), nuclear fuel cladding experiences a multiaxial loading state in which the Pellet-Cladding Mechanical Interaction (PCMI) produces a strain biaxiality ratio (εzz/εθθ) of between 0 and 1. This study examines the effect of strain state and loading path on the fracture of Zircaloy-4 fuel cladding. A new mechanical test has been developed in order to apply several levels of strain biaxiality. Digital Image Correlation (DIC) is used to measure the strain field, to identify the onset of cladding failure, and to obtain the corresponding critical loading state. Results show that the strain biaxiality has a significant effect on the hoop strain at failure, and the smallest fracture strain is obtained for nearly plane strain conditions. A model is proposed for the accumulation of damage due to plastic deformation, based on the Lode parameter of the plastic strain rate tensor. The model is used to predict the failure hoop strain as a function of the strain biaxiality ratio (εzz/εθθ), and good agreement is found between the experimentally measured and predicted failure strains.

Introduction

A Reactivity Initiated Accident (RIA) caused by the ejection of a control rod is a design basis accident that can potentially occur in a Pressurized Water Reactor (PWR). During the first phase of a RIA, the zirconium alloy fuel cladding tubes may be subjected to Pellet-Cladding Mechanical Interaction (PCMI) due to the thermal expansion of the pellet. In this scenario, the expanding pellet imposes a hoop strain and a tensile axial strain on the fuel cladding. Data from previous experiments in research reactor programs simulating the PCMI phase of RIAs show that, at high burnup of the fuel, the cladding undergoes a multiaxial loading state characterized by a strain biaxiality ratio (εzz/εθθ) between plane-strain (no axial strain in the cladding) and equi-biaxial tension (equal tensile strains in the hoop and axial directions). These values were obtained from measurements of the residual axial and hoop strains on cladding tubes [START_REF] Papin | Summary and Interpretation of the CABRI REP-Na Program[END_REF].

Additional experiments are necessary to characterize the failure conditions of fuel cladding in accidental scenarios. However, due to the high cost and complexity of tests in experimental reactors, it would be preferable to obtain the necessary information from laboratory tests that allow the effects of various parameters to be studied individually. To understand the effect of the strain biaxiality ratio (εzz/εθθ) on the hoop strain at fracture of the cladding, laboratory tests should accurately reproduce reactor conditions that affect this parameter. The objective of this study is to develop a laboratory test that reproduces the strain biaxiality imposed on the cladding in the reactor, and to analyze its effect on the hoop strain at failure of non-irradiated Cold-Worked Stress-Relieved (CWSR) Zircaloy-4 cladding. Tests using different strain ratios were coupled with a ductile failure model to examine how the failure strain and failure mechanism are affected by the strain biaxiality ratio.

Previous researchers have used various testing methods to evaluate the effect of strain biaxiality on the fracture behavior of zirconium alloys. The effect of the strain biaxiality ratio for sheet metals has been studied by Yunchang and Koss [START_REF] Yunchang | The Influence of Multiaxial States of Stress on the Hydrogen Embrittlement of Zirconium Alloy Sheet[END_REF], Seo et al. [START_REF] Seo | Forming Limit Diagrams of Zircaloy-4 and Zirlo Sheets for Stamping of Spacer Grids of Nuclear Fuel Rods[END_REF] and Kim et al. [START_REF] Kim | Analytical determination of forming limit curve for zirlo and its experimental validation[END_REF]. Yunchang and Koss [START_REF] Yunchang | The Influence of Multiaxial States of Stress on the Hydrogen Embrittlement of Zirconium Alloy Sheet[END_REF] performed biaxial tests that showed a decrease in the ductility of hydrided Zircaloy-2 sheets as the strain biaxiality increases. The ductility decrease is more pronounced for high hydride concentrations. However, Seo et al. [START_REF] Seo | Forming Limit Diagrams of Zircaloy-4 and Zirlo Sheets for Stamping of Spacer Grids of Nuclear Fuel Rods[END_REF] and Kim et al. [START_REF] Kim | Analytical determination of forming limit curve for zirlo and its experimental validation[END_REF] found either no variation or a slight increase of the ductility when the biaxiality increases from 0 to 1 for Zircaloy-4 and Zirlo sheets, respectively. Each of these studies [START_REF] Yunchang | The Influence of Multiaxial States of Stress on the Hydrogen Embrittlement of Zirconium Alloy Sheet[END_REF][START_REF] Seo | Forming Limit Diagrams of Zircaloy-4 and Zirlo Sheets for Stamping of Spacer Grids of Nuclear Fuel Rods[END_REF][START_REF] Kim | Analytical determination of forming limit curve for zirlo and its experimental validation[END_REF] found that the failure strain of sheets is much higher in pure tension (strain biaxiality ≈ -0.5) than in plain strain (strain biaxiality ≈ 0). These studies have revealed that the strain biaxiality ratio is an important parameter with a significant effect on the strain at failure. However, these results were obtained for sheet materials and cannot be directly applied to cladding tubes, which have different crystallographic texture and mechanical properties.

The effect of biaxiality on the failure of cladding specimens has been studied by many authors with different loading modes. Kaplar et al [START_REF] Kaplar | Mechanical Properties of Unirradiated and Irradiated Zr-1% Nb Cladding[END_REF] examined the impact of biaxiality in plane strain and equi-biaxial conditions for Zr-1%Nb alloy. In their tests, fuel cladding tubes were loaded with a combination of internal pressure and axial tensile force to obtain different strain biaxialities. Their study shows that an increase of the strain ratio (εzz/εθθ) from 0 to 1 leads to a 50% reduction of the hoop strain at failure. In these tests, the tube fails by the onset of plastic instabilities in the form of ballooning followed by localized necking [START_REF] Khan | Ballooning Deformation of Zircaloy-4 Fuel Sheath[END_REF]. This type of loading does not reproduce the displacement-controlled strain imposed on the cladding by the expanding fuel pellet, and thus the test does not reproduce PCMI conditions.

Very few mechanical tests can load tubular samples with a high-biaxiality displacementcontrolled load. One of the tests that best represents PCMI involves placing a pellet-shaped cylindrical device inside a sample of fuel cladding. The device is compressed axially, causing it to expand in diameter, thereby creating a hoop strain in the cladding sample (Fig. 1a). Mishima [START_REF] Mishima | An Improved Testing Device for Measuring Tangential Expandability of Beryllium Tube[END_REF] developed a test in which a ductile, pellet-shaped lead cylinder was compressed by two pistons to measure the failure hoop strain of beryllium cladding tubes. Researchers have also used a segmented expanding mandrel with an axial core of ductile material [START_REF] Garlick | Fracture of zircaloy cladding under simulated power ramp conditions[END_REF]. Compression of the ductile core causes the segmented mandrel to expand, thereby imposing a hoop strain on the sample. Expansion Due to Compression (EDC) tests have been the subject of various recent numerical and experimental studies [START_REF] Grigoriev | Experimental Evaluation of Critical Strain Energy Density for Irradiated Cladding under Simulated RIA Condition[END_REF][START_REF] Dufourneaud | Numerical simulation of the expansion due to compression test[END_REF][START_REF] Dufourneaud | Elastic-Plastic Deformation of a Nuclear Fuel Cladding Specimen under the Internal Pressure of a Polymer Pellet[END_REF][START_REF] Ménager | Effect of External Hydride Layer on Fuel-Cladding Deformation during EDC-Testing[END_REF][START_REF] Saux | High Temperature Expansion Due to Compression Test for the Determination of a Cladding Material Failure[END_REF][START_REF] Latourte | Full field measurements used for assessing industrial issues -Two examples[END_REF] due to their ability to characterize the loading and the failure mode encountered during the PCMI stage of RIAs. In the original version of the EDC test, the ends of the tubular sample are free and unconstrained. During this type of free-end EDC test, the sample shrinks axially as it expands in diameter (Fig. 1a). Therefore, the strain biaxiality in the sample, (εzz/εθθ), is not representative of that experienced by fuel cladding during PCMI. In order to produce a more representative biaxial strain state in the sample, Nobrega [START_REF] Nobrega | Improvements in the design and analysis of the segmented expanding mandrel test[END_REF] introduced two axial retainers at the ends of the tube to restrain it from axial shrinkage (Fig. 1b). The end restraining fixture prevents any change in length of the sample, so the strain biaxiality is much more representative than that of a freeend test. Hellouin de Menibus et al [START_REF] Hellouin De Menibus | Fracture of Zircaloy-4 cladding tubes with or without hydride blisters in uniaxial to plane strain conditions with standard and optimized expansion due to compression tests[END_REF] used this type of fixed-end EDC test on Zircaloy-4 cladding at 25°C and found that the fracture strain is half of that obtained in free-end EDC tests.

Shinozaki et al. [START_REF] Shinozaki | Failure Behavior of the Cladding with Outer Surface Pre-crack in Biaxial Stress Test[END_REF][START_REF] Shinozaki | Improved-EDC tests on the Zircaloy-cladding tube with an outer surface pre-crack[END_REF] and Li et al. [START_REF] Li | The effect of hydride morphology on the failure strain of stress-relieved Zircaloy-4 cladding with an outer surface pre-crack under biaxial stress states[END_REF] developed a type of EDC test in which the sample experiences tensile axial strain as well as tensile hoop strain. Therefore, the strain biaxiality is very similar to that experienced by fuel cladding during PCMI. This test uses two actuators to load the sample. The first applies an axial tensile load to the ends of the tube, and the second compresses a media inside the tube to impose a hoop strain. They performed tests on pre-cracked cladding samples with an outer hydride rim, and they measured the strains using strain gauges mounted on the side of the tube opposite the pre-crack. Their research revealed that the rupture strain of a pre-cracked cladding tube decreases as strain biaxiality increases. In their case, the rupture is governed by the propagation of an existing crack, which they analyzed by calculating the 𝐽 integral.

In order to be able to test the rupture strain of irradiated cladding tubes in conditions representative of the PCMI phase of an RIA, a research program is ongoing at the French Institut Tripartite CEA-EDF-Framatome. An experimental device that can apply a displacementcontrolled biaxial load is being developed at the CEA for application in a hot cell. This Expansion Due to Compression with Tension (EDCT) device is used in the current study to determine the effect of the biaxiality on the fracture strain, and the results are coupled with a numerical procedure to analyze the ductile damage. The EDCT device was used in this study to perform tests on CWSR fuel cladding samples. The initial tests were performed on nonirradiated samples at room temperature with a slow strain rate. The current paper describes the Piston Pellet Cladding Piston results and analysis of these tests. Based on the results, future efforts will be directed at adapting the device to allow testing at higher temperatures, greater strain rates, and implementation in a hot cell.

The new EDCT device and the strain measurement technique are described in Section 2. In Section 3, the results of the mechanical tests are presented. The influence of biaxiality on the fracture strain is studied using free-end and fixed-end EDC tests in addition to EDCT tests. Scanning electron microscopy (SEM) observations of the fracture surfaces are also presented. In Section 4, the tests are simulated using the finite element method (FEM). Section 5 presents a fracture initiation model that reproduces the experimental results by post-processing elastoplastic finite element simulations.

Experimental procedures

The material used in this study was non-irradiated CWSR Zircaloy-4 fuel cladding, which is commonly used in PWRs. The tubes have an outer diameter of 9.5 mm, an inner diameter of 8.36 mm, a wall thickness of 0.57 mm, and a length of 27 mm. In order to ensure that fracture occurred in the central part of the sample, some samples were prepared with a wall thickness that was reduced to 0.52 mm in this region. The pellet inside the sample is made of polytetrafluoroethylene (PTFE) and has a length of 14 mm and a diameter of 8.34 mm.

Three types of EDC tests were used to obtain different strain biaxialities, which consequently affects the stress state. The tests were designed to obtain a mean hoop strain rate of approximately 10 -4 s -1 in the area of interest. All tests were conducted up to failure, except the free-end EDC test, which did not rupture. The tests were carried out at room temperature on an electro-mechanical tensile testing machine equipped with a 100 kN load cell and displacement sensors.

Free-end and fixed-end EDC tests

In this study, EDC tests were performed using a compression cage on a tensile testing machine. The free-end EDC test consists of axially compressing, between two pistons, a ductile PTFE pellet inside a sample of cladding tube with free ends, as shown in Fig. 1a. The axial compression causes the pellet to expand in diameter, thereby imposing a hoop strain on the cladding in the area in contact with the pellet, the area of interest. The test reproduces the displacement-controlled loading mode of PCMI, in which the expanding fuel pellet imposes a strain on the fuel cladding. However, at the beginning of the test, the sample shrinks axially as it expands in diameter, such that 𝜀 𝑧𝑧 /𝜀 𝜃𝜃 ≈ -0.6 and 𝜎 𝑧𝑧 /𝜎 𝜃𝜃 ≈ 0. Thus, the biaxiality created in this test is not representative of the PCMI phase.

The fixed-end EDC test (Fig. 1b) is an improved version of the free-end EDC test. The tubular sample passes through the plates of an end restraining fixture, and the ends are fixed in place. The samples used in the current study are similar to those used by Hellouin de Menibus et al. [START_REF] Hellouin De Menibus | Fracture of Zircaloy-4 cladding tubes with or without hydride blisters in uniaxial to plane strain conditions with standard and optimized expansion due to compression tests[END_REF]. Each sample consists of a tube with a threaded cap welded to each end. Prior to welding the caps in place, a ductile pellet and two stainless steel disks are inserted inside the sample. After passing the ends of the sample through the holes in the restraining fixture, a nut is screwed onto each threaded cap. Each threaded cap contains a central hole so that a piston can pass through it in order to compress the pellet. At the beginning of the test, the strain state resembles that of a free-end EDC test. However, as the test progresses, the strain biaxiality approaches plane-strain (εzz/εθθ ≈ 0). Thus, the biaxiality is more representative than that of a free-end test, and it approaches the lower limit of the PCMI phase near the end of the test.

EDC test with Tensile load (EDCT)

In order to reach higher biaxiality levels, a type of EDC test with an axial tensile load on a zirconium alloy sample was developed in this study. Fig. 2 presents a schematic illustration of the test.

Fig. 2. EDC test with a tensile load

During this test, in addition to the compression of the pellet, four deformable tie rods (two on each side), which are attached to the outer plates of the compression cage, apply a tensile load to the ends of the sample. Therefore, the sample is subjected to simultaneous tensile hoop and axial loads during the test. The strain biaxiality, εzz/εθθ, in this test is controlled by the mechanical behavior of the tie rods. These tie rods are made of stainless steel and are designed to deform plastically during the test. The force that they apply to the sample varies throughout the test and is a function of their geometry and the mechanical tensile behavior of the stainless steel, namely the yield strength and the strain hardening characteristics. Finite element modeling was used to identify a suitable material and geometry for the tie rods in order to obtain a strain biaxiality close to 0. Thereafter, the cross section of the ties rods was altered slightly to either increase or decrease the biaxiality.

Strain measurement technique

Strain measurements on the outer surface of the tube were performed with a digital image correlation (DIC) system (Fig. 3). The surfaces of the specimens were spray-painted with a random black speckle pattern on a uniform white background, and the samples were illuminated with two spotlights. The DIC system is equipped with 2 cameras with 4 megapixel resolution, a focal length of 50 mm, and an aperture of 𝑓/22. The specimen diameter corresponds to about 300 pixels. The frequency of image acquisition was set to 2 frames per second, and the correlation algorithm was performed with the VIC 3D® commercial software. A Gaussian filter was applied to reduce noise. A facet size of 27 pixels with a grid spacing of 7 pixels was used. The 3-D DIC technique can calculate out-of-plane and in-plane full-field displacements of the cladding throughout the test.

Fig. 3. Schematic representation of the EDC test setup and the 3D DIC system

Fig. 4 shows the hoop and axial Green-Lagrange strain components in the tube measured by the 3D DIC system. The strain at failure is recorded in the area of interest of the specimen, at point A in Fig. 4, where the axial and hoop strains reach their maximum absolute values. The axial strain gradient is high due to bulging, but the hoop strain gradient is relatively small, which allows the failure hoop strain to be determined reliably. The Green-Lagrange strains calculated by VIC-3D are transformed into logarithmic strains, which will be used throughout the remainder of this paper.

Fig. 4. Measured hoop and axial Green-Lagrange strain fields for an EDCT test

Experimental results

Mechanical tests

Two specimen shapes were used in this study. The first one has the thickness of industrial fuel cladding (0.57 mm). With this configuration, a strain ratio of 0.4 was achieved. Beyond this value, the excessive axial force causes the sample to rupture in tension outside of the area of interest. To reach higher biaxiality values, the thickness of the tube was reduced using a lathe to 0.52 mm (91% of the initial thickness) over the central 14 mm of the sample. This solution leads to rupture in the area of interest at greater biaxiality values. It was verified that for a fixed-end EDC test and an EDCT test with a biaxiality equal to 0.21, machining to this thickness did not affect the failure strain. The various EDCT tests performed with a wall thickness of 0.52 mm are listed in Table 1. Fig. 5 presents the evolution of the hoop strain obtained from free-end, fixed-end, and EDCT tests as a function of strain biaxiality. All of the strain measurements correspond to point A, shown in Fig. 4. The last point of each curve corresponds to the rupture of the specimen, with the exception of the free-end EDC test, which reached a hoop strain of 0.53 without rupturing. The loading path of the free-end EDC test (black curve) is close to pure tension. Due to plastic anisotropy, the theoretical biaxiality in pure tension is not equal to -0.5. At the beginning of this test, the strain ratio is close to -0.6, and it increases to -0.3 at a hoop strain of 0.53. For the fixed-end EDC test (green curve), the strain ratio is initially -0.3, and it increases to -0.07 at a hoop strain of 0.17. Although the ends of this specimen were fixed, an axial shrinking was measured in the area of interest, and a positive axial strain was measured close to the ends of the sample.

In order to characterize the material ductility at low biaxiality, several EDCT tests (red curves) were performed at strain ratios less than zero. During these tests, the tie rods had a small cross section and allowed the specimen to shrink in the axial direction, but not as much as in the free-end EDC test. Using this design, rupture was observed at a hoop strain of 0.316 for a biaxiality of -0.20. It was not possible to obtain a lower biaxiality in this study, because the tie rods would need to be thinner than could be machined with a uniform diameter. The minimum observed failure strain was 0.095 for a biaxiality of 0.04, which is close to plane strain. As the biaxiality increased from 0.04 to 0.75, the failure strain increased from 0.095 to 0.137.

Fig. 5. Evolution of strain ratio for free-end, fixed-end, and EDCT tests

These same data are presented in a different manner in Fig. 6, in which the hoop strain is plotted versus the axial strain. Using this format, the data can be compared to Forming Limit Diagrams (FLD) that have been reported in the literature. The figure exhibits a characteristic V-shape with a strongly decreasing ductility for biaxialities less than 0 and a slightly increasing ductility for biaxialities between 0 and 0.75. This type of curve was found by [START_REF] Jain | Fracture limit prediction using ductile fracture criteria for forming of an automotive aluminum sheet[END_REF] for AA6111-T4 aluminum sheets, and by [START_REF] Graf | Effect of Changing Strain Paths on Forming Limit Diagrams of AI 2008-T4[END_REF] for aluminum 2008-T4 sheets using the hemispherical punch stretching method. FLDs were obtained for aluminum tubes in [START_REF] He | A modified Marciniak-Kuczynski model for determining the forming limit of thin-walled tube extruded with initial eccentricity[END_REF] and [START_REF] Magrinho | On the determination of forming limits in thin-walled tubes[END_REF]. In [START_REF] He | A modified Marciniak-Kuczynski model for determining the forming limit of thin-walled tube extruded with initial eccentricity[END_REF], the tubes were loaded using a combination of internal pressure and axial force to vary the biaxiality. The authors observed a marked increase in ductility as the biaxiality increased from 0.1 to 1. In [START_REF] Magrinho | On the determination of forming limits in thin-walled tubes[END_REF], EDC tests were performed using rubber as the compression medium. The biaxiality was varied using a clamping system to decrease the effective free length of the specimen, which leads to an increase in biaxiality. The FLD at necking obtained by the authors exhibits the characteristic V-shape. The authors also measured the local strain at fracture by measuring the thickness of the specimen at the crack. It is interesting to note that this local strain at fracture decreases as the biaxiality increases up to 0.35, but this trend should not be directly compared to the macroscopic behavior represented in Fig. 6. In the case of zirconium alloys, Yunchang and Koss [START_REF] Yunchang | The Influence of Multiaxial States of Stress on the Hydrogen Embrittlement of Zirconium Alloy Sheet[END_REF] measured a decrease in the principal strain at failure when biaxiality increases from plain strain to equi-biaxial strain for hydrided Zircaloy-2 sheets using hemispherical punch tests, and they noted that void nucleation is strongly influenced by the presence of the hydrides. FLDs obtained from punch tests on Zirlo and Zircaloy-4 sheets [3-4] exhibit a V-shape similar to that of Fig. 6. Thus, the dependence of the failure strain on the biaxiality observed in this study is similar to trends exhibited by FLDs reported by several previous researchers.

Fig. 6. Major (hoop) strain as a function of the minor (axial) strain.

Observation of damage

In the EDCT tests, failure occurs by axial propagation of a through-thickness crack, whose axial length at the end of the test is equal to the height of the media. The though-thickness crack initiates in the middle of the specimen, where the hoop strain is greatest, and propagates through the thickness of the sample at an angle close to 45° to the principal direction of loading (the hoop direction). This failure scenario closely resembles the macroscopic failure mechanism observed during RIA experimental reactor tests [START_REF] Nea | Nuclear Fuel Behaviour Under Reactivity-initiated Accident (RIA) Conditions[END_REF][START_REF] Fuketa | NSRR/RIA experiments with highburnup PWR fuels[END_REF]. The fracture surfaces of three of the EDCT tests were observed with a Scanning Electron Microscope (SEM). The samples were subjected to biaxialities, εzz/εθθ, of 0.75 (Fig. 7a), 0.04 (Fig. 7b) and 0.20 (Fig. 7c). The fracture surfaces consist of very fine dimples (between 0.5 and 3 μm) that initiated at Zr(Fe,Cr)2 strengthening particles (Laves phase) [START_REF] Saux | Behavior and failure of uniformly hydrided Zircaloy-4 fuel claddings between 25 °C and 480 °C under various stress states, including RIA loading conditions[END_REF], some of which are circled in Fig. 7d. The dimples have approximately the same size for all biaxialities. However, the dimples for tests in which the absolute value of the biaxiality is large (-0.20 and 0.75) appear deeper, whereas they appear to be sheared for a biaxiality close to 0. It therefore can be concluded that failure of Zircaloy-4 cladding is essentially controlled by void nucleation at fine strengthening particles followed by rapid coalescence by internal necking at high biaxiality [START_REF] Thomason | A three-dimensional model for ductile fracture by the growth and coalescence of microvoids[END_REF] and possible coalescence by void shearing at low biaxiality [START_REF] Torki | A mechanism of failure in shear bands[END_REF].

Plastic behavior and modeling of the tests 4.1 Finite Element Analysis

Finite element simulations were carried out to provide insight into the deformation and the failure modes observed during the tests. The simulations were performed using the Cast3M software, which was developed by the CEA [29]. Quadrangular quadratic finite elements with 8 nodes and 4 Gauss points are used for the calculations, which are based on a large displacement and large rotation formulation. The geometry used for the finite element simulations of the EDCT tests is shown in Fig. 8, along with the boundary conditions. In order to simplify the calculations, only a half of the pellet and the cladding was modeled, as was done by previous researchers [START_REF] Grigoriev | Experimental Evaluation of Critical Strain Energy Density for Irradiated Cladding under Simulated RIA Condition[END_REF][START_REF] Dufourneaud | Numerical simulation of the expansion due to compression test[END_REF][START_REF] Dufourneaud | Elastic-Plastic Deformation of a Nuclear Fuel Cladding Specimen under the Internal Pressure of a Polymer Pellet[END_REF]. The left edge of the model is the axis of symmetry and is constrained in the radial (horizontal) direction. The bottom edge of the model is a plane of symmetry at the mid-plane of the sample and is constrained in the axial (vertical) direction. The tie rods are modeled as a single cylindrical structure with an appropriate cross-sectional area. The bottom surface of the tie rod is fixed to the upper surface of the fuel cladding. During the simulation, an axial (vertical) displacement is applied to the upper surface of the tie rod. A displacement of equal magnitude but opposite direction is applied to the upper surface of the piston that compresses the pellet. Both the tie rod and the fuel cladding experience axial tensile stress as a result of this load, and at the same time, the pellet expands as it is compressed to apply a tensile hoop strain to the cladding.

The simulations of the fixed-end and free-end EDC tests use this same geometry, but without the tie rods. For the free-end EDC simulations, the top surface of the cladding tube is free and unconstrained. For the fixed-end EDC simulations, the top surface of the cladding tube is constrained in the axial (vertical) direction.

Fig. 8. Geometry and boundary conditions used for the finite element simulations of the EDCT tests

During the tests, frictional contact occurs between the piston and the pellet, and between the pellet and the cladding. The friction is modeled by Coulomb's friction law in the Cast3m software. It is not possible to experimentally measure the coefficients of friction between the different components. Therefore, simulations were used to identify suitable values.

A series of simulations revealed that the coefficient of friction between the pistons and the pellet has negligible influence on the simulated force-displacement response and on the local deformation of the specimen. Therefore, to facilitate the simulations, a coefficient of 0 was used for this interface. However, the coefficient of friction between the pellet and the cladding has a significant effect on the simulated strain field, strain biaxiality, and forcedisplacement response. Values of 0.1, 0.2, and 0.4 were used in the simulations for this interface. A coefficient of friction of 0.4 leads to the closest agreement between the simulations and the experimental results for all types of EDC tests considered. Therefore, this value was used for all of the simulations in this study. It is interesting to note that Desquines et al [START_REF] Desquines | The issue of stress state during mechanical tests to assess cladding performance during a reactivity-initiated accident (RIA)[END_REF] also used a value of 0.4 for the coefficient of friction between the pellet and the cladding in their simulations.

Material behavior

To simplify the model describing the behavior of the polymer pellet, the effect of the strain rate is not accounted for. Because polymers exhibit a strong strain rate dependence, the model parameters given below should only be used for the strain rate used in this study (10 -4 s -1 ). In the simulations, the flow stress of the polymer pellet, 𝜎 ̅ pel (based on von Mises plasticity), is expressed as:

𝜎 ̅ pel = 𝜎 𝑦 + 𝐾 𝜖̅ 𝑝 𝑛 ( 1
)
𝜖̅ 𝑝 is the equivalent plastic strain, 𝜎 𝑦 is the yield strength equal to 8 MPa, 𝐾 is the strength coefficient equal to 30 MPa, and 𝑛 is the strain hardening exponent equal to 0.5. Young's modulus is 0.6 GPa and Poisson's ratio is 0.47.

A unified anisotropic viscoplastic model for the cladding tube, which was specifically identified in RIA-PCMI conditions by Le Saux et al. [START_REF] Saux | A model to describe the anisotropic viscoplastic mechanical behavior of fresh and irradiated Zircaloy-4 fuel claddings under RIA loading conditions[END_REF], was used in this work. The complete model accounts for strain hardening, temperature, strain rate and irradiation conditions. A simplified version is used in this work. The plastic orthotropy of the material is described by a Hill quadratic yield criterion [START_REF] Hill | The Mathematical Theory of Plasticity[END_REF] with isotropic hardening. The equivalent stress 𝜎 𝐻 is given in the orthotropic radial-circumferential-axial (𝑟, 𝜃, 𝑧) coordinate system of the tube as:

𝜎 𝐻 = [𝐻 𝑟𝑟 (𝜎 𝜃𝜃 -𝜎 𝑧𝑧 ) 2 + 𝐻 𝜃𝜃 (𝜎 𝑧𝑧 -𝜎 𝑟𝑟 ) 2 + 𝐻 𝑧𝑧 (𝜎 𝑟𝑟 -𝜎 𝜃𝜃 ) 2 + 2 𝐻 𝑟𝜃 . 𝜎 𝑟𝜃 2 + 2 𝐻 𝑟𝑧 . 𝜎 𝑟𝑧 2 + 2 𝐻 𝜃𝑧 . 𝜎 𝜃𝑧 2 ] 1/2 (2)
The following coefficients are used following Le Saux et al. [START_REF] Saux | A model to describe the anisotropic viscoplastic mechanical behavior of fresh and irradiated Zircaloy-4 fuel claddings under RIA loading conditions[END_REF] : 𝐻 𝑟𝑟 = 0.58 𝐻 𝜃𝜃 = 0.42 𝐻 𝑧𝑧 = 0.41 𝐻 𝑟𝜃 = 𝐻 𝑟𝑧 = 𝐻 𝜃𝑧 = 1.5

The yield surface is expressed as: 𝜙 = 𝜎 𝐻 -𝜎 ̅ zir = 0 where 𝜎 ̅ zir is the flow stress of the cladding material. The plastic strain rate tensor is calculated by the normality rule and is given by:

𝜀̇𝑝 = 𝑝̇ 𝜕𝜎 𝐻 𝜕𝜎 (3)
𝑝 is the equivalent strain associated with the Hill criterion. 𝜎 ̅ zir is plotted as a function of 𝑝 in Fig. 9. Young's modulus and Poisson's ratio for zirconium are 𝐸 = 95 GPa and 𝜈 = 0.34. Tensile tests were performed to identify the behavior of the 304L steel used for the tie rods, as shown in Fig. 9. Von Mises plasticity is used to describe its behavior. Young's modulus and Poisson's ratio for 304L steel are taken as: 𝐸 = 200 GPa and 𝜈 = 0.3. Finite element simulations were then used to select the cross section and the length of the tie rods to obtain the desired strain biaxiality in the sample during the EDCT tests.

Fig. 9. True stress-strain curves for 304L steel and Zircaloy-4

The simulated hoop strain versus biaxiality curves are compared with the experimental data in Fig. 10. The simulations are in good agreement with the tests, and the different levels of biaxiality are well reproduced. At the beginning of each of the EDCT tests, the biaxiality is very high due to the fact that the tie rods start pulling on the cladding before the media begins to impose a hoop strain. The biaxiality then settles to its final value at a hoop strain of a few percent. The fixed-end EDC test (dark green symbols and curve) is also well reproduced by the simulation. For the cases where the biaxiality is less than 0, the pellet becomes severely deformed and tends to flow between the piston and the cladding. This large distortion leads to convergence difficulties, and the simulations terminate before reaching the rupture point. These simulations are indicated with a star marking the last computed step before the simulation diverged. Each of the simulations can then be post-processed to create a rupture initiation criterion, as discussed in the following section.

Model

Ductile fracture is often described by three distinct phases [START_REF] Garrison | Ductile fracture[END_REF]. In the first phase, cavities nucleate on defects (Laves phase in the case of Zircaloy-4). If the stress triaxiality is high, the cavities can then grow and eventually coalesce to form a macroscopic crack (Fig. 7a and Fig. 7c). Therefore, many micromechanical models of void growth and coalescence [START_REF] Beremin | Cavity formation from inclusions in ductile fracture of A508 steel[END_REF][START_REF] Roy | A model of ductile fracture based on the nucleation and growth of voids[END_REF] consider stress triaxiality as the driving force for damage and rupture. Unfortunately, these models are not very accurate for certain situations in which the stress triaxiality is low, which is the case for the thin cladding tubes used in this study. For example, these models predict that for shear loading, cavities do not increase in size and do not coalesce, which leads to unrealistically high predictions of the ductility. In shear loading, damage and void coalescence may occur via fine localization bands within the material, via changes in the shape and orientation of the cavities that can lead to coalescence [START_REF] Torki | A mechanism of failure in shear bands[END_REF], or via coalescence triggered by void sheeting [START_REF] Haynes | Temperature-dependent void-sheet fracture in Al-Cu-Mg-Ag-Zr[END_REF][START_REF] Tanguy | Plastic and damage behavior of a high strength X100 pipeline steel : experiments and modelling[END_REF][START_REF] Bron | Ductile rupture in thin sheets of two grades of 2024 aluminum alloy[END_REF] or void shearing [START_REF] Rogers | The tensile fracture of ductile metals[END_REF][START_REF] Pineau | Failure of metals I: Brittle and ductile fracture[END_REF][START_REF] Noell | The mechanisms of ductile rupture[END_REF], similar to that seen in Fig. 7b. In fact, researchers have observed that ductility in pure shear can be less than that in pure tension [START_REF] Bai | A new model of metal plasticity and fracture with pressure and Lode dependence[END_REF][START_REF] Bai | A comparative study of three groups of ductile fracture loci in the 3D space[END_REF], despite the lower stress triaxiality. To be able to represent these experimental trends, shear effects can be considered by including the Lode parameter in the fracture model [START_REF] Bai | A new model of metal plasticity and fracture with pressure and Lode dependence[END_REF][START_REF] Dunand | On the predictive capabilities of the shear modified Gurson and the modified Mohr-Coulomb fracture models over a wide range of stress triaxialities and lode angles[END_REF][START_REF] Cao | A Lode-dependent enhanced Lemaitre model for ductile fracture prediction at low stress triaxiality[END_REF].

In this work, a macroscopic fracture initiation model is proposed to describe failure for all levels of biaxiality. The proposed model is "uncoupled" so that the effect of damage on the elasto-plastic behavior is neglected. Uncoupled models for ductile fracture are based on the definition of a damage indicator, 𝐷, which is computed by post-processing elasto-plastic finite element simulations. For each material point (Gauss points in the finite element simulations) 𝐷 is computed for each time step, 𝑡, as:

𝐷(𝑡) = ∫ 𝑊(𝑠𝑡𝑟𝑒𝑠𝑠 𝑠𝑡𝑎𝑡𝑒)𝑑𝑝 𝑝(𝑡) 0 (4)
Fracture is assumed to initiate when 𝐷 reaches 1. The fracture strain 𝑝 𝑓 is the value of the cumulated plastic strain at the fracture initiation point. However, due to the stress state dependence of the damage evolution law, the fracture strain 𝑝 𝑓 does not necessarily correspond to the maximum value of cumulated strain in the entire structure. Note that for a constant stress state, 𝑝 𝑓 = 1/𝑊.

Several forms for the 𝑊 function have been proposed in the literature. For high stress triaxialities, 𝑊 can take the form of the Rice & Tracey model, in which void growth rate is affected only by the triaxiality, 𝑊 = 𝑓(𝑇) = 𝛼 exp (𝛽𝑇) [START_REF] Rice | On the ductile enlargement of voids in triaxial stress fields[END_REF] where 𝛼 and 𝛽 are model parameters, and 𝑇 is the stress triaxiality ratio defined as:

𝑇 = 1 3 𝜎 𝑘𝑘 𝜎 𝑒𝑞 ( 5
)
𝜎 𝑒𝑞 is the von Mises equivalent stress. In the case of plastically anisotropic materials, the use of ''isotropic" quantities such as stress triaxiality or Lode parameter might not be relevant to describe failure [START_REF] Benzarga | Plastic potentials for anisotropic porous solids[END_REF]. In this case, it was recently proposed (based on unit cell simulations) to use a modified triaxiality ratio (𝑇 * ) to describe damage growth [START_REF] Shinohara | Anisotropic ductile failure of a high-strength line pipe steel[END_REF]. It is defined as: An optimization software is used to tune the model parameters (𝑎, 𝑏, 𝑐, 𝑑, 𝑊 0 ) using the Levenberg-Marquardt algorithm [START_REF] Ranganathan | The Levenberg-Marquardt Algorithm[END_REF]. Elasto-plastic simulations for several biaxiality levels were post-processed up to the experimental failure point (hoop strain at failure) as recorded by stereo-correlation. The maximum value of damage over specimen 𝑖, 𝐷 𝑖 max , is computed for this time step. The cost function is then defined as:

𝑇 * = 1 3 𝜎 𝑘𝑘 𝜎 𝐻 (6)
𝛿 = ∑(𝐷 1 𝑚𝑎𝑥 -1)² 𝑁 𝑖=1 (10)
𝐷 max = 1 corresponds to failure of the specimen. The 𝑁 specimens considered in the fit are only those for which 𝛽 > -0.1. The optimized damage parameters are listed in Table 2. Table 2. Parameters of the damage model

Macroscopic failure point

The predicted and experimental macroscopic failure points are compared in Fig. 12. For strain biaxialities between 0 and 0.75, the failure point can be evaluated by post-processing finite element simulations. The agreement for these tests is very good. In particular, strain to failure is minimum for a biaxiality equal to 0 and slowly increases as biaxiality increases. Below 0, it was not possible to carry out finite element simulations up to the experimental failure point, because the simulations encountered convergence problems at a hoop strain of approximately 0.18. In these cases, the failure strain was calculated assuming a fixed strain biaxiality (red symbols). In practice, the simulations were carried out using a single element with a prescribed deformation gradient corresponding to the desired biaxiality. Nonetheless, a very good agreement is found when comparing the experimental and simulated strains to failure. Contour plots of the damage indicator are presented in Fig. 13 for levels of biaxiality at failure equal to 0.04 and 0.75. For both levels of strain biaxiality, damage was localized near the mid-plane of the cladding. For plane strain conditions (biaxiality of 0.04), fracture is predicted to initiate at the external wall of the cladding. At higher biaxiality, it initiates at the inner wall. As failure occurs very rapidly, these trends could not be verified experimentally. This result is the outcome of a complex load history. Regardless of the biaxiality, the maximum cumulated plastic strain is reached at the inner diameter, mainly because the maximum hoop strain is reached at this location. Friction and bending generate a Lode parameter gradient, which causes the initiation point to move from the outer diameter (at low biaxiality) to the inner diameter (at high biaxiality).

Conclusions

During a Reactivity Initiated Accident (RIA), nuclear fuel cladding experiences a multiaxial loading state in which the Pellet-Cladding Mechanical Interaction (PCMI) produces a strain biaxiality ratio εzz/εθθ of between 0 and 1. The effect of loading path on the onset of ductile fracture of non-irradiated Cold-Worked Stress Relieved (CWSR) Zircaloy-4 fuel cladding is studied by submitting tubular specimens to strain states with biaxiality ratios between -0.20 and 0.75. The results show that as the biaxiality ratio increases from -0.20 to 0.04, the hoop strain at fracture decreases sharply from 0.316 to 0.095. However, for biaxiality ratios between 0.04 and 0.75, the hoop strain at fracture increases from 0.095 to 0.137. Examination of fracture surfaces shows that failure is essentially controlled by the coalescence of very small voids that nucleate at strengthening particles. Rupture appears to occur by coalescence via void shearing at low biaxialities, and by coalescence via internal necking at high biaxialities. Simulations using elasto-plastic behaviors of the cladding, pellets and tie rods are used to determine the loading path to fracture for all biaxiality levels. A ductile failure initiation model is proposed, based on the Lode parameter of the plastic strain rate tensor to account for the changes in failure mechanisms and strains to failure for the different levels of biaxiality. The model shows very good agreement between experimental and simulated strains at failure.

Fig. 1 .

 1 EDC test configuration, (a) free-end EDC test, (b) fixed-end EDC test

  EDC test Fixed-end EDC test EDC test with tensile load

Fig. 7 (

 7 Fig. 7 (a) Fracture surfaces observed for EDCT tests with biaxiality 𝜀 𝑧𝑧 /𝜀 𝜃𝜃 = 0.75, 𝜀 𝜃𝜃 (𝑟𝑢𝑝𝑡𝑢𝑟𝑒) = 0.137, (b) with biaxiality 𝜀 𝑧𝑧 /𝜀 𝜃𝜃 = 0.04, 𝜀 𝜃𝜃 (𝑟𝑢𝑝𝑡𝑢𝑟𝑒) = 0.095, (c) with biaxiality 𝜀 𝑧𝑧 /𝜀 𝜃𝜃 = -0.20, 𝜀 𝜃𝜃 (𝑟𝑢𝑝𝑡𝑢𝑟𝑒) = 0.316 and (d) Laves phase precipitates.

Fig. 10 .

 10 Fig. 10. Comparison of the simulated loading paths (solid lines) and experimental data (points) at the mid-plane of the external surface of the sample

Fig. 11 .

 11 Fig. 11. Dependence between triaxiality and Lode parameter for a von Mises material (isotropic) and a Hill material (with the parameters used in this study).

Fig. 12 .

 12 Fig. 12. Numerical predictions and experimental results

  due to convergence problems Failure strain predicted by the model Failure strain predicted for a fixed biaxiality points Experimental data

Fig. 13 .

 13 Fig. 13. Damage distribution computed for the experimental failure for (a) biaxiality = 0.04 and (b) biaxiality = 0.75

  ε zz /ε θθ = 0.04 (b) ε zz /ε θθ = 0.75

Table 1 .

 1 Measured hoop strain at fracture of Zircaloy-4 tubes of thickness 0.52 mm for different strain biaxiality ratios

	A	A

Acknowledgements

The authors would like to thank the project GAINE of the French Institut Tripartite CEA-EDF-Framatome for supporting this study.

𝜎 𝐻 is the Hill equivalent stress defined in eq. ( 2). In the case of moderate stress triaxiality, the 𝑊 damage function can be written as a function of both the stress triaxiality and the Lode parameter using an additive form, 𝑊 = 𝑓(𝑇) + 𝑔(𝐿) as proposed in [START_REF] Defaisse | Ductile fracture of an ultra-high strength steel under low to moderate stress triaxiality[END_REF] to represent failure of an ultra-high strength steel, or using a multiplicative form 𝑊 = 𝑓(𝑇) × 𝑔(𝐿) proposed in [START_REF] Cortese | Prediction of ductile failure in materials for onshore and offshore pipeline applications[END_REF] and [START_REF] Cortese | A J2-J3 approach in plastic and damage description of ductile materials[END_REF] to represent failure of pipeline steels. More complex forms for 𝑊 have been proposed in the literature. Most of them are reviewed in [START_REF] Bai | A new model of metal plasticity and fracture with pressure and Lode dependence[END_REF]. In this work, the following definition of the Lode parameter will be used for any second order tensor 𝑡. Let 𝑡 ′ be the deviator of 𝑡 and 𝑡 1 ′ ≥ 𝑡 2 ′ ≥ 𝑡 3 ′ its eigenvalues. The Lode parameter of tensor 𝑡 is then defined as [START_REF] Defaisse | Ductile fracture of an ultra-high strength steel under low to moderate stress triaxiality[END_REF][START_REF] Vadillo | A modified Gurson model to account for the influence of the Lode parameter at high triaxialities[END_REF]:

Note that the value of Lode(𝑡) is between -1 and 1. To define the 𝑊 function, the Lode parameter of the stress tensor is often used [START_REF] Bai | A new model of metal plasticity and fracture with pressure and Lode dependence[END_REF]. The Lode parameter is, however, used to describe the fact that ductility is minimum under generalized plane strain conditions, which can lead to easy band localization [START_REF] Rice | The localization of plastic deformation[END_REF]. Based on this consideration, it is proposed in this work to express the failure criterion as a function of the Lode parameter of the plastic strain rate tensor 𝜀ṗ. Note that in the case of von Mises plasticity, the Lode parameter of the stress tensor (𝐿 in the following) is equal to that of the plastic strain rate tensor (referred to as 𝐿 * in the following).

In this work, the thickness of the tube is such that a plane stress state is almost achieved. In this case, the Lode parameter and the triaxiality are linked. In the case of von Mises plasticity with a strain biaxiality between -0.5 and 1, this relation is expressed as:

In the case of Hill plasticity, the relation between 𝑇 * and 𝐿 * cannot be explicitly expressed. This relation is plotted in Fig. 11 for the investigated strain biaxiality range and the indicated Hill parameters.

In this study, the damage 𝑊 is expressed as a function of the Lode parameter 𝐿 * as:

𝑎, 𝑏, 𝑐, 𝑑 and 𝑊 0 are parameters to be adjusted. sign(𝐿 * ) is equal to -1 if 𝐿 * < 0 and to +1 if 𝐿 * > 0. Note that the function has an even part (∝ |𝐿 * | 𝑏 ) and an odd part (∝ sign( 𝐿 * )|𝐿 * | 𝑑 ) to distinguish equibiaxial tension (𝐿 = -1) and pure tension (𝐿 = 1).


Loïc Paulevé
email: loic.pauleve@irccyn.ec-nantes.fr

Adrien Richard
email: richard@i3s.unice.fr

Static Analysis of Boolean Networks Based on Interaction Graphs: a Survey ⋆

Boolean networks are discrete dynamical systems extensively used to model biological regulatory networks. The dynamical analysis of these networks suffers from the combinatorial explosion of the state space, which grows exponentially with the number n of components. To face this problem, a classical approach consists in deducing from the interaction graph of the network, which only contains n vertices, some information on the dynamics of the network. In this paper, we present results in this topic, mainly by focusing on the influence of positive and negative feedbacks.

Introduction

Biological regulatory networks model the influence of biological components between each other. While such relations can be deduced through experiments by biologists, such systems can raise complex global dynamical behaviours and motivate the use of formal verification techniques to validate a model, and predict and understand particular behaviours, such as proteins productions dynamics [START_REF] Endy | Modelling cellular behaviour[END_REF][START_REF] De | Modeling and simulation of genetic regulatory systems: A literature review[END_REF].

Boolean Networks (BNs) are a typical formalism used to model biological regulatory networks dynamics (see the seminal work of Stuart Kauffman [START_REF] Kauffman | Metabolic stability and epigenesis in randomly connected nets[END_REF][START_REF] Kauffman | Origins of Order Self-Organization and Selection in Evolution[END_REF] and René Thomas [START_REF] Thomas | Boolean formalization of genetic control circuits[END_REF][START_REF] Thomas | Biological Feedback[END_REF]). A BN associates to each component a Boolean value and a local logical function describing its evolution in front of the value of other components of the network. The global dynamics is then a Boolean map applying the local functions on component values, with different settings of iteration schemes.

Because the size of the state space is exponential with the number of components, classical analysis, such as model-checking techniques, may be intractable with large systems. To cope with this combinatorial explosion, an amount of work establish relationships between the interaction graph of a BN and some properties on the dynamics. An interaction graph references the positive and negative influences between components into a signed directed graph having one vertex per component; and hence is a compact and static abstraction of BN dynamics.

This paper aims at giving an insight of the different results on the relationships between the interaction graph and the dynamical properties of BNs. These kinds of results are of particular interest in the modelling of regulatory networks: those interaction graphs are generally the most reliable information on the biological system, while few knowledge on their precise behaviours are available.

This paper is structured as follows. Sect. 2 presents the main definitions used in this paper. Sect. 3 details the major results extracting dynamical properties of a BN from the topology of its interaction graph. Finally, we discuss these state-of-the-art approaches, and mention other static analysis techniques for BNs in Sect. 4.

Definitions

Notations. Let n be a positive integer, [n] = {1, . . . , n}, and B = {0, 1}. The ith component of a point x ∈ B n is denoted x i . If I ⊆ [n], then x I denotes the point y ∈ B n such that y i = 1x i if i ∈ I, and y i = x i otherwise. x is an abbreviation for x [n] , and if i ∈ [n], we write x i instead of x {i} . For all x, y ∈ B n , we denote by ∆(x, y) the set of i ∈ [n] such that

x i = y i . Let f : B n → B n . The ith component of f (x) is denoted f i (x).
Hence, f i may be seen as a function from B n to B. For all x ∈ B n , we set ∆f (x) = ∆(x, f (x)). Observe that f (x) = x ∆f (x) .

Boolean networks

A Boolean network is a discrete dynamical system that consists in a collection of n Boolean components evolving along a discrete time by mutual interactions. Dynamics of such a system are usually described by a directed graph on B n that is constructed from a Boolean function f : B n → B n . Several iteration graphs can be considered. In the scope of this paper, we focus on the synchronous, asynchronous and generalized iteration graphs.

Definition 2.1 [Iteration graphs]

• The synchronous iteration graph of f , denoted by sig(f ), is the directed graph on B n that contains an arc

x → y if and only if y = f (x) = x. • The asynchronous iteration graph of f , denoted by aig(f ), is the directed graph on B n that contains an arc

x → y if and only if there exists i ∈ ∆f (x) such that y = x i . • The generalized iteration graph of f , denoted by gig(f ), is the directed graph on B n that contains an arc

x → y if and only if there exists a non empty set I ⊆ ∆f (x) such that y = x I .

On the one hand, Kauffman [START_REF] Kauffman | Metabolic stability and epigenesis in randomly connected nets[END_REF][START_REF] Kauffman | Origins of Order Self-Organization and Selection in Evolution[END_REF] proposed to model gene (regulatory) networks with synchronous iteration graphs. One the other hand, Thomas [START_REF] Thomas | Boolean formalization of genetic control circuits[END_REF][START_REF] Thomas | Biological Feedback[END_REF][START_REF] Thomas | Multistationarity, the basis of cell differentiation and memory. ii. logical analysis of regulatory networks in terms of feedback circuits[END_REF] argued that it is unlikely that several components change simultaneously their value, and he proposed to model gene networks with asynchronous iteration graphs. In constrast with the synchronous description, the asynchronous description is indeterministic: a state x may have several successors (the number of successors of x is precisely |∆f (x)|). An other interesting feature of the asynchronous description is that it can been seen has a relevant approximation of continuous descriptions based on piece-wise linear differential systems [START_REF] Glass | Logical analysis of continuous, non-linear biochemical control networks[END_REF][START_REF] Glass | Combinatorial and topological methods in nonlinear chemical kinetics[END_REF][START_REF] Snoussi | Qualitative dynamics of a piecewise-linear differential equations : a discrete mapping approach[END_REF][START_REF] Thomas | Biological Feedback[END_REF][START_REF] Snoussi | Logical identification of all steady states : the concept of feedback loop caracteristic states[END_REF]. Note that the generalized iteration graph contains, as subgraphs, both the synchronous and asynchronous iteration graphs.

In these three cases, the interpretation of f : B n → B n as a gene network is the same: [n] corresponds to a set of genes, and at state x ∈ B n : the protein encoded by the ith gene is present if x i = 1 and absent otherwise, and the ith gene is expressed ("on") if f i (x) = 1 and unexpressed ("off") otherwise. Hence, the gene expressions depend on the "proteins concentration", and the evolution of these concentrations depends on the gene expressions.

Example 2.2 [Running example] n = 3 and f is defined by:

f 1 (x) = x 3 ∧ (x 1 ∨ x 2 ), f 2 (x) = x 3 ∧ x 1 , f 3 (x) = x 1 ∨ x 2 ∨ x 3 .
The synchronous, asynchronous and generalized iteration graphs of f are:

sig(f ) aig(f ) gig(f )

Fixed points and attractors

A point x ∈ B n is a fixed point of f if x = f (x). Fixed points are of particular interest, since they correspond to stable states (in each iteration graph, x is a fixed point of f if and only if it has no successor). In the context of gene networks, fixed points correspond to stable patterns of gene expressions and are often associated to particular biological functions. In the next definition, the notion of fixed point is extended to the one of attractor.

Definition 2.3 [Attractors]

Let Γ be the synchronous, asynchronous or generalized iteration graph of f . The attractors of Γ are the smallest non empty subsets A ⊆ B n without leaving arc, i.e. such that for every arc x → y in Γ, if x ∈ A then y ∈ A. An attractor of cardinality at least two is said cyclic.

It derives that a cyclic attractor does not contain any fixed point. As a consequence, once a the system is inside a cyclic attractor, it cannot reach a fixed point, and it necessarily describes sustained oscillations. Other easy observations on attractors follow. Proposition 2.4 A point x is a fixed point of f if and only if {x} is an attractor of Γ. The non-cyclic attractors of Γ can thus be identified to the fixed points of f . Also, for every x in B n , there exists at least one path starting from x that reaches an attractor. Consequently, Γ has always at least one attractor.

Example 2.5 [Running example continued] The attractors of aig(f ) and gig(f ) are the fixed point 000 and the cyclic attractor {001, 011, 101, 111}. The attractors of sig(f ) are the fixed point 000 and the cyclic attractor {011, 101, 111}.

Interaction graphs

Inherent to the notion of Boolean network is the one of interaction graph. Such a graph depicts, in a very coarse way, the qualitative interactions between components, and is usually represented as a directed graph with [n] as vertex set. An arc from j to i then means that the evolution of the ith component depends on the evolution of the jth one. Here, we consider two kinds of interactions: the positive and negative ones. This hence leads us to consider signed directed graphs.

Definition 2.6 [Signed directed graph] A signed directed graph on [n] is a graph G whose the set of vertices is [n],
and whose the set of arcs is a subset of

[n] × {+, -} × [n]. If (j, s, i
) is an arc of G, we say that G has an arc from j to i of sign s. G is simple if, for every i, j ∈ [n], there exists at most one arc from j to i. A positive (resp. negative) cycle of G is an elementary directed cycle that contains an even (resp. odd) number of negative arcs. The length of a cycle is the number of arcs it involves.

We are now in position to define, in a natural way, the interaction graph of a network whose dynamics are described from f : B n → B n . Definition 2.7 [Interaction graph of f ] The interaction graph of f , denoted by G(f ), is the signed directed graph on [n] defined by: for all i, j ∈ [n], there exists a positive (resp. negative) arc from j to i if and only if there exists x ∈ B n with x j = 0 such that

f i (x) < f i (x j ) (resp. f i (x) > f i (x j )).
In other words, there exists a positive (resp. negative) arc from j to i if, in at least one state x, an increase of x j induces an increase (resp. decrease) of f i (x), i.e. the appearance of the jth protein turns on (resp. off) the ith gene.

Example 2.8 [Running example continued] The interaction graph of f follows. It has a negative cycle of length one, and a positive cycle of length one, (two of) two and three. Note also that it is simple.

Absence of cycle and convergence toward a unique fixed point

The following fundamental fixed point theorem of Robert [START_REF] Robert | Iterations sur des ensembles finis et automates cellulaires contractants[END_REF][START_REF] Robert | Discrete iterations: a metric study[END_REF] shows that the presence of a cycle in the interaction graph is necessary to prevent a simple synchronous convergence toward a unique fixed point. This convergence result was latter proved for asynchronous iteration graphs by Robert [START_REF] Robert | Les systèmes dynamiques discrets[END_REF], and for generalized iteration graphs by Bahi and Michel [START_REF] Bahi | Convergence of discrete asynchronous iterations[END_REF] (actually, these authors proved the convergence for a much more general class of iteration schemes that include delays of communication between components). Theorem 3.2 If G(f ) has no cycle, then f has a unique fixed point, and every path of gig(f ) reaches this fixed point.

Absence of positive cycle and upper bound on the number of fixed points

In 1980, Thomas stated the following general rule: a necessary condition for a dynamical system to admit several stable states is the presence of a positive cycle in its interaction graph [START_REF] Thomas | On the relation between the logical structure of systems and their ability to generate multiple steady states or sustained oscillations[END_REF]. In the context of gene networks, one can interpret this rule as follow: positive cycles are necessary for differentiation processes [START_REF] Thomas | Biological Feedback[END_REF]. This rule was first formally stated and proved in continuous frameworks [START_REF] Plathe | Feedback loops, stability and multistationarity in dynamical systems[END_REF][START_REF] Snoussi | Necessary conditions for multistationarity and stable periodicity[END_REF][START_REF] Gouzé | Positive and negative circuits in dynamical systems[END_REF][START_REF] Cinquin | Positive and negative feedback: strinking a balance between necessary antagonists[END_REF][START_REF] Soulé | Graphical requirements for multistationarity[END_REF][START_REF] Soulé | Mathematical approaches to differentiation and gene regulation[END_REF]. A Boolean version of Thomas' rule have then been stated and proved by Remy, Ruet and Thieffry [START_REF] Remy | Graphic requirements for multistability and attractive cycles in a boolean dynamical framework[END_REF].

Theorem 3.3 If G(f ) has no positive cycle, then f has at most fixed point.
The following theorem, established in [START_REF] Richard | Necessary conditions for multistationarity in discrete dynamical systems[END_REF], shows that positive cycles are more generally necessary for the presence of several asynchronous attractors. This theorem is not valid in the synchronous case. Also, since aig(f ) has always at least one attractor (cf. Proposition 2.4), the theorem remains valid if "at most one" is replaced by "a unique". Finally, using it as a base case, one can obtain an upper bound on the number of asynchronous attractors that only depends on the structure of connections between positive cycles [START_REF] Aracena | Positive and negative circuits in discrete neural networks[END_REF][START_REF] Aracena | Maximum number of fixed points in regulatory boolean networks[END_REF][START_REF] Richard | Positive circuits and maximal number of fixed points in discrete dynamical systems[END_REF]. If G(f ) has no positive cycle, then I = ∅ satisfies the conditions of the statement, and the resulting upper bound is one (hence, Theorem 3.5 generalizes Theorem 3.4). The bound shows that it is not the number of positive cycles which is important to obtain several fixed points (or asynchronous attractors), but much more the structure of the connections between positive cycles. For instance, if G(f ) has a thousand of positive cycles, and if all these cycles share a same vertex i, then I = {i} satisfies the conditions of the statement, and the resulting upper bound is only two. In other words, in order for f to have a lot of fixed points, G(f ) has to contains a lot of "rather disconnected" positive cycles.

We have seen that if G(f ) has no positive cycle, then f has at most one fixed point. The following theorem of Aracena [START_REF] Aracena | Positive and negative circuits in discrete neural networks[END_REF][START_REF] Aracena | Maximum number of fixed points in regulatory boolean networks[END_REF] shows that, under a very weak additional condition, f has actually no fixed point.

Theorem 3.6 If G(f ) has a minimal in-degree at least one and has no positive cycle, then f has no fixed point.

Aracena proved this theorem under the additional assumption that G(f ) is simple, but this condition can be easily removed. Now, observe that this theorem implies quite directly Theorem 3.3. Indeed, suppose that G(f ) has no positive cycle. If G(f ) has a minimal in-degree at least one and has no positive cycle, then by Theorem 3.6, f has no fixed point. Otherwise, there exists a vertex, say n without loss generality, of in-degree zero in G(f ). So

f n = cst = c. Let h : B n-1 → B n-1 defined by h i (x) = f i (x, c) for all x ∈ B n-1 and i ∈ [n -1].
It is easy to see that G(h) is a subgraph of G(f ). So G(h) has no positive cycle, and proceeding by induction on the number of components, one obtains that h has at most one fixed point. Since it is clear that h and f have the same number of fixed points, we deduce that f has at most one fixed point.

Absence of negative cycle and lower bound on the number of fixed points

Thomas also stated a general rule on negative cycles: a necessary condition for a dynamical system to produce sustained oscillations is the presence of a negative cycle in its interaction graph [START_REF] Thomas | On the relation between the logical structure of systems and their ability to generate multiple steady states or sustained oscillations[END_REF]. In the context of gene networks, one can interpret this rule as follow: negative cycles are necessary for homoeostasis phenomena [START_REF] Thomas | Biological Feedback[END_REF]. The next theorem, proved in [START_REF] Richard | Negative circuits and sustained oscillations in asynchronous automata networks[END_REF], may be seen as a Boolean version of second Thomas' rule (and it is not valid in the synchronous case).

Theorem 3.7 If G(f ) has no negative cycle, then aig(f ) has no cyclic attractor. Hence, if G(f ) has no negative cycle, then f has at least one fixed point.

The second assertion is a consequence of the first one, since if G(f ) has no negative cycle, then aig(f ) has no cyclic attractor, and since aig(f ) has always at least one attractor, we deduce that aig(f ) has an attractor of size one, which is a fixed point of f (cf. Proposition 2.4). In [START_REF] Goles | Sequential operator for filtering cycles in boolean networks[END_REF], Goles and Salinas provide an algorithm with running time O(n 2 ) to compute, under the absence of negative cycle in G(f ), a fixed point of f (in the general case, the problem of finding a fixed point is NP-complete).

Note that Theorems 3.3 and 3.7 provides a nice "proof by dichotomy" of (the first assertion in) the theorem of Robert (Theorem 3.1): if G(f ) has no cycle, then following Theorem 3.3 (resp. 3.7), f has at most (resp. at least) one fixed point, so f has indeed a unique fixed point.

While Theorem 3.6 provides a sufficient condition for the absence of fixed point in term of absence of positive cycle, the following theorem, also established by Aracena [START_REF] Aracena | Positive and negative circuits in discrete neural networks[END_REF][START_REF] Aracena | Maximum number of fixed points in regulatory boolean networks[END_REF], provides a sufficient condition for the presence of multiple fixed points in term of absence of negative cycle. Theorem 3.8 Let G(f ) * be the signed directed graph on [n] that contains an arc from j to i of sign s if G(f ) has an arc of sign s from j to i or from i to j. If G(f ) is simple, has a minimal in-degree at least one, and if G(f ) * has no negative cycle, then there exists x ∈ B n such that x and x are fixed points of f . It is easy to see that if G(f ) is strongly connected and has no negative cycle, then G(f ) is simple, has a minimal in-degree at least one, and G(f ) * has no negative cycle. Hence, the next theorem is a consequence of the previous one. Theorem 3.9 If G(f ) is strongly connected and has no negative cycle, then f has at least two fixed points.

The following theorem, established in [START_REF] Paulevé | Topological Fixed Points in Boolean Networks[END_REF], generalizes Theorem 3.8 by characterizing the number of topological fixed points of f , i.e. the number of fixed points of f that only depend on G(f ). Theorem 3.10 Let G(f ) # be the signed directed graph on [n] that contains an arc from j to i of sign s if and only if G(f ) has an arc from j to i of sign s and no arc from j to i of sign -s. Let G(f ) * be as in Theorem 3.8. Let us say that x is a topological fixed point of f if it is a fixed point of every function h : B n → B n such that G(h) = G(f ). We have the following two properties:

• Let p be the number of connected components in G(f ) # . If the following three conditions are satisfied then f has exactly 2 p topological fixed points, and otherwise, f has 0 topological fixed point: (i) the minimal in-degree of G(f ) # is at least one; (ii) G(f ) * has no negative cycle; (iii) for every i ∈ [n], there exists at most one j ∈ [n] such that G(f ) has both a positive and a negative arc from j to i.

• If x is a topological fixed point of f , then x is a topological fixed point of f . Note that if G(f ) is simple, then G(f ) # = G(f )
and the third condition is trivially satisfied. We then recover Theorem 3.8.

Comparison of iteration graphs

As said above, there are two principal iteration graphs used to model BNs dynamics, the synchronous and asynchronous ones. It is then of great interest to try to compare them. There are very few results in this topics. However, Noual [START_REF] Noual | General iteration graphs and boolean automata circuits[END_REF] recently established the following fundamental theorem, which shows that under some conditions on the cycles of G(f ), every path of the generalized iteration graph, and in particular every path of the synchronous iteration graph, can be "imitated" by at least one path of the asynchronous iteration graph. Other results concerning comparison of iteration schemes are provided in [START_REF] Tchente | Cycles generated by sequential iterations[END_REF][START_REF] Robert | Les systèmes dynamiques discrets[END_REF][START_REF] Goles | Comparison between parallel and serial dynamics of boolean networks[END_REF][START_REF] Aracena | On the robustness of update schedules in boolean networks[END_REF] Theorem 3.11 Assume that G(f ) is simple, has no positive cycle of even length, and has no negative cycle of odd length. If gig(f ) has an arc x → y, then aig(f ) has a path from x to y of length |∆(x, y)|. Hence, the number of attractors in gig(f ) and sig(f ) is at least the number of attractors in aig(f ).

Network reduction

Since the number of possible states of a BN increases exponentially with the number n of components, it is interesting to try to reduce the system, typically by removing components, without losing too many dynamical properties. The following theorem, stated with other notations in [START_REF] Naldi | A reduction of logical regulatory graphs preserving essential dynamical properties[END_REF], shows that if the nth component does not interact with itself, it can be suppressed in a natural way, without affecting the number of fixed points, the reachability by asynchronous paths, and the main features of the interaction graph. (For sake of simplicity, this theorem only deals with the suppression of the nth component, but it can be very easily adapted to the suppression of any component.) Theorem 3.12 Let f : B n-1 → B n-1 be defined from f : B n → B n by:

∀x ∈ B n-1 , ∀i ∈ [n -1], fi (x) = f i (x), x = (x, f n (x, 0)) ∈ B n .
If G(f ) has no arc from n to itself, then the following three properties hold:

• A point x is a fixed point of f if and only if x is a fixed point of f . As a consequence, f and f have the same number of fixed points.

• aig( f ) has a path from x to y if and only if aig(f ) has a path from x to ỹ. As a consequence, the number of attractors in aig( f ) is at least the number of attractors in aig(f ).

• If G( f ) has a path from j to i with an even (resp. odd) number of negative arcs, then G(f ) has a path from j to i with an even (resp. odd) number of negative arcs. As a consequence, if G(f ) has no positive (resp. negative) cycle, then G( f ) has no positive (resp. negative) cycle.

Note that this theorem gives another way to prove the (first assertion in the) theorem of Robert (Theorem 3.1). Indeed, suppose that G(f ) has no cycle. Then it has no arc from n to itself, so f and f have the same number of fixed points, and G( f ) has no cycle (according to the third point). Using an induction on the number of components, we deduce that f has a unique fixed point, so that f has a unique fixed point.

Discussion

We have presented several results related to the static analysis of BNs through the topological analysis of their interaction graphs. First, we have seen that cycles are necessary to obtain "complex behaviours" (Theorems 3.1 and 3.2). Then, by focusing on positive and negative cycles, we have seen that necessary or sufficient conditions for the presence of several attractors or the presence of cyclic attractors can be obtained (Theorems 3.3-3.10). We have also seen that conditions on cycles permit to obtain results about the comparison of iteration schemes or the model reduction (Theorems 3.11 and 3.12).

For the sake of simplicity, we restricted ourself to the Boolean case and results concerning the global interaction graph. However, Theorems 3.1-3.5, 3.7 and 3.12 have extensions to the non-Boolean discrete case and/or generalizations in terms of local interaction graphs [START_REF] Robert | Iterations sur des ensembles finis et automates cellulaires contractants[END_REF][START_REF] Robert | Les systèmes dynamiques discrets[END_REF][START_REF] Bahi | Convergence of discrete asynchronous iterations[END_REF][START_REF] Shih | Solution of the Boolean Markus-Yamabe problem[END_REF][START_REF] Shih | A combinatorial analogue of the jacobian problem in automata networks[END_REF][START_REF] Remy | Graphic requirements for multistability and attractive cycles in a boolean dynamical framework[END_REF][START_REF] Richard | Necessary conditions for multistationarity in discrete dynamical systems[END_REF][START_REF] Richard | An extension of a combinatorial fixed point theorem of shih and dong[END_REF][START_REF] Richard | Positive circuits and maximal number of fixed points in discrete dynamical systems[END_REF][START_REF] Richard | Negative circuits and sustained oscillations in asynchronous automata networks[END_REF][START_REF] Richard | Local negative circuits and fixed points in non-expansive boolean networks[END_REF].

All the topics regarding the static analysis based on interaction graphs have not been covered in this short survey. In particular, a lot of works have been dedicated to the length of synchronous attractors [START_REF] Goles | Comportement itératif des fonctions á multiseuil[END_REF][START_REF] Tchente | Cycles generated by sequential iterations[END_REF][START_REF] Goles | Dynamical behaviors of Kauffman networks with and-or gates[END_REF][START_REF] Colón-Reyes | Boolean Monomial Dynamical Systems[END_REF][START_REF] German | Large attractors in cooperative bi-quadratic Boolean networks[END_REF][START_REF] Just | Large attractors in cooperative bi-quadratic Boolean networks[END_REF][START_REF] Goles | Comparison between parallel and serial dynamics of boolean networks[END_REF].

With the aim at making tractable the computation of more precise dynamical properties from BNs, recent work investigate the use of other static analysis techniques. For instance, Naldi et al. [START_REF] Naldi | Decision diagrams for the representation and analysis of logical models of genetic networks[END_REF] propose to represent each function g i defined by g i (x) = 1 iff f i (x) = x i with a binary decision diagram; leading to an efficient enumeration of all fixed points of f through algebraic operations on the n resulting decision diagrams. Paulevé et al. developed the Process Hitting framework into which over-approximations of asynchronous dynamics of BNs can be encoded. By using abstract interpretation techniques, very efficient over-and under-approximations of particular reachability properties within Process Hitting models have been built [START_REF] Paulevé | Abstract Interpretation of Dynamics of Biological Regulatory Networks[END_REF][START_REF] Paulevé | Static Analysis by Abstract Interpretation of Biological Regulatory Networks Dynamics[END_REF], having hence applications to reachability (static) analysis within BNs.

Overall, a large variety of properties can be statically derived from BNs. The extension of the presented topological analysis of the interaction graph to other modelling frameworks for biological regulatory networks is yet challenging. While, as mention above, several important results have already been extended to discrete networks, the new degree of freedom brought by discrete modelling prevents the straightforward extension of many results on BNs without imposing constraints on possible logical expressions (such as Theorems 3.6 and 3.8). Finally, the efficient quantification of dynamical properties within hybrid models extending discrete networks with time or stochastic features [START_REF] Siebert | Incorporating time delays into the logical analysis of gene regulatory networks[END_REF][START_REF] Ahmad | Hybrid modelling and dynamical analysis of gene regulatory networks with delays[END_REF][START_REF] Twardziok | Stochasticity in reactions: a probabilistic boolean modeling approach[END_REF] calls for innovative static analysis techniques.

3

 3 Static analysis based on the interaction graph

Theorem 3 . 1

 31 If G(f ) has no cycle, then f has a unique fixed point x, and every path of sig(f ) reaches this fixed point (in at most n steps).

Theorem 3 . 4

 34 If G(f ) has no positive cycle, then aig(f ) has at most one attractor.

Theorem 3 . 5

 35 Let I be a subset of [n]. If every positive cycle of G(f ) has a vertex in I, then aig(f ) has at most 2 |I| attractors.

Acknowledgements. This work has been partially supported by the French National Agency for Research (ANR-10-BLANC-0218 BioTempo project).