B Tanguy

J Besson

R Piques

A Pineau

Ductile to brittle transition of an A508 steel characterized by Charpy impact test. Part-II: modeling of the Charpy transition curve

Keywords: Ductile to brittle transition, Charpy test, Local approach to fracture, Finite element simulation

A finite element simulation of the Charpy test is developed in order to model the ductile to brittle transition curve of a pressure vessel steel. The material (an A508 steel) and the experimental results are presented in a companion paper (Part I [1]).

The proposed simulation includes a detailed description of the material viscoplastic behavior over a wide temperature range. Ductile behavior is modeled using modified Rousselier model. The Beremin model is used to describe brittle fracture. The Charpy test is simulated using a full 3D mesh and accounting for adiabatic heating and contact between the specimen, the striker and the anvil.

The developed model is well suited to represent ductile tearing. Using brittle failure parameters identified below -150 • C, it is possible to represent the transition curve up to -80 • C assuming that the Beremin stress parameter σ u is independent of temperature. Above this temperature, a temperature dependent Beremin stress parameter, σ u , must be used to correctly simulate the transition curve. Quasi-static and dynamic tests can then be consistently modeled.

Charpy V-notch impact tests are still widely used to study the fracture properties of steels. This test was first proposed more than a century ago [START_REF] Russell | Experiments with a new machine for testing materials by impact[END_REF][START_REF] Charpy | Note sur l'essai des métaux à la flexion par choc de barreaux entaillés[END_REF]. The test is simple to perform and has been widely used since then. This testing mode was initially used to determine macroscopic characteristics such as the Charpy fracture energy CVN. Instrumented Charpy devices are now used which allow to measure the whole force-displacement curve from which more information can be gained. However these data are representative of the material and of the specimen so that it remains difficult to transfer these data to actual large structures.

In this paper, a micromechanical analysis of the Charpy test is proposed. The aim is to determine model parameters which are intrinsic to the material. The parameters are then used to simulate the Charpy test so that macroscopic characteristics can be predicted and compared to actual tests in order to validate the model. This strategy allows one to transfer the results of the tests to larger structures with an increased level of confidence.

The first numerical simulation of the Charpy test was proposed by Norris [START_REF] Norris | Computer simulation of the Charpy V-notch toughness test[END_REF] using plane strain conditions. Effects of material rate sensitivity [START_REF] Tvergaard | Effect of material rate sensitivity on failure modes in the Charpy V-notch test[END_REF], of temperature dependence [START_REF] Tvergaard | An analysis of the temperature and rate dependence of Charpy V-Notch energies for a high nitrogen steel[END_REF], of the 3D geometry [START_REF] Mathur | Dynamic 3D analysis of the Charpy V-notch test[END_REF][START_REF] Mathur | 3D analysis of failure modes in the Charpy impact test[END_REF] and of specimen size [START_REF] Benzerga | Size effects in the Charpy V-notch test[END_REF] were studied in a series of theoretical papers lacking of comparisons with experiments.

Such comparisons were carried out by other authors for steels similar to the material investigated in the present study [START_REF] Böhme | Application of micromechanical material models to the evaluation of Charpy tests[END_REF][START_REF] Schmitt | Recent advances in the application of the Gurson model to the evaluation of ductile fracture thoughness[END_REF][START_REF] Rossoll | Determination of the fracture toughness of a low alloy steel by the instrumented Charpy impact test[END_REF][START_REF] Folch | Application of coupled brittle ductile model to study correlation between Charpy energy and fracture toughness values[END_REF]. Simulations of sub-size Charpy specimens were also proposed in [START_REF] Sainte-Catherine | Finite element simulations of Charpy-V and sub-size tests for a low alloy RPV ferritic steel[END_REF] The material used in the present study is fully described in a companion paper [START_REF] Tanguy | Ductile-brittle transition of a A508 steel characterized by the Charpy impact test. Part-I: experimental results[END_REF]. This steel (A508 Cl.3) is used to manufacture the pressure vessel of French nuclear pressurized water reactors. The aging behavior of this material under irradiation is monitored using Charpy V-notch specimens located in capsules inside the pressure vessel. The present work aims at proposing a computational strategy what could ultimately be used to predict the fracture toughness of irradiated materials from Charpy data.

The paper details the material models needed for the simulation of the ductile to brittle transition (DBT). They include the description of: (i) the viscoplastic behavior over a wide temperature range, (ii) ductile tearing, (iii) brittle failure.

The numerical procedures are then described together with the identification of the material parameters. Results of the simulation of ductile crack extension and subsequent brittle failure are finally presented which allow the simulation of the whole Charpy transition curve for an unirradiated material.

Material models

This study was performed on an A508 (16MND5) steel presented in the companion paper [START_REF] Tanguy | Ductile-brittle transition of a A508 steel characterized by the Charpy impact test. Part-I: experimental results[END_REF]. In this section, the different material models needed to simulate the Charpy test are detailed. This includes: (i) the viscoplastic behavior of the undamaged material, (ii) the ductile damage behavior, (iii) a criterion to describe the initiation of brittle failure.

Viscoplastic behavior

Ductile damage in the material is mainly initiated on MnS inclusions whose volume fraction, given by the Franklin formula [START_REF] Franklin | Comparison between a quantitative microscope and chemical methods for assessment of non-metallic inclusions[END_REF] as well as image analysis, is equal to 1.75 10 -4 . As the interface between these inclusions and the surrounding matrix is very weak, the MnS volume fraction can be considered as the initial porosity used in the model for ductile fracture. Further void nucleation may occur on iron carbides (Fe 3 C) for high levels of deformation as shown in [START_REF] Tanguy | Ductile-brittle transition of a A508 steel characterized by the Charpy impact test. Part-I: experimental results[END_REF]. One consequence of the low initial porosity is that coupling between work hardening and void growth can be neglected during tensile tests. These tests can therefore be directly used to identify the viscoplastic behavior of the undamaged material in the transition temperature range.

Experiments conducted at various temperatures (T ∈ [-196 • C, 200 • C]) and strain rates ( ε ∈ [10 -4 s -1 , 5. 10 3 s -1 ]) are detailed elsewhere [START_REF] Tanguy | Ductile-brittle transition of a A508 steel characterized by the Charpy impact test. Part-I: experimental results[END_REF][START_REF] Tanguy | Mechanical behaviour of A508 steel based on double nonlinear viscoplactic constitutive equation[END_REF]. This included quasi-static tensile tests as well as dynamic tests using Hopkinson split bars. The main results are summarized as follows.

The flow stress of the material is expressed as a function of temperature T and plastic strain, p, as:

R(p, T ) = R 0 + Q 1 (1 -exp(-b 1 p)) + Q 2 (1 -exp(-b 2 p)) (1)
where R 0 , Q 1 and b 1 depend on T whereas Q 2 and b 2 are chosen constant. The equivalent plastic strain rate ṗ is given by the viscoplastic flow function ṗ = F(σ -R) which is expressed as:

1 ṗ = 1 F = 1 ε1 + 1 ε2 with εi = σ -R K i n i i = 1, 2 (2)
where σ is the tensile stress. The strain rates ε1 and ε2 are each representative of a deformation mechanism: (1) Peierls friction, (2) phonon drag [START_REF] Campbell | The temperature and strain-rate dependance of the shear strength of mild steel[END_REF]. The mechanism with the smallest deformation rate controls deformation. In practice, the phonon drag mechanism only prevails at very high strain rates (> 1000 s -1 ). In eq. 1, K 1 and n 1 are temperature dependent whereas K 2 and n 2 are constant.

Adjusted model parameters for plastic hardening and strain rate effects are given on fig. 1. Fig. 2 shows that these parameters lead to a good agreement between the experimental conventional yield stress (R 0 ) with the predicted one values. The good agreement is observed over the whole experimental strain rate and temperature ranges which were investigated, as already shown in [START_REF] Tanguy | Ductile-brittle transition of a A508 steel characterized by the Charpy impact test. Part-I: experimental results[END_REF].

Ductile failure

Ductile failure is commonly represented by models using a single damage parameter which represents the void volume fraction or porosity, f [START_REF] Tvergaard | Material failure by void growth to coalescence[END_REF][START_REF] Rousselier | Ductile fracture models and their potential in local approach of fracture[END_REF]. The Gurson model extended by Tvergaard and Needleman [START_REF] Tvergaard | Material failure by void growth to coalescence[END_REF] (so called GTN model) is based on micromechanical considerations and is commonly used. The Rousselier [START_REF] Rousselier | Ductile fracture models and their potential in local approach of fracture[END_REF] model is based on thermodynamics [START_REF] Lemaitre | A course on damage mechanics[END_REF][START_REF] Lemaitre | A continuous damage mechanics model for ductile fracture[END_REF]. However the original formulation of this model is not suitable for temperature and rate dependent materials as shown in [START_REF] Tanguy | An extension of the Rousselier model to viscoplastic temperature dependent materials[END_REF]. To overcome these difficulties, a modified Rousselier model was proposed in [START_REF] Tanguy | An extension of the Rousselier model to viscoplastic temperature dependent materials[END_REF][START_REF] Besson | Modeling of crack growth in round bars and plane strain specimens[END_REF]. This modified Rousselier model will be used in the following. All models assume an additive decomposition of the strain rate tensor ε as:

ε = εe + εp + εth (3)
where ε e is the elastic strain tensor, ε p the plastic strain tensor and ε th the thermal strain tensor. The elastic strain tensor is related to the stress tensor σ by the Hooke's law: σ = C : ε e where C is the elastic fourth order tensor. The dependence of the elastic properties on damage is neglected as in many other studies devoted to ductile rupture. Thermal expansion (caused by adiabatic heating) can also be neglected as it occurs only in areas where plastic strain is very large.

The model is based on the implicit definition of an effective stress (σ ) which is a function of both the macroscopic Cauchy stress (σ) and the porosity:

σ eq (1 -f )σ + 2 3 f D R exp q R 2 σ kk (1 -f )σ -1 def. σ = 0 (4)
where σ eq is the von Mises equivalent stress and σ kk the trace of the stress tensor, while q R and D R are material parameters which need to be adjusted.

The viscoplastic yield function is written as:

φ = σ -R(p). (5)
Viscoplastic flow occurs when φ ≥ 0. The flow stress R is given by eq. 1 as a function of the effective plastic strain, p. The plastic strain rate tensor is given by the normality rule as:

εp = (1 -f ) ṗ ∂φ ∂σ = (1 -f ) ṗ∂σ ∂σ . (6)
Noting that σ is an homogeneous differentiable function of σ of degree 1, Euler's theorem applies so that: (∂σ /∂σ) : σ = σ . Consequently one obtains:

εp : σ = (1 -f ) ṗσ . ( 7
)
This corresponds to the fact that the macroscopic plastic work (left handside) is equal to the microscopic one (right handside) [START_REF] Shima | Plasticity theory for porous metals[END_REF]. The multiplier ṗ is given by eq. 2 as ṗ = F(σ -R). The evolution of the damage variable is governed by mass conservation modified to account for nucleation of new voids [START_REF] Tvergaard | Material failure by void growth to coalescence[END_REF][START_REF] Chu | Void nucleation effects in biaxially stretched sheets[END_REF]:

ḟ = (1 -f )trace( εp ) + A n ṗ. (8)
In this expression, the first right handside term corresponds to void growth and the second one to strain controlled nucleation. A n is a material parameter representing the nucleation rate [START_REF] Zhang | A complete Gurson model approach for ductile fracture[END_REF]. All model parameter directly related to ductile fracture (i.e.

q R , D R and A n ) are assumed to be temperature independent.

Under quasi-static loading, isothermal conditions are assumed. Under rapid loading corresponding to adiabatic conditions, the temperature T increases due to plastic deformation. In the following, the temperature evolution is written as:

C p Ṫ = β εp : σ = (1 -f )β ṗσ (9)
where C p is the volume heat capacity and β a constant factor. Note that C p equals

(1 -f )C 0 p where C 0 p is the heat capacity of the undamaged material so that C 0 p Ṫ = β ṗσ . The value of C 0 p was obtained from the literature [START_REF] Tigges | Nocivité des défauts sous revêtement des cuves de réacteurs à eau sous pression[END_REF].

Brittle failure

Brittle fracture is described using the Beremin model [START_REF] Beremin | A local criterion for cleavage fracture of a nuclear pressure vessel steel[END_REF] which accounts for the random nature of brittle fracture. The model is based on the Weibull weakest link theory. The model uses the following definition of a local effective stress for brittle failure σ Ip :

σ Ip =    σ I exp(-p/k) if ṗ > 0, p > p c 0 otherwise ( 10
)
where σ I is the maximum principal stress of σ. The condition ṗ > 0 expresses the fact than failure can only occur when plastic deformation occurs. The exp(-p/k) term has been proposed in [START_REF] Beremin | A local criterion for cleavage fracture of a nuclear pressure vessel steel[END_REF] to phenomenologically account for grain shape changes or crack blunting effect induced by plastic deformation. p c is the critical strain over which cleavage can occur. This brittle failure model can be applied as a post-processor of calculations including ductile tearing. In this case, care must be taken while computing the failure probability as ductile crack advance leads to unloading of the material left behind the crack front. Considering that each material point is subjected to a load history, σ(t), p(t) (t = time) the probability of survival of each point at time t is determined by the maximum load level in the time interval [0, t]. An effective failure stress σIp is then defined as [START_REF] Tanguy | Numerical modelling of Charpy V-notch tests[END_REF][START_REF] Lefevre | A modified Beremin model to simulate the warm pre-stress effect[END_REF]:

σIp (t) = max t ∈[0,t] σ Ip (t ) (11)
The failure probability, P R , is obtained by computing the Weibull stress, σ w :

σ w = V σm Ip dV V 0 1/m P R = 1 -exp - σ w σ u m ( 12
)
where the volume integral is taken over the whole volume of the specimen. V 0 is a reference volume which can be arbitrarily fixed. Model parameters (σ u , m, k, p c ) must be adjusted. They were considered as constant as in the original Beremin model.

Numerical simulation

The model for ductile rupture has been implemented in a general purpose implicit finite element code [START_REF] Besson | Large scale object-oriented finite element code design[END_REF]. An implicit scheme is used to integrate the constitutive equations. The consistent tangent matrix is computed using the method proposed in [START_REF] Simo | Consistent tangent operators for rate-independent elastoplasticity[END_REF]. Details of the implementation can be found in [START_REF] Besson | Modeling of crack growth in round bars and plane strain specimens[END_REF]. Finite strains are treated using a corotational reference frame [START_REF] Ladevèze | Sur la théorie de la plasticité en grandes déformations[END_REF][START_REF] Besson | Mécanique non-linéaire des matériaux[END_REF] defined so that the stress rate corresponds to the Jauman rate.

Depending on the geometry different element types were used. This includes plane strain (PE) elements used to model the Charpy test in a simplified way (section 6), axisymmetric (AXI) elements used to model round notched bars and 3D elements used to model the Charpy tests. In all cases except in the contact zones in 3D cases where linear elements (i.e. 8 nodes bricks) have been used, caculations are done using quadratic elements with reduced integration. In regions where the crack propagates, 8 nodes (PE and AXI) and 20 nodes (3D) elements are used. The contact between the Charpy sample and the striker and the anvil was modeled using Coulomb friction with a friction coefficient equal to 0.1. The elements in the contact regions are linear elements in the case of 3D simulations. To join the regions, nodes belonging to quadratic elements and lying at the middle of an egde between a quadratic and a linear element are constrained so that they remain at the middle of the edge.

The finite element mesh used to simulate the Charpy specimen is shown on fig. 3. Only 1 4 of the specimen is meshed due to symmetries. The specimen, striker and anvil sizes and geometries are those of the AFNOR standard [START_REF] Afnor | Essai de flexion par choc sur éprouvette Charpy. Partie 1 : méthode d'essai[END_REF].

Tests under dynamic conditions (striker speed v 0 = 5 m.s -1 ) and quasi-static conditions (striker speed v 0 = 1 µm.s -1 ) corresponding to those of the companion paper [START_REF] Tanguy | Ductile-brittle transition of a A508 steel characterized by the Charpy impact test. Part-I: experimental results[END_REF] were simulated. Inertial effects are not accounted for in both cases as they were shown to affect failure only at very low temperatures and not in the transition regime [START_REF] Tvergaard | An analysis of the temperature and rate dependence of Charpy V-Notch energies for a high nitrogen steel[END_REF][START_REF] Tahar | Applications de l'approche locale de la rupture fragile à l'acier 16MND5 : Corrélation résilience-ténacité, Probabilité de rupture bimodale (clivage et intergranulaire)[END_REF][START_REF] Rossoll | Mechanical aspects of the Charpy impact test[END_REF]. Dynamic tests were simulated using adiabatic heating conditions as discussed in [START_REF] Tanguy | Ductile-brittle transition of a A508 steel characterized by the Charpy impact test. Part-I: experimental results[END_REF]. Quasi-static tests were simulated under isothermal conditions.

Once the porosity has reached f c at a given Gauss point, this Gauss point is considered as "broken". In that case, the behavior is replaced by an elastic behavior with a very low stiffness (Young's modulus: 1 MPa). Once a sufficiently high number of Gauss points are broken (4 for 20 nodes bricks, 2 for 8 nodes quadrangles), the element is automatically removed by checking this condition after global convergence for each time increment (therefore the elements are removed in an explicit way).

In the case of brittle fracture, the volume integral in eq. 12 is evaluated using a standard Gauss quadrature. Before computing the maximum principal stress, σ I , the stress tensor is smoothed in order to obtain a constant stress tensor trace in each element. At each Gauss point, the computed stress tensor σ is therefore replaced by the following tensor: σ = σ -1 3 (σ kk -σ kk )1 where 1 is the unity tensor and σ kk the mean stress trace in the element (computed using the Gauss quadrature). This smoothing procedure differs from the one proposed by ESIS [START_REF] Schwalbe | ESIS P6 98: Procedure to measure and calculate material parameters for the local approach to fracture using notched tensile specimens[END_REF] which recommends to compute the Weibull stress using the average stress tensor for each element. The present procedure allows to smooth the pressure variations caused by the quasi-incompressibility of the material (for low levels of porosity) while preserving a more accurate stress profile. In any case, differences between both procedures remain very limited.

Adjustment of the material parameters

The identification of material parameters relative to the viscoplastic behavior is relatively straightforward as tests can be performed on simple specimens in which stress and strain can be considered as uniform. This is not the case in the case of ductile and brittle rupture. The identification procedure is then based on the comparison between the experimental response of a small structure and the simulated one. This simulation must be performed using finite element simulation.

In order to simplify the identification procedure, the ductile and brittle rupture models are identified at temperatures where a single failure mechanism is active. Brittle rupture was identified for T ∈ [-196 • C, -150 • C] and ductile rupture for

T = 20 • C (see fig. 3 in [1]
). Fitted parameters are assumed to remain constant over the whole investigated temperature range [-196 • C, 200 • C]. The validity of the assumption is discussed in section 5.2. On the other hand, parameters describing hardening and strain rate dependence are temperature dependent (section 2.1). Adjusted parameters are gathered in tab. 1.

Ductile failure

Axisymmetric notched bars [START_REF] Mackenzie | On the influence of state of stress on ductile failure initiation in high strength steels[END_REF][START_REF] Besson | Mécanique et Ingénierie des Matériaux -Essais mécaniques -Éprouvettes axisymétriques entaillées[END_REF] are used to identify the Rousselier model. These are small structures which can be rapidly computed as many simulations with different sets of parameters are required by the identification algorithms. Different notch radii are used: ρ = 6.0, 2.4 and 1.2 mm for a minimum sample diameter Φ 0 equal to 6.0 mm (fig. 4-(a)). Tests are performed measuring the minimum diameter variation as explained in [START_REF] Tanguy | Ductile-brittle transition of a A508 steel characterized by the Charpy impact test. Part-I: experimental results[END_REF]. In these samples, the values of the stress triaxiality ratio ( 13 σ kk /σ eq ), which plays an important role in ductile fracture, are close to those encountered in Charpy specimens. In order to obtain a parameter identification representative of ductile crack extension, it is necessary to include in the experimental data base tests where stable crack growth is obtained. This can be done by conducting the tests on notched bars with a controlled diameter variation rate. With a standard axial displacement control, the test usually becomes unstable after crack initiation so that the rapid load drop (fig. 4-(a)), which corresponds to crack growth, cannot be monitored.

As the model parameters are numerous, it is desirable to obtain some of them directly from the metallurgical analysis of the material [START_REF] Besson | Modeling of scatter and size effect in ductile fracture: Application to thermal embrittlement of duplex stainless steels[END_REF][START_REF] Grange | An anisotropic Gurson model to represent the ductile rupture of hydrided Zircaloy-4 sheets[END_REF]. In the present case, the initial porosity is set equal to the MnS volume fraction (i.e. f 0 = 1.75 10 -4 ). Cross section examinations of notched bars and Charpy specimen, show that secondary voids are nucleated on cementite particles (Fe 3 C) at high levels of plastic strain. Nucleation starts for p ≈ 0.5 [START_REF] Tanguy | Ductile-brittle transition of a A508 steel characterized by the Charpy impact test. Part-I: experimental results[END_REF]. The carbide volume fraction is obtained from the chemical composition of the material: 2.3%. This sets the maximum value of the nucleated porosity. Nucleation is assumed to end for p ≈ 1.1 which is the level of plastic deformation in a tensile bar close to the rupture surface. Consequently the nucleation rate parameter A n (eq. 8) is set to 0.038 for 0.5 ≤ p ≤ 1.1 et to 0.0 otherwise. It can be noted that secondary nucleation plays a limited role in notched bars but that it must be accounted for in the case of Charpy specimen in order to correctly model crack initiation at the root of the notch and crack extension at the free surfaces of the samples. It should also be noted that the nucleation law used in this study strongly differs from the one used by other authors for similar materials [START_REF] Rossoll | Determination of the fracture toughness of a low alloy steel by the instrumented Charpy impact test[END_REF][START_REF] Bernauer | Modifications of the Beremin model for cleavage fracture in the transition region of a ferritic steel[END_REF]. In these cases, the phenomenological law proposed in [START_REF] Chu | Void nucleation effects in biaxially stretched sheets[END_REF] was used with material parameters such that nucleation is almost exhausted for p = 0.5, which is not in agreement with the observations reported in [START_REF] Tanguy | Comment on the "Effect of carbide distribution on the fracture toughness in the transition temperature region of an SA 508 steel[END_REF].

Mesh size and mesh design play an important role in calculations involving the simulation of crack initiation and propagation [START_REF] Rousselier | Ductile fracture models and their potential in local approach of fracture[END_REF][START_REF] Tvergaard | Analysis of the cup-cone fracture in a round tensile bar[END_REF][START_REF] Xia | A computational approach to ductile crack growth under large scale yielding conditions[END_REF][START_REF] Liu | Mesh-dependence and stress singularity in finite element analysis of creep crack growth by continuum damage mechanics approach[END_REF]. Due to the softening behavior of the material, mesh size influences the fracture energy and it is often argued that this quantity should be adjusted on experiments involving stable crack growth in order to fit the crack propagation energy [START_REF] Rousselier | Ductile fracture models and their potential in local approach of fracture[END_REF][START_REF] Gullerud | Simulation of ductile crack growth using computational cells: numerical aspects[END_REF]. The mesh size is then interpreted as some microstructural characteristic length (e.g. interparticle spacing or grain size). Experiments on nodular cast iron having the same volume fraction of nodules but different sizes support this interpretation as it was shown that ductility decreases with decreasing interparticle spacing and that this effect can be modeled using finite element sizes proportional to particle spacing [START_REF] Steglich | Micromechanical modelling of damage and fracture of ductile materials[END_REF].

The element size is characterized by three quantities [START_REF] Rivalin | Ductile tearing of pipeline-steel wide plates -II.: Modeling of in-plane crack propagation[END_REF]: h ⊥ the mesh size in the direction perpendicular to the crack plane, h l the mesh size along the propagation direction and h t the mesh size tangent to the crack front (see fig. 5). h ⊥ plays the most significant role on crack propagation. It should be considered as the material parameter to be adjusted. h l must be of the same order of magnitude as h ⊥ in order to consistently capture the stress-strain fields in the propagation direction. In particular it must be small enough to capture the stress maximum which develops ahead of the notch root so that brittle fracture can be correctly modeled. Finally h t is used to describe the gradients along the thickness of the specimen. For this reason, h t is usually smaller at the free boundaries. The values of these mesh size parameters are given in Table 1. The value of 100 µm retained for h ⊥ corresponds roughly to the mean distance between MnS inclusions.

Brittle failure

Model parameters (σ u , m and k) are first identified using tensile round notched bars (NT) tested at low temperature [-196 • C, -150 • C] in order to deal with brittle fracture only. 28 tests were performed. The parameters are identified using an automatic identification procedure which minimizes the quadratic difference between the theoretical and experimental failure probabilities. For a given set of model parameters, the Weibull stress at experimental failure and the associated theoretical failure probability (eq. 12) are computed for each sample. Using the same Weibull stress, the different tests are ordered (increasing Weibull stress) to compute the experimental failure probability as:

P i r = (i - 1
2 )/N where N is the total number of specimens and i = 1 . . . N the rank of the specimen.

The reference volume, V 0 , is taken equal to 0.001 mm 3 . The values of σ u and m were identified for k = 0, 2, 4 (see eq. 10). The best fit was obtained with k = 4. Note however that the influence of k is small as plastic deformation remains limited in the brittle failure temperature range (see fig. 3 in [START_REF] Tanguy | Ductile-brittle transition of a A508 steel characterized by the Charpy impact test. Part-I: experimental results[END_REF]). Experimental and theoretical failure probabilities are compared on fig. 4-(b) as a function of the computed Weibull stress. In this figure, a good agreement between experimental and theoretical values is observed when adopting the following values for the Weibull parameters: m = 17.8 and σ u = 2925 MPa. These values are typical of those found in other studies devoted to cleavage fracture of A508 steels (see e.g. [START_REF] Beremin | A local criterion for cleavage fracture of a nuclear pressure vessel steel[END_REF]).

Results of Charpy V-test simulation

Simulation of ductile tearing

Simulated load-striker displacement curves are compared with experiments on fig. 6 for dynamic1 tests and on fig. 7 for quasi-static tests. The effects of test temperature and striker speed are well reproduced after significant displacement has been applied. The load decrease caused by crack growth is also reproduced.

Load-displacement curves at T = -60 • C and T = 20 • C computed using quasistatic, dynamic with adiabatic heating and dynamic with constant temperature are compared on fig. 8-(a). It can be noted that the isothermal dynamic simulation results in much higher forces than the adiabatic one and overestimates the actual experimental data. On the other hand, under quasi-static conditions, isothermal calculations appear to be adequate. Fig. 8-(b) compares the energy-crack advance curves obtained using the three previous hypothesis. In this figure, ∆a max represents the maximum crack advance which is measured at the center of the specimen. Experimental data at T = -60 • C were obtained from the examination of the fracture surfaces where the transition between ductile crack growth and brittle fracture could easily be observed. In this figure it is noted that a dynamic calculation with constant temperature largely overestimates the Charpy energy for a given crack advance. This result is a direct consequence of the load overestimation for a given displacement. With a more realistic assumption corresponding to adiabating heating, a good agreement is found between tests and simulation. In particular, the slightly smaller energy needed to propagate a crack under quasi-static conditions is well reproduced at -60 • C. This effect can easily be explained by the temperature effect on the constitutive laws of the material. At temperatures above 0 • C, the ductile upper shelf energy (USE) is computed accurately provided that the calculation is carried out for deflection as large as 18 mm. Comparison with experiments are shown on fig. 12-(b) and 12-(d). Some experimental dispersion on the USE can be observed which has been related to material heterogeneities [START_REF] Besson | Modeling of scatter and size effect in ductile fracture: Application to thermal embrittlement of duplex stainless steels[END_REF][START_REF] Decamp | Size and geometry effects on ductile rupture of notched bars in a C-Mn steel: Experiments and modelling[END_REF].

Opening stresses ahead of the crack tip at the center of the Charpy specimen are plotted on fig. 9 for a maximum crack advance, ∆a max , of 1.0 and 2.0 mm at T = -60 • C. The smallest stress maximum is obtained for quasi-static conditions while the highest value corresponds to dynamic/isothermal conditions. It is also shown that the stress maximum is located at about 1.0 mm for the crack tip under quasi-static conditions and only at about 500 µm otherwise. These results show that the mesh size used along the propagation direction h l = 150 µm is large enough to capture the stress field. Stresses tend to slightly increase with increasing crack length. It was noted that the opposite effect is obtained when the crack length is about 4.0 mm because the crack tip is closer to the free surface as this causes a strong constraint loss. Fig. 10-(b) shows crack advance and the crack front shape for ∆a max = 1.0 mm (CVN=97 J) and ∆a max = 3.6 mm (CVN=157 J). The simulation compares well with experiments (fig. 10-(a)): (i) tunneling effect is observed as the crack length is higher at the center of the specimen, (ii) the stress maximum is located at the center of the specimen ahead of the crack tip; this result is in agreement with the position of observed cleavage initiation sites [START_REF] Tanguy | Ductile-brittle transition of a A508 steel characterized by the Charpy impact test. Part-I: experimental results[END_REF]; (iii) the specimen width reduction is well simulated. Note however that slanted shear lips are not reproduced. It is also interesting to note that for large deflections (δ > 7.5 mm) the contact line between the striker and the specimen is no longer located on the plane of symmetry as observed experimentally. Fig. 11 represents the temperature increase, ∆T , at the center of the specimen ahead of the notch tip for two crack advances (2.0 and 3.6 mm). It is shown that the temperature rise is maximum (350 • C) at the notch tip; this is caused by the larger plastic strain needed to reach failure as the stress triaxiality is smaller close to a free surface. Ahead of this zone, the temperature increase at failure lies between 200 and 250 • C. The side opposite to the notch is initially plastically deformed under compression. This leads to a temperature increase which occurs well before this part is actually cracked: this may possibly cause brittle cracks to arrest in this zone. Comparisons with measurements made at the outer free surface of the specimen (fig. 6 in [START_REF] Tanguy | Ductile-brittle transition of a A508 steel characterized by the Charpy impact test. Part-I: experimental results[END_REF]) shows that the simulation overestimates the experiments by a factor 2. On the one hand, this may be caused by experimental problems such as the rapid cooling of the thermocouple at the free surface or the difficulty to put the thermocouple at the location of maximum heat production. On the other hand, this may also be caused by the model itself as failure at the free surface occurs by the formation of shear lips which are not accounted for.

Results: modeling of the Charpy transition curve

In order to determine the Charpy transition curve, the previous calculations including ductile crack extension are post-processed using the Beremin model with the material coefficients optimized in section 4.2. The Charpy energy CVN corresponding to a failure probability of 10%, 50% and 90% is plotted as a function of the test temperature. Results are shown on fig. 12-(a) for dynamic adiabatic conditions and fig. 12-(c) for quasi-static conditions. It is shown that the prediction is satisfactory up to T = -80 • C but that, above this temperature the model largely underestimates the rupture energies. This effect is more pronounced for dynamic tests than for quasi-static tests. This result is nevertheless important as the Weibull parameters were identified in the [-196 • C, -150 • C] temperature range showing that data extrapolation is possible up to -80 • C. Similar results have been reported in the literature [START_REF] Rossoll | Determination of the fracture toughness of a low alloy steel by the instrumented Charpy impact test[END_REF][START_REF] Bernauer | Modifications of the Beremin model for cleavage fracture in the transition region of a ferritic steel[END_REF]. Bernauer et al. [START_REF] Bernauer | Modifications of the Beremin model for cleavage fracture in the transition region of a ferritic steel[END_REF] explained this result considering that void nucleation at carbides eliminates potential sites for cleavage. Therefore they proposed a modified version of the Beremin model to take into account the decrease of the potential defect population. This situation prevails only at high temperature where high plastic strains are needed to reach cleavage stresses. These authors used a nucleation law such that 50% of the carbides have nucleated void for a plastic strain equal to 0.3. This assumption cannot be applied to the present steel since it was experimentally observed that void nucleation at carbides initiates for plastic strains larger than 0.5 [START_REF] Tanguy | Ductile-brittle transition of a A508 steel characterized by the Charpy impact test. Part-I: experimental results[END_REF][START_REF] Tanguy | Modélisation de l'essai Charpy par l'approche locale de la rupture[END_REF].

The transition curve was adjusted using a temperature dependent value for the σ u parameter as already proposed in [START_REF] Tanguy | Numerical modelling of Charpy V-notch tests[END_REF][START_REF] Lefevre | A modified Beremin model to simulate the warm pre-stress effect[END_REF]. This is similar to using a temperature dependent critical cleavage stress as in [START_REF] Margolin | Modeling for ductile-to-brittle transition under ductile crack growth for reactor pressure vessel steels[END_REF][START_REF] Carassou | Modelling of the ductile to brittle transition of a low alloy steel[END_REF][START_REF] Rossoll | Charpy impact test modelling and local approach to fracture[END_REF]. In the case of dynamic tests, as the conventional Charpy V-notch test, the temperature increase ahead of the crack tip should not be neglected when using this solution. Calculations show that the temperature increase, when the stress maximum is reached, lies between 20 • C and 30 • C, in first half of the specimen (i.e. ∆a max < 4.0 mm). This temperature increase is even larger on the specimen side opposite to the notch as the crack propagates in a region which has been predeformed in compression. Consequently the brittle failure model accounting for load history (section 2.3) must be slightly modified. The rupture probability of a volume element is no longer represented by σ Ip but by ω = (σ Ip /σ u ) m where both σ u and m may be temperature dependent and must be evaluated at each location. The load history integrating stress variations but also temperature changes is represented by:

ω(t) = max t ∈[0,t] ω(t ) (13)
Finally the failure probability is given by:

P R (t) = 1 -exp(-Ω(t)) with Ω(t) = V ω(t) dV V 0 . ( 14
)
It is indeed no longer possible to define a Weibull stress and to use experimental results obtained at different temperatures to adjust the model parameters.

The value of σ u was adjusted as a function of temperature in order to match dynamic experiments in the whole temperature range. It was then applied to predict quasi-static tests. The value of m was assumed to be constant and equal to 17. 12-(d) for quasi-static conditions showing a good agreement for both loading conditions. Note that a less good agreement was obtained for dynamic tests when post-processing the calculation with an uniform value for σ u corresponding to the initial test temperature, i.e. not taking into account the local temperature increase (fig. 14). As shown in fig. 12 the values of T K 7 given in [START_REF] Tanguy | Ductile-brittle transition of a A508 steel characterized by the Charpy impact test. Part-I: experimental results[END_REF] (Table 3) are well reproduced using a temperature dependent value for σ u . In particular, the model shows that the DBT temperature is slightly affected by the impact velocity; this results from a complex interaction, which can only be described by a detailed modeling, between stress increase due to viscosity, temperature increase and ductile crack propagation.

As previously noted, some scatter is also observed in the purely ductile failure region. This scatter is not accounted for in the present modeling as it would require numerous simulations integrating material heterogeneities [START_REF] Besson | Modeling of scatter and size effect in ductile fracture: Application to thermal embrittlement of duplex stainless steels[END_REF][START_REF] Decamp | Size and geometry effects on ductile rupture of notched bars in a C-Mn steel: Experiments and modelling[END_REF]. In the transition region where some ductile crack extension occurs, some of the experimental dispersion is likely due to the ductile failure process and not only to the brittle fracture.

Discussion: Simulating the Charpy test

As shown in this study, the Charpy test is a complex test which involves many effects. This complexity makes the Charpy test difficult to simulate and difficult to interpret as a fracture test as most of the effects are not related to damage processes leading to rupture. In the following, these effects are reviewed and discussed in order to determine if they can be neglected:

Inertial effects. Based on previous studies [START_REF] Tvergaard | An analysis of the temperature and rate dependence of Charpy V-Notch energies for a high nitrogen steel[END_REF][START_REF] Tahar | Applications de l'approche locale de la rupture fragile à l'acier 16MND5 : Corrélation résilience-ténacité, Probabilité de rupture bimodale (clivage et intergranulaire)[END_REF][START_REF] Rossoll | Mechanical aspects of the Charpy impact test[END_REF], inertial effects can be neglected as soon as large scale yielding is reached. This is always the case in the transition regime. Only simulations of the Charpy tests below -150 • C would require to account for inertia.

2D/3D calculations. Performing 2D plane strain (PE) simulations instead of 3D calculations represents a tremendous reduction of the computation time in particular in the case of implicit 3D codes. Unfortunately, the difference between both calculations is important as shown on fig. 15. Compared to the 3D calculation, the PE simulations largely overestimate the load (computed assuming a beam thickness of 10 mm) and consequently the Charpy energy and overestimate crack advance. These results are in agreement with other studies [START_REF] Schmitt | Recent advances in the application of the Gurson model to the evaluation of ductile fracture thoughness[END_REF][START_REF] Rossoll | Determination of the fracture toughness of a low alloy steel by the instrumented Charpy impact test[END_REF][START_REF] Eberle | The calculation of dynamic J R -curves from 2D and 3D finite element analyses of a Charpy test using a rate-dependent damage model[END_REF] Accordingly, due to the higher stress level in the fracture process zone the brittle failure probability increases faster using PE calculations, compared to 3D simulations.

Contact.

Neglecting contact can be done by imposing a fixed displacement on the opposite side of the notch and a zero displacement on the initial contact line between the anvil and the specimen. This technique has been used in [START_REF] Tvergaard | An analysis of the temperature and rate dependence of Charpy V-Notch energies for a high nitrogen steel[END_REF][START_REF] Benzerga | Synergistic effects of plastic anisotropy and void coalescence on fracture mode in plane strain[END_REF]. Comparison of PE and 3D calculations with or without contact are shown on fig. 15. Calculations were performed using the same mesh for the Charpy specimen. In the case where displacements are imposed, the boundary conditions are as follows (see fig. 15): u y = 0 at point A and u y = v i t at point B where v i is the impact velocity and t the time. In case of 3D meshes points A and B represent a row of nodes. These boundary conditions slightly differ from those used in [START_REF] Tvergaard | An analysis of the temperature and rate dependence of Charpy V-Notch energies for a high nitrogen steel[END_REF][START_REF] Benzerga | Synergistic effects of plastic anisotropy and void coalescence on fracture mode in plane strain[END_REF] where the displacements v i t is prescribed over a surface and not a line. A significant difference is observed for both plane strain and 3D calculations. Calculations without contact develop higher stresses (and higher stress triaxiality) in the notch plane causing earlier ductile crack initiation and earlier brittle failure. In the case of large striker displacements, the location of the contact zone between the striker and the specimen changes from the center to the side of the striker (fig. 10). This shows that accounting for contact is necessary in particular when a comparison with actual experimental data is needed. Neglecting contact can be used when investigating trends as done in many papers (see e.g. [START_REF] Tvergaard | An analysis of the temperature and rate dependence of Charpy V-Notch energies for a high nitrogen steel[END_REF][START_REF] Mathur | Dynamic 3D analysis of the Charpy V-notch test[END_REF][START_REF] Benzerga | Size effects in the Charpy V-notch test[END_REF]).

Adiabatic heating and heat diffusion.

As shown on fig. 8, the difference between computed loads using isothermal or adiabatic conditions is large as soon as crack growth begins. The difference regarding crack extension and Weibull stress is also important. Consequently, the production of heat by plastic deformation cannot be neglected. Another problem arises in the transition region as the duration t r of the test (about 2 ms) could be large enough so that significant heat diffusion could take place. Typical values for the heat diffusivity D of ferritic steels are equal to 2. 10 -5 m 2 .s -1 [START_REF] Tigges | Nocivité des défauts sous revêtement des cuves de réacteurs à eau sous pression[END_REF]. The corresponding diffusion length is equal to x D = √ Dt r = 0.2 mm which is about the mesh size used in this study. This indicates that neglecting heat diffusion is still possible but that accounting for this phenomenon would slightly improve the accuracy of the simulation but would largely increase the difficulties in the numerical calculations.

Viscoplasticity. Accounting for viscoplasticity is important as the difference between quasi-static and dynamic tests is large. At high impact velocities, stresses are increased, which causes earlier brittle failure (fig. 12) but also increases the ductile upper shelf energy (about 40 J). Both effects are correctly represented by the model developed in this work. As shown in [START_REF] Tanguy | An extension of the Rousselier model to viscoplastic temperature dependent materials[END_REF], the original Rousselier model cannot be used to model damage growth in the case of a rate dependent material. The present study uses a modified version of the model [START_REF] Tanguy | An extension of the Rousselier model to viscoplastic temperature dependent materials[END_REF]. Another solution was proposed by Sainte-Catherine et al. [START_REF] Sainte Catherine | Finite element simulations and empirical correlation for Charpy-V and sub-size Charpy tests on an unirradiated low alloy RPV ferritic steel[END_REF].

Stress history.

In the literature, the calculation of the failure probability after ductile tearing is performed either using the current stress σ Ip (e.g. [START_REF] Moran | Crack growth and cleavage in mismatched welds : a micromechanics study using a cell model[END_REF][START_REF] O'dowd | Prediction of cleavage failure probabilities using the weibull stress[END_REF]) or the maximum stress over the load history σIp (e.g. [START_REF] Tanguy | Numerical modelling of Charpy V-notch tests[END_REF][START_REF] Lefevre | A modified Beremin model to simulate the warm pre-stress effect[END_REF][START_REF] Rossoll | Mechanical aspects of the Charpy impact test[END_REF][START_REF] Bernauer | Modifications of the Beremin model for cleavage fracture in the transition region of a ferritic steel[END_REF][START_REF] Bakker | Prediction of cleavage fracture events in the brittle-ductile transition region of a ferritic steel[END_REF][START_REF] Gao | Analysis of ductile to cleavage transition in partthrough cracks using a cell model incorporating statistic[END_REF]). Computed fracture probability using σ Ip or σIp are compared on fig. 16. Using σIp leads to significantly higher failure probability for a given displacement. In addition for very large crack advance, the probability computed using σ Ip may decrease illustrating the large stress drop that occurs when the crack reaches the free surface. This also clearly shows that this way of computing the failure probability is not suitable. Note that in the case of small scale yielding plasticity, the difference between both methods can remain very small as shown in [START_REF] Ruggieri | A tranferability model for brittle fracture including constraint and ductile tearing effects: a probabilistic approach[END_REF].

Shear lips.

In actual tests, shear lips are created on both free sides of the notch (fig. 10). It remains however difficult to model this phenomenon and a very large number of elements would be required [START_REF] Besson | Modeling of crack growth in round bars and plane strain specimens[END_REF]. The question of the importance of the shear lip formation on the rupture energy still remains an open question.

Conclusions

A model of the behavior of ferritic steels including a description of viscoplasticity, ductile damage and brittle failure has been used to simulate the Charpy transition curve for a pressure vessel steel. An accurate simulation of the load-displacement curve and a good representation of ductile crack growth under both dynamic and quasi-static conditions is obtained using ductile damage parameters identified at 20 Some minor refinements (e.g. heat diffusion, inertia effects) could still improve the simulation. However the main difficulty encountered in interpreting the Charpy test is that a temperature dependent σ u parameter must be used to correctly predict the sharp increase of CVN in the transition regime. Similar solutions have been adopted by other authors (see e.g. [START_REF] Margolin | Modeling for ductile-to-brittle transition under ductile crack growth for reactor pressure vessel steels[END_REF][START_REF] Wallin | Statistical model for carbide induced brittle fracture in steel[END_REF][START_REF] Martin-Meizoso | Modelling cleavage fracture of bainitic steels[END_REF]) to model both CVN and toughness. The micromechanical causes of the need for an increasing σ u parameter remain however unclear. The explanation based on the elimination of cleavage nucleation sites at carbides caused by ductile void nucleation proposed by Bernauer etal. [START_REF] Bernauer | Modifications of the Beremin model for cleavage fracture in the transition region of a ferritic steel[END_REF] cannot be invoked for this material as strain levels at the onset of nucleation are too high. A change in cleavage controlling mechanisms (see e.g. [START_REF] Martin-Meizoso | Modelling cleavage fracture of bainitic steels[END_REF]) or in default population are more likely to be responsible for the observed phenomenon but further detailed investigations are needed.

List of Tables

1 Adjusted model parameters. 24

List of Figures

1 Material coefficients R 0 , Q 1 , b 1 , K 1
and n 1 as functions of the temperature. Other coefficients are constant and equal to: Adiabatic heating β 0.9

Q 2 =472 MPa, b 2 = 1.7, K 2 = 0.18 MPa.s 1/n 2 , n 2 = 1.1 25
C 0 p 3.2 MPa.s -1 n 1 b 1 K 1 Q 1 R 0 Temperature ( • C) n 1 , b 1 R 0 , Q 1 , (MPa) K 1 (MPa.s 1/n 1 ) 60
Q 2 =472 MPa, b 2 = 1.7, K 2 = 0.18 MPa.s 1/n 2 , n 2 = 1.1 +300 • C +200 • C • +100 • C 20 • C + -30 • C ♦ -60 • C -100 • C R 0 (MPa) ε (s -1
) 1000 800 600 400 29

+20 • C 0 • C -30 • C -60 • C -80 • C -100 • C displacement (mm) F (kN)
+20 • C -20 • C -60 • C -80 • C -100 • C -140 • C displacement (mm) F (kN)

C 0 pF

 0 volume heat capacity of the undamaged material (MPa K -1 ) C elastic fourth order tensor D heat diffusivity (m 2 .s -1 ) viscoplastic flow functionF load K i , n i viscoplastic law material coefficients P R total failure probability Q 1 , Q 2 ,b 1 , b 2 plastic hardening parameters R(p, T ) flow stress dependent on cumulative plastic strain and temperature (MPa) R 0 quasi-static yield stress (MPa) T K 7 temperature index corresponding to 70J/cm 2 (56 Joules) for a standard V-notch Charpy specimen V 0 reference volume f void volume fraction f 0 initial void volume fraction h ⊥ mesh size in the direction perpendicular to the crack plane h l mesh size along the propagation direction h t mesh size tangent to the crack front m exponent of the Weibull distibution p equivalent plastic strain of undamaged material (s -1 ) p c critical strain over which cleavage can occur q R , D R modified Rousselier model parameters x D diffusion length δ specimen deflexion ε strain rate (s -1 ) ε strain rate tensor εe elastic strain rate tensor εp plastic strain rate tensor εth thermal strain tensor φ viscoplastic yield surface of the porous material σ tensile stress (MPa) σ effective stress (MPa) σ macroscopic Cauchy stress σ eq von Mises equivalent stress σ kk trace of the stress stensor σ I maximum principal stress of σ σ Ip local effective stress for brittle failure σ u brittle failure model parameter σ w Weibull stress ∆a max maximum crack advance measured at the center of the specimen 1 Introduction

8 .

 8 The adjusted function for the value of σ u is shown on fig. 13. The increase of σ u with temperature implies that the cleavage initiation sites remain the same but are more difficult to trigger. Comparisons with experiments are shown on fig. 12-(b) for dynamic conditions and fig.

Fig. 1 .

 1 Fig. 1. Material coefficients R 0 , Q 1 , b 1 , K 1 and n 1 as functions of the temperature. Other coefficients are constant and equal to:Q 2 =472 MPa, b 2 = 1.7, K 2 = 0.18 MPa.s 1/n 2 ,

Fig. 2 .

 2 Fig. 2. Variation of the yield stress as a function of strain rate for different temperatures (lines: model, points: experiments).

Fig. 3 .

 3 Fig.3. Finite element mesh used for the analysis.

NT 1. 2 NT 2 NT 1 . 2 -Fig. 4 .

 22124 Fig. 4. Identification of the model parameters on notched axisymmetric bars (Φ 0 = 6. mm, Φ = 10.8 mm, r =1.2, 2.4, 6. mm: (a) ductile rupture (line: model, dots: experiment), horizontal thick bars indicate experimental scatter; (b) brittle rupture.

Fig. 5 .

 5 Fig. 5. Mesh design along the crack path.

Fig. 6 .

 6 Fig. 6. (a) Comparison of computed (dots) and experimental (lines) Load-Displacement curves for the dynamic Charpy test at different temperatures. (c) Deformed Charpy specimen at 20 • C for different values of the striker displacement δ.

Fig. 7 .

 7 Fig. 7. Comparison of computed (dots) and experimental (lines) Force-Displacement curves for the quasi-static Charpy tests at different temperatures.

Fig. 8 .Fig. 9 .δFig. 10 .Fig. 11 .Fig. 12 . 3 m = 17 . 8 TFig. 13 .Fig. 14 .Fig. 15 .

 891011123178131415 Fig. 8. Comparison of computations carried out using different conditions (static/dynamic, isothermal/adiabatic): (a) Load-displacement curves; (b) Energy-crack advance curves; dots represent data obtained using low blow tests. Data points are taken from [1] fig. 5.

  • C and a temperature dependent viscoplastic law. Using brittle failure model parameters identified in the fully brittle domain (T < -150 • C) allows the simulation of the lower part of the transition curve up to -80 • C. Above this temperature, the original model gives a pessimistic prediction of fracture energies. It then becomes necessary to use a temperature dependent σ u parameter in the Beremin model to represent the end of the transition curve. A consistent extension of the Beremin model is proposed for this case. Using this modification it becomes possible to represent experimental results under both dynamic and quasi-static conditions which gives some confidence in the applicability of the modified model.

  Temperature increase as a function of the position at the center of the Charpy specimen after a crack advance equal to 2. mm and 3.6 mm. The position is given using the initial (undeformed) coordinates as broken elements are removed during the calculation. The thick line represents points for which the elements are broken.[START_REF] Afnor | Essai de flexion par choc sur éprouvette Charpy. Partie 1 : méthode d'essai[END_REF] 12 Prediction of the Charpy transition curve under quasi-static and dynamic conditions assuming a constant value for σ u and a temperature dependent σ u (USE: ductile upper shelf energy).[START_REF] Tahar | Applications de l'approche locale de la rupture fragile à l'acier 16MND5 : Corrélation résilience-ténacité, Probabilité de rupture bimodale (clivage et intergranulaire)[END_REF] 13 σ u as a function of temperature 3714 Comparison of the Charpy transition curves predicted using the local value of the temperature to compute σ u (solid lines) and using the initial uniform value of the temperature to compute σ u

	10 Fracture at -60 • C under dynamic conditions: (a) Experimental Table 1 Adjusted model parameters. fracture surface, (b) Simulation of ductile crack propagation (contour plots indicate σ I ). Elastic properties Young's modulus E 210 GPa	34
	2 3 4 5 6 7 16 Effect of the failure probability computation mode: Comparison Variation of the yield stress as a function of strain rate for different temperatures (lines: model, points: experiments). Finite element mesh used for the analysis. Identification of the model parameters on notched axisymmetric bars (Φ 0 = 6. mm, Φ = 10.8 mm, r =1.2, 2.4, 6. mm: (a) ductile rupture (line: model, dots: experiment), horizontal thick bars indicate experimental scatter; (b) brittle rupture. Mesh design along the crack path. Load-Displacement curves for the dynamic Charpy test at different temperatures. (c) Deformed Charpy specimen at 20 • C for different values of the striker displacement δ. different temperatures. Force-Displacement curves for the quasi-static Charpy tests at Comparison of computed (dots) and experimental (lines) σ Ip (thin lines) (T = -60 • C). dashed lines) or 3D elements (solid lines) with σIp (thick lines) or of calculations of the rupture probability using plane strain (PE, (a) Comparison of computed (dots) and experimental (lines) Poisson's ratio ν 0.3 Plastic hardening R 0 (T ), Q 1 (T ), b 1 (T ) see fig. 1 Q 2 , b 2 472 MPa, 1.7 Strain rate effect K 1 (T ), n 1 (T ) see fig. 1 K 2 , n 2 0.18 MPa.s 1/n 2 , 1.1, Rousselier model f 0 1.75 10 -4 q R , D R , f c 0.89, 2.2, 0.4 A n 0.038 for 0.5 ≤ p ≤ 1.1 Beremin model V 0 , σ u , 0.001 mm 3 , 2925 MPa m, k, p c 17.8, 4., 10. -4 11 (dashed lines). Mesh size h ⊥ 100µm accounted for. how the boudary conditions are applied when contact is not without contact (thin lines) (T = -60 • C). The sketch illustrates lines) or 3D elements (solid lines) with contact (thick lines) or (Load-displacement) carried out using plane strain (PE, dashed 15 Effect of loading conditions: comparison of calculations h l , h t (3D) 150µm, 800. . . 1200µm	26 27 28 29 30 31 40 38 39
	8	Comparison of computations carried out using different conditions (static/dynamic, isothermal/adiabatic): (a) Load-displacement curves; (b) Energy-crack advance curves; dots represent data obtained using low blow tests. Data points are taken from [1] fig. 5.	32
	9	Stress ahead of the crack tip for ∆a max = 1 mm (a) and ∆a max = 2 mm (b) assuming quasi-static, dynamic isothermal and dynamic adiabatic conditions. Figures indicate the energy (CVN) reached in the simulation. T = -60 • C.	33
		22

In this paper dynamic means at high speed (v 0 = 5m.s -1 ) without taking into account the inertial effects

Acknowledgments

Financial support from Direction de la Sureté des Installations nucléaires (DSIN) and Electricité de France (EdF) is acknowledged. Technical support from B. Marini (CEA-Saclay) and Y. Grandjean (EdF-Chinon) is also acknowledged. Special thanks to Dr. L. Laiarinandrasana for fruitful discussion about numerical simulation.


Laurent Baratchart

Martine Olivi

Fabien Seyfert

A unified approach to Nevanlinna-Pick interpolation problems

This work deals with Complex-valued interpolation by a Schur rational function of given degree at a set of nodes located in the closed lower half-plane, with prescribed maximum points for the modulus (i.e. points where it is equal to 1) on the real axis. The motivation comes from broadband matching, for which the technique we develop offers a new tool.

I. INTRODUCTION

From the seventies on, Nevanlinna-Pick interpolation has been widely studied in connection with various problems in signal processing and control [START_REF] Ball | Robust Control, Multidimensional Systems and Multivariable Nevanlinna-Pick Interpolation[END_REF]. The theoretical approach developed in [START_REF] Sarason | Generalized interpolation in H ∞[END_REF], in which interpolants are viewed as operator dilations, is at the origin of a wide literature including constructive results [START_REF] Zames | H ∞ -optimal sensitivity theory for siso feedback systems[END_REF], [START_REF] Georgiou | A topological approach to Nevanlinna-Pick interpolation[END_REF]. Among the most relevant (and also historical) application domains, we mention sensitivity minimization ( [START_REF] Zames | Feedback, minimax sensitivity, and optimal robustness[END_REF], [START_REF] Zames | H ∞ -optimal sensitivity theory for siso feedback systems[END_REF]), filter design and broadband matching [START_REF] Youla | Interpolation with positive real functions[END_REF].

Pick's theorem yields necessary and sufficient conditions on a set of interpolation data in order that there exists a complex analytic function bounded by 1 in modulus (a socalled Schur function) to meet these conditions on the disc or the half-plane. Composing with the Cayley transform z → (1z)/(1 + z), the problem may be equivalently be stated in terms of Carathéodory functions rather than Schur functions (a Carathéodory function is an analytic function with nonnegative real part).

Nevanlinna's theorem gives a description of all solutions to this interpolation problem [9, Ch. I, sec. 2, Ch. IV, sec. 6]. However, in applications, it is essential to control the degree of the interpolant. Interpolation with a degree constraint was considered in connection with the covariance extension problem in [START_REF] Georgiou | A topological approach to nevanlinna-pick interpolation[END_REF], in which the spectral density plays a key role in controlling the degree. Subsequently in [START_REF] Georgiou | A topological approach to Nevanlinna-Pick interpolation[END_REF], the following interpolation problem was studed: describe all Caratheodory functions Y (z) of degree at most l solving l +1 interpolation conditions Y (z j ) = w j , j = 0, . . . , l, where z 0 , z 1 , z 2 . . . z l lie interior to the analyticity domain and the interpolation values w 0 , w 1 , w 2 . . . , w l are such that the associated Pick matrix P is positive definite; the last condition corresponds to the indeterminate case (i.e. the problem has more than one solution [START_REF] Dym | J-contractive matrix functions, reproducing kernel spaces and interpolation[END_REF], [START_REF] Ahlfors | Conformal invariants[END_REF]). An existence result was proved [START_REF] Georgiou | A topological approach to Nevanlinna-Pick interpolation[END_REF]Theorem 5.2], showing that to each polynomial d(z) of degree at most 2l there is a pair of polynomials (π, χ)

such that Y = π χ is a Carathodory function of degree at most l satisfying the above interpolations conditions, and having d(z) as "dissipation polynomial":

d(z) = π(z)ξ * (z) + π * (z)ξ (z) (1)
where the superscript * indicates the paraconjugate polynomial (see definitions below). In fact, after normalization, the correspondence between dissipation polynomials and rational solutions π/χ is one-to-one. This result was proved in [START_REF] Gusev | A complete parameterization of all positive rational extensions of a covariance sequence[END_REF] for the covariance extension problem when d(z) is strictly positive on the boundary of the analyticity domain, and in [START_REF] Georgiou | The interpolation problem with a degree constraint[END_REF] in the general case where d(z) ≥ 0 there. Later, the theory was further adapted to issues in control and filter design, where freedom in the choice of d(z) can been used to shape the interpolants in order to meet additional design specifications (see e.g. [START_REF] Georgiou | The inverse problem of analytic interpolation with degree constraint and weight selection for control synthesis[END_REF]).

In the present work, we analyze a still more general version of this problem where both the interpolation conditions and the zeros of the dissipation polynomial may lie on the boundary. Our motivation comes from broadband matching, and it is more natural to formulate the problem in terms of the scattering matrix representation of a twoport, which is a 2 × 2 matrix of complex-valued functions of the frequency. Both for technical and practical reasons, we choose to work on the lower half-plane C -rather than the disk or the more classical right-half-plane. In this framework, the paraconjugate of a polynomial p is p * (s) = p(-s). The scattering matrix of a lossless filter is thus a 2 × 2 function, analytic in C -, and assuming unitary values on the real line (real points represent frequencies). Such a matrix is called inner. The filter is finite-dimensional if its scattering matrix is rational, and finite-dimensional filters are those realizable in practice. Hereafter, we identify two-ports with their scattering matrix. The problem of synthesizing the filter S in Figure 1 so that the reflexion coefficient G 11 of the scattering matrix of the global system is smallest possible on some bandwidth is a very old one. It corresponds to the need of conveying energy to the system L in the bandwidth, rather than having the energy bouncing back. When the filter is finite-dimensional, the so-called matching theory of Fano and Youla [START_REF] Carlin | Wideband Circuit Design[END_REF] provides one with a parametrization of all G that can be realized as in Figure 1 for fixed L and varying S of given degree. However, it is still unknown at present how to set up matching filtering characteristics from this parametrization, as soon as the load has degree greater that one. Another approach was proposed by J. Helton [START_REF] Helton | Broadbanding:gain equalization directly from data[END_REF] in an infinite dimensional setting, where the matching problem gets reformulated as a H ∞ approximation problem of Nehari type, the solution to which is elegantly produced in terms of the norm and maximizing vectors of a Hankel operator. However, the "infinite degree" of this optimal filter makes it hardly realizable or even computable in practice.

We propose here an intermediate approach [START_REF] Baratchart | Generalized Nevanlinna-Pick interpolation on the boundary. Application to impedance matching[END_REF] where a finite-dimensional filter response of prescribed degree is being synthesized by imposing matching and stopping frequencies with respect to some general frequency varying load. A classical computation show that the reflection coefficient G 11 of the overall system at frequency ω is

G 11 (ω) = det(S(ω)) S 22 (ω) -L 11 (ω) 1 -S 22 (ω)L 11 (ω) .
A matching frequency, that is, a frequency ω for which G 11 (ω) = 0, satisfies (whenever

|L 11 (ω)| < 1) that S 22 (ω) = L 11 (ω). (2)
In contrast, a stopping frequency ω, defined by the property that

|G 11 (ω)| = 1, satisfies (whenever |L 11 (ω)| < 1) that |S 22 (ω)| = 1 or equivalently that S 12 (ω) = S 21 (ω) = 0. (3)
Writing S (of McMillan degree l) in Belevitch form [START_REF] Belevitch | Elementary applications of scattering formalism to network design[END_REF], normalized so that lim s→∞ S(s) = I 2 (see remark 2.1), we have that

S = 1 q p * -r r * p , (4)
where p, q are monic complex polynomials of degree l and r is a complex polynomial of degree at most l -1, satisfying the Feldtkeller equation:

qq * = pp * + rr * . (5)
If (x 1 , x 2 , . . . , x l ) denotes a set of matching frequencies and if we set γ j = L 11 (x j ) for j = 1, . . . , l, the matching problem amounts to solve [START_REF] Carlin | Wideband Circuit Design[END_REF] whith fixed transmission polynomial r by imposing l -1 stopping frequencies. This problem looks quite similar to the one considered in [START_REF] Georgiou | A topological approach to Nevanlinna-Pick interpolation[END_REF], however there is a major difference: interpolation takes place on the boundary. In [START_REF] Baratchart | Generalized Nevanlinna-Pick interpolation on the boundary. Application to impedance matching[END_REF], existence and uniqueness of the solutions were proved under the restrictive assumption that r has no root on the real axis, which is not suited for the application to broadand matching just described. In this paper, we remove this requirement. Moreover, we consider two different normalizations and deal with mixed type of interpolation conditions (in the analyticity domain and on the boundary). We propose a topological approach that provides an algorithmic way to compute the solutions. The tricky part is to address both interpolation points on the real line and non strictly positive dissipation polynomials. This is possible upon moving back and forth from Schur to Carathéodory functions and using representation formulas in the Hardy space H 2 (C -) for Carathéodory functions.

II. MAIN RESULTS.

Let X = (x 1 , x 2 . . . x m ) be m distinct real points associated to m interpolation values γ = (γ 1 , γ 2 . . . γ m ) in D m , with D the unit disk. Let further Z = (z 1 , z 2 . . . z l ) ∈ (C -) l be l distinct complex numbers associated to l complex interpolation values β = (β 1 , β 2 . . . , β l ). We define P(Z, β ) to be the so-called Pick matrix associated with the interpolation data (z k , β k ), namely the Hermitian l × l matrix defined by:

P k, j (Z, β ) = 1 -β k β j i(z k -z j ) . (6)
We denote by P + Z (or simply P + ), the set of those β ∈ C l such that (Z, β ) is positive definite. The total number of interpolation conditions is thus equal to N = m + l.

We consider the following two problems, each of which corresponds to a specific normalization : Given X ∈ R m , γ ∈ D m , Z ∈ (C -) l , β ∈ P + and r = 0 a complex polynomial of degree at most N -1, such that r(x k ) = 0, k = 1, . . . , m, to find (p, q) a pair of complex polynomials such that, p q

(x k ) = γ k , for k = 1, .., m (7)
p q (z k ) = β k , for k = 1, .., l (8)
qq *pp * = rr *

where Problem P: p and q are monic of degree N and q is stable in the broad sense (no root in the open lower half-plane C -). Problem P: p and q have degree at most N -1, q is stable in the broad sense, and verifies the normalization q(-i) > 0. Remark 2.1: Problem P imposes a normalization at infinity on the interpolant p/q, since p, q are monic of degree N hence (p/q)(∞) = 1. In connection with the matching problem discussed in Section I, where p/q is thought of as the entry S 22 of the scattering matrix S of a filter, this is justified in that the low-pass equivalent model of LC-resonant filters behaves like an open circuit at infinite frequency, which translates into the normalization S(∞) = Id. However, if for example we add a transmission line in front of the filter, this line can be modeled in the narrow band approximation by a reflection coefficient which is unimodular with free phase at infinity. This extra design parameter can be used to meet an additional interpolation condition or, in a dual way, to reduce the degree of p and q while keeping the interpolation properties of p/q. This leads us naturally to problem P.

Let P N designate the set of complex polynomials of degree at most N and PM N the subset of monic polynomials. To each Q ∈ P 2N which is non-negative on R, we associate by spectral factorization a unique polynomial φ(Q) ∈ P N of degree half the degree of Q such that φ(Q)( φ(Q)) * = Q, which is stable MTNS 2016, July 11-15, 2016 Minneapolis, MN, USA in the broad sense and satisfies φ(Q)(-i) > 0. If in addition Q is monic of degree 2N, we put ϕ(Q) ∈ PM N for the corresponding polynomial, normalized this time so as to be monic. We also put N = N -1. Now, for fixed r, Equation (5) associates

• to each p ∈ PM N a unique monic polynomial q = ϕ(pp * + rr * ) ∈ PM N , • to each p ∈ PN a unique polynomial q = φ(pp * + rr * ) ∈ PN. The two spectral factorizations ϕ and φ allow us to define two evaluation maps which play an essential role in the study of problem P and problem P:

ψ : p ∈ PM N →           p(x 1 )/q(x 1 )
. . . p(x m )/q(x m ) p(z 1 )/q(z 1 ) . . .

p(z l )/q(z l )           , ( 10
)
where q ∈ PM N is computed from p using the map ϕ, and

ψ : p ∈ PN →           p(x 1 )/q(x 1 ) . . . p(x m )/q(x m ) p(z 1 )/q(z 1 ) . . . p(z l )/q(z l )           , (11)
where q ∈ PN is computed from p using the map φ. Since |p| 2 ≤ |p| 2 + |r| 2 = |q| 2 on R in both cases, the rational function p/q has no real pole and no pole in C - since q is stable in the broad sense. Thus, by the maximum principle, we conclude that |p/q| ≤ 1 on C -. In addition, since no x k is a root of r by assumption, we have that |p(x k )/q(x k )| < 1 hence p/q is not a Blaschke product. Therefore the Pick matrix associated with the interpolation data (z k , p(z k )/q(z k )) must be positive definite whence ψ and ψ take their values in D m × P + Z . Below, we say that a polynomial p ∈ PN has n zeros at infinity if p has degree N -n. Zeros at infinity are considered to lie on the real line. Hereafter the degree of p is abbreviated as deg p. The main result of the paper may now be stated as follows.

Theorem 2.1: 1) ψ is a homeomorphism from PM N onto the product space D m × P + . If p ∈ PM N has no real root in common with r, then ψ is continuously differentiable in a neighborhood of p with invertible derivative, hence it is a local diffeomorphism at p. 2) ψ is a homeomorphism from PN onto D m × P + . If p ∈ PN has no common real root with r (including at infinity), then ψ is differentiable at p and its differential is invertible, hence it is a local diffeomorphism at p. Remark 2.2: Theorem 2.1 shows that, if no root of r coincides with an interpolation point on the real line, Problems P and P have one and only one solution. If p has a common real root with r, then the associated scattering matrix (4) drops in degree.

We now come to a Proposition that enables us to use continuation techniques in order to practically solve for problem P. We define PM N (r) to be the open subset of PM N comprised of those polynomials that have no common real root with r.

Proposition 2.1:

ψ(PM N (r)
) is an open, dense and connected subset of D m × P + . Suppose that ν 0 , ν 1 both lie in ψ(PM N (r)), and that γ is a continuous path from ν 0 to ν 1 in D m × P + . Then, for every ε > 0 there exists a continuous path γ from ν 0 to ν 1 in ψ(PM N (r)) such that sup

t∈[0,1] γ(t) -γ(t) ≤ ε,
where . designates an arbitrary but fixed norm on

R 2N ∼ C N ⊃ D m × P + .
We also point out an interesting relation between problem P and P. In the statement we write ψ r and ψr to stress the dependency of ψ, ψ with respect to the polynomial r. Proposition 2.2: Suppose ν ∈ D m × P + , and (α k ) is a sequence of real numbers tending to +∞. Let p k = ψ -1 α k r (ν) and write the leading term of φ(α 2 k rr * + p k p * k ) in the form e iβ k , noting that it is unimodular. Then, we have that

lim k→∞ e iβ k p k α k = ψ-1 r (ν).
III. A TASTE OF PROOFS. We give an idea of the proof of Theorem 2.1, part 1). A complete proof may be found in [START_REF] Baratchart | Generalized Nevanlinna-Pick interpolation on the boundary[END_REF]. The proof of part 2) must be adapted, but follow the same lines.

Step 1. The map ψ is continuous at every p and if p has no real common root with r, then ψ is differentiable at p. This mainly relies on analogous properties for the spectral factorization ϕ.

Step 2. ψ is injective.

Step 3. ψ is proper: the pre-image of each compact set in D m × P + is compact.

f -1 : f (Ω) → Ω is continuous,
that is: f is a homeomorphism onto its image. Finally, the fact that the image of PM N is all of D m × P + rests on the properness of ψ which implies that ψ(PM N ) is closed (and open according to the above) in the connected space D m × P + .

Step 5. The map ψ is a diffeomorphism from PM N (r) onto its image.

Our proof of the injectivity of ψ as well as the fact that the differential of ψ is locally invertible in a neighborhood of p ∈ ψ(PM N (r)) relies on representation formulas in the Hardy space H 2 (C -). Every f ∈ H 2 (C -) is the Cauchy integral of the non-tangential limit of its real part:

f (z) = - 1 iπ R ℜ f (t) t -z dt, z ∈ C -, (12)
MTNS 2016, July 11-15, 2016 Minneapolis, MN, USA and the non-tangential limit of its imaginary part is the Hilbert transform of the nontangential limit of its real part [9, Ch. III, sec. 2]:

ℑ f (x) = 1 π lim ε→0 + |x-t|>ε ℜ f (t) t -x dt, a.e. x ∈ R. (13)
Consequently, the nontangential limit of f can be recovered from its real part as

f (x) = ℜ f (x) + i π lim ε→0 + |x-t|>ε ℜ f (t) t -x dt, a.e. x ∈ R. (14)
To use these formulas, we move from the Schur class S to the Carathéodory class C . The map f → (1f )/(1 + f ) maps C to S and back. However, the image of f ∈ S under this map may not be in H 2 (C -). This has to do with the fact that real parts of Carathéodory functions, unlike those of Schur functions, may be Poisson integrals of measures rather than functions, and therefore cannot in general be recovered from their pointwise trace on R.

We now briefly explain how injectivity is proved. Assume that there exist distinct polynomials p 1 (z) and p 2 (z) in PM N such that ψ(p 1 ) = ψ(p 2 ). For j = {1, 2}, let q j = ϕ(p j ), we have

p 1 q 1 (x k ) = p 2 q 2 (x k ), 1 ≤ k ≤ m, (15)
p 1 q 1 (z ) = p 2 q 2 (z ), 1 ≤ l ≤ l. (16)
Since, there are at most deg r real numbers t for which |p j (t)/q j (t)| = 1, we can find a complex number ξ of modulus 1, distinct from -1, such that 1 + ξ p j /q j is never zero on R for j = {1, 2}. Consider the Cayley transforms of the Schur functions ξ p j /q j , that is

1 -ξ p j (z) q j (z) 1 + ξ p j (z) q j (z) = 1 -ξ 1 + ξ +Y j (z), (17)
where

Y j := 2ξ 1 + ξ q j (z) -p j (z) q j (z) + ξ p j (z) (18)
is a Carathéodory function which belongs to H 2 (C -). By a straightforward computation we obtain

ℜ(Y j ) = Y j +Y * j = d g j , (19)
where

d = 2rr * |1 + ξ | 2 , g j = (q j + ξ p j )(q j + ξ p j ) * |1 + ξ | 2 , ( 20
)
are non-negative polynomials, g j being monic. At each real interpolation point x k , we rewrite (15) using (13):

J(x k ) = ∞ -∞ d(t) g 2 (t) -g 1 (t) g 1 (t)g 2 (t) dt t -x k = 0, (21)
where we omitted the principal value in the integral (21) because the integrand is in fact nonsingular. At each complex interpolation point z k , we rewrite ( 16) using ( 12)

I(z k ) = ∞ -∞ d(t) g 2 (t) -g 1 t) g 1 (t)g 2 (t) dt t -z k = 0, (22)
and taking conjugates for any polynomial Q in P 2l+m-1 . Observing now that g 2g 1 vanishes at the x k , and choosing Q(t) = (g 2 (t)g 1 (t))/ ∏ m 1 (tx k ) (note that g 2 and g 1 being monic g 2g 1 has degree at most 2l + m -1), we conclude since the integrand is nonnegative that g 1 = g 2 . But, given d, g j non negative polynomials such that d/g j belongs to L 2 (R), there exist a unique Y j ∈ H 2 (C -) satisfying [START_REF] Zames | H ∞ -optimal sensitivity theory for siso feedback systems[END_REF]. Thus Y 1 = Y 2 , which in turn implies p 1 = p 2 and q 1 = q 2 .

I(z k ) = ∞ -∞ d(

Fig. 1 .

 1 Fig. 1. Filter plugged to a load L with reflection coefficient L 11

  Step 4. ψ is a homeomorphism from PM N onto D m × P + . This result is established using a famous result by Brouwer, known as invariance of the domain [15, chap. 10, sect. 62]. If Ω ⊂ R n is open and f : Ω → R n is continuous and injective, it says that f is an open map, meaning that it maps open sets to open sets. Hence f (Ω) is open and the inverse map

  t) g 2 (t)g 1 (t) g 1 (t)g 2 (t)

				dt t -z k	= 0.	(23)
	We now combine linearly equations (22), (23) and (21) using
	2l +m arbitrary complex coefficients. Putting everything over
	a common denominator yields
	∞ -∞	d(t)	g 2 (t) -g 1 (t) g 1 (t)g 2 (t) ∏ l k=1 |t -z k | 2	Q(t)dt ∏ m k=1 (t -x k )	= 0, (24)

Laurent Baratchart, Martine Olivi and Fabien Seyfert are with Inria, the French National Institute for computer science and applied mathematics, BP 93, 06 902 Sophia-Antipolis Cedex, France,Laurent.Baratchart@inria.fr, Martine.Olivi@inria.fr, Fabien.Seyfert@inria.fr

22nd International Symposium on Mathematical Theory of Networks and Systems July 11-15, 2016. Minneapolis, MN, USA