Jean-Marie Mirebeau

Numerical resolution of Euler equations, through semi-discrete optimal transport

Keywords:

Semi-discrete optimal transport

The optimal transport problem, formulated by Monge in 1781, has been the subject of continued research ever since. The objective is find a map T pushing a distribution of mass µ onto another one ν, and which minimizes the average transport cost, see [Vil08] for an overview. We present in this section the framework of semi-discrete optimal transport, in which one of the measures ν is discrete, consisting of a finite weighted sum of Dirac masses, and the second measure µ is absolutely continuous with respect to the Lebesgue measure.

Semi-discrete optimal transport was introduced by Pogorelov and Alexandrov as a theoretical tool for the study of generalized solutions to optimal transport [START_REF] Abgrall | Construction of Simple, Stable, and Convergent High Order Schemes for Steady First Order Hamilton-Jacobi Equations[END_REF], which are not solutions to the Monge-Ampere equation in either the classical or the viscosity sense. Its numerical implementation yields some of the most reliable, fast and scalable methods for computational optimal transport [START_REF] Aurenhammer | Minkowski-Type Theorems and Least-Squares Clustering[END_REF]Mer11,[START_REF] Lévy | A numerical algorithm for L 2 semi-discrete optimal transport in 3D[END_REF].

We give in §1.1 a brief overview of the alternative computational approaches to optimal transport, and proceed in the following subsections §1.2, §1.3, to the description of the semidiscrete method. A prerequisite is to define the mathematical counterpart of (possibly nonoptimal) mass transport.

1 Figure 1: A map T sends domain X (blue boundary), equipped with a density µ = ρ Leb X (grayscale) onto a finite set Y (red dots). The push-forward ν = T # µ is defined by ν({y}) = µ(T -1 ({y})) (mass of the green region), for all y ∈ Y .

Figure 2: The three paradigms of numerical optimal transport. I: Fully discrete, II: Fully continuous, III: Semi-Discrete. All illustrations of this section courtesy of Q. Merigot.

Definition 1.1 (Push forward of a measure µ by a map T ). Let X, Y be metric spaces. Let µ be a Borel measure on X, and T : X → Y be a Borel measurable map. The push-forward of µ by T is the Borel measure ν = T # µ on Y defined for any Borel subset A ⊆ Y by ν(A) := µ(T -1 (A)).

Mass transport preserves total mass by construction: ν(Y ) = µ(T -1 (Y )) = µ(X) if ν = T # µ, in particular µ ∈ Prob(X) ⇒ T # µ ∈ Prob(Y ) where Prob denotes the collection of Borel probability measures. Note also that T # µ depends linearly on µ, but non-linearly on T . In the following, X, Y denote metric spaces which closed balls are compact. This allows to equip Prob(X) with the weak* topology dual to continuous bounded functions C 0 b (X). All subsets, maps and measures are implicitly assumed to be Borelian in the following.

The optimal transport problem, as originally stated by Monge, involves two probability measures µ ∈ Prob(X), ν ∈ Prob(Y ), and a lower-semi-continuous transport cost c : X × Y → R + ∪ {+∞}. The objective is to find a map T : X → Y pushing µ onto ν and minimizing the average transport cost:

min T #µ=ν ˆX c(x, T (x))dµ(x).
(1)

Under this form, the optimal transport problem may have no solution, or several solutions. Adequate problem relaxations, or structural assumptions on the cost and the measures, make the problem well posed. They also allow for numerical approaches, see the next subsection.

Computational optimal transport

The main inputs to the optimal transport problem are the measures µ ∈ Prob(X) and ν ∈ Prob(Y ). We choose to categorize numerical approaches to optimal transport by the structure of this input data, which always belongs to one of the following two categories:

• Discrete measures, given as a finite sum µ = x∈X µ x δ x of Dirac masses.

• Absolutely continuous measures µ = ρ Leb X with respect to Lebesgue on X ⊆ R d .

Discrete and absolutely continuous probabilities both form dense affine subspaces of Prob(X), so there is in principle no loss of generality in assuming either form. Surprisingly enough, this input form has deep implications on the algorithm structure, the numerical cost, and the generality of the (weak) solutions to the optimal transport problem that can be extracted. A method is said fully discrete (resp. fully continuous) if both input probabilties are discrete (resp. continuous), and semi-discrete otherwise.

I Fully discrete optimal transport, via the linear relaxation.

Kantorowitch [START_REF] Kantorovich | Kantorovich: On the transfer of masses[END_REF] introduced a convex and linear relaxation of the optimal transport problem (1), for which the existence of a minimizer is guaranteed under the assumptions of §1 on the spaces X, Y and the cost c. The relaxed transport cost from µ ∈ Prob(X) to

ν ∈ Prob(Y ) is W (µ, ν) := min{ ˆX×Y c(x, y)dπ(x, y); π ∈ Prob(X × Y ), P X # π = µ, P Y # π = ν}, (2)
where P X , P Y denote projections onto first and second coordinates respectively. Monge's original formulation (1) corresponds to couplings of the special form π = (Id, T ) # µ, where T is the unknown transport map. Kantorowitch's formulation is especially needed when both measures µ, ν are discrete, due to the lack of measure preserving1 maps T : (X, µ) → (Y, ν).

When µ and ν both consist of N weighted Dirac masses, Kantorowitch's formulation (2) amounts to minimizing a linear objective depending on N 2 non-negative variables (characterizing the coupling π), subject to 2N -1 linear constraints. General purpose linear program solvers such as the simplex method, or primal dual interior point methods, are directly applicable. Linear program solvers specialized in this specific instance also exist, such as the Hungarian algorithm, and Bertseka's auction's algorithm.

The performance of these methods is however disappointing on large problem instances N ≥ 5000 Dirac masses due to the quadratic number of optimization variables. Since the support of the optimal π is usually a small subset of X × Y , concentrated along the graph of a deterministic map T , multiscale strategies have been proposed [OR15, Sch15] to eliminate inactive variables. Another approach [BCC + 15] consists in adding a small multiple of the entropy of π with respect to µ ⊗ ν to the linear objective in (2), which makes the problem smooth and strictly convex.

II Fully continuous optimal transport, via the Monge-Ampere PDE formulation.

The PDE approach to optimal transport requires some structural assumptions: X, Y are convex and compacts domains of R d , the cost c(x, y) = 1 2 xy 2 is quadratic, and the measures µ = ρ Leb X , ν = σ Leb Y have continuous densities bounded above and below. Monge's problem (1) then admits a unique solution, of the form T = ∇u where u solves the following Monge-Ampere equation in the sense of viscosity solutions:

     u : X → R, convex, ∇u(X) = Y, det(∇ 2 u(x)) = f (x)/g(∇u(x)), ∀x ∈ Ω.
(3) Each line of this formulation raises challenging discretization problems. First, the constraint of convexity can be discretized in various ways, none of which is particularly simple Finally, the Monge-Ampere operator det(∇ 2 u) can be discretized through standard finite differences [START_REF] Loeper | Numerical solution of the Monge-Ampère equation by a Newton's algorithm[END_REF][START_REF] Benamou | Two numerical methods for the elliptic Monge-Ampere equation[END_REF]. However, it is desirable to preserve its monotony at the discrete level, in the spirit of [START_REF] Kuo | Discrete Methods for Fully Nonlinear Elliptic Equations[END_REF], which requires more complex implementations [START_REF] Brittany | Convergent Finite Difference Solvers for Viscosity Solutions of the Elliptic Monge-Ampère Equation in Dimensions Two and Higher[END_REF][START_REF] Benamou | Monotone and Consistent discretization of the Monge-Ampere operator[END_REF].

Other approaches to optimal transport have been proposed, such as the fluid dynamics formulation of Benamou and Brenier [START_REF] Benamou | The Monge-Kantorovitch mass transfer and its computational fluid mechanics formulation[END_REF]. The next section introduces the semi-discrete approach.

An economic metaphor

We illustrate semi-discrete optimal transport through an economic metaphor. The source domain X ⊆ R d , connected and bounded, is regarded as a large city on which population is distributed according to a density µ = ρ Leb X . The target domain Y ⊆ R d is finite, and is regarded as the location of bakeries in the city, see Figure 3 (left). The function c : X × Y → R + ∪ {+∞} represents the cost for agent x ∈ X to walk to bakery y ∈ Y .

It is natural for potential customers to walk to the closest bakery. The basin of attraction of bakery y ∈ Y is called its Voronoi region, see Figure 3 (center):

Vor(y) := {x ∈ X; ∀z ∈ Y, c(x, y) ≤ c(x, z)}.
However, if the bakeries are allowed to set different prices ψ : Y → R for their products, then people x ∈ X will make a compromise and minimize the sum c(x, y) + ψ(y) of transport cost and purchase price. In this refined model, each bakery y ∈ Y attracts a members of its Laguerre region, see Figure 3 (center)

Lag ψ (y) := {x ∈ X; ∀z ∈ Y, c(x, y) + ψ(y) ≤ c(x, z) + ψ(z)}. ( 4
)
When the cost is quadratic c(x, y) := 1 2 |x -y| 2 , the Voronoi and Laguerre regions are convex polytopes. Their computation is a non-trivial yet classical problem of discrete geometry, for which mature software packages are available such as [cga].

The destination of each customer x ∈ X is given by the price dependent transport map

T ψ (x) := argmin y∈Y c(x, y) + ψ(y). (5)
Agents lying in the intersection of two distinct Laguerre regions Lag ψ (y) ∩ Lag ψ (z) do not have a well defined destination. We thus assume in the following, including Th 1.3, that these interesections are Lebesgue-negligible, hence also µ-negligible, for any price function ψ : Y → R and any distinct y, z ∈ Y . This property holds under mild non-degeneracy assumptions on the cost c, and is satisfied by the quadratic cost c(x, y) := 1 2 |x -y| 2 . The frequentation of the bakeries is encoded in the push-forward

T ψ# µ = y∈Y µ(Lag ψ (y))δ y .
Customers x ∈ X, by their individual cost minimizations, construct an optimal transport plan.

Lemma 1.2. For any ψ : Y → R, the map T ψ is an optimal transport from µ to T ψ# µ.

Proof. Let T : X → Y such that T # µ = T ψ# µ. The price dependent transport plan T ψ has the optimality property:

∀x ∈ X, c(x, T ψ (x)) + ψ(T ψ (x)) ≤ c(x, T (x)) + ψ(T (x)).
Integrating with respect to µ yields

ˆX c(x, T ψ (x))dµ(x) + ˆX ψ(T ψ (x))dµ(x) ≤ ˆX c(x, T (x))dµ(x) + ˆX ψ(T (x))dµ(x). (6)
Denoting by ν ∈ Prob(Y ) the common value of

T # µ = T ψ# µ, we obtain ˆX ψ(T ψ (x))dµ(x) = ˆY ψdν = ˆX ψ(T )dµ(x). (7)
Subtracting from both sides of (6) the equal terms (7), namely the total revenue of the bakeries, yields as announced ˆX c(x, T ψ (x))dµ(x) ≤ ˆX c(x, T (x))dµ(x).

The Kantorowitch functional

We continue with the economical metaphor of §1.2, and consider the situation where each bakery y ∈ Y only has the capacity to cater for proportion ν y of the population, and y∈Y ν y = 1. It turns out that it is possible to fix the prices ψ so that the bakery frequentations T ψ# µ match their capacities ν. Doing so amounts to an optimal transport problem.

Theorem 1.3 (Aurenhammer, Hoffman, Aranov [START_REF] Aurenhammer | Minkowski-Type Theorems and Least-Squares Clustering[END_REF]). Let µ = ρ Leb X be a probability with density on the bounded domain X ⊆ R d , and let ν = y∈Y ν y δ y be a probability on the finite set Y . The following are equivalent:

(a) Finding an optimal transport from µ to ν.

(b) Finding prices ψ on Y such that T ψ# µ = ν.

(c) Maximizing the concave function

Φ(ψ) := y∈Y ˆLag ψ (y) (c(x, y) + ψ(y))dµ(x) - y∈Y ψ(y)ν y . (8)
Interestingly, Monge's formulation (1) regards mass transport T # µ = ν as a constraint and then optimizes the cost, whereas point (b) of Theorem 1.3 suggests exploring a family of costoptimal mappings T ψ until the image measure matches ν. The principles underlying these formulations of optimal transport seem opposite, and indeed Theorem 1.3 is a by-product of the dual formulation of the linear relaxation (2) of optimal transport, identified by Kantorowitch:

W (µ, ν) = sup{ ˆX ϕdµ -ˆY ψdν; ϕ, ψ obey ∀(x, y) ∈ X × Y, ϕ(x) -ψ(y) ≤ c(x, y)}, (9)
where implicitly we also require ϕ ∈ L 1 (µ) and ψ ∈ L 1 (ν). If Y is finite, then (9) amounts to maximize over the bakery prices ψ : Y → R the following quantity, actually equal to (8)

ˆX ϕ ψ dµ - y∈Y ψ(y)ν y where ∀x ∈ X, ϕ ψ (x) := min y∈Y c(x, y) + ψ(y). (10)
Note that ϕ ψ is the optimal choice of ϕ in (9) for a given ψ.

The semi-discrete approach to numerical optimal transport is nothing else than the maximization of the Kantorowitch functional (8). Coordinate-wise increments were the first introduced optimization method [OP89], and have cubic complexity O(N 3 (ln N )/ε) where N := #(Y ) and ε is the desired precision. Gradient based methods such as LBFGS are much more efficient in practice [Mer11, Lév14], and are themselves outperformed by Newton methods [START_REF] Fernando De Goes | Blue noise through optimal transport[END_REF]. Their implementation requires to identify the first and second derivatives of the Kantorowitch functional. Differentiation with respect to the price ψ(y) at bakery y ∈ Y is denoted by ∂ y or ∂/∂ y .

For any prices ψ : Y → R, and for each point x ∈ X not at the intersection of two Laguerre cells, one has

∂ y ϕ ψ (x) = 1 if x ∈ Lag ψ (y), 0 otherwise, hence ∂ y Φ(ψ) = µ(Lag ψ (y)) -ν y . (11)
Maximizing Φ amounts to solving ∇Φ = 0, which is thus equivalent to the discrete Monge-Ampere equation ∀y ∈ Y , µ(Lag ψ (y)) = ν y . The expressions of the second derivatives, first identified in [START_REF] Fernando De Goes | Blue noise through optimal transport[END_REF], involve integrating the density ρ over Laguerre cell boundaries. Assuming a quadratic cost c(x, y) = 1 2 |x -y| 2 , the Kantorowitch functional is twice differentiable, and for all distinct y, z ∈ Y

∂ 2 Φ ∂y∂z = ˆLag ψ (y,z) ρ(s) y -z ds, ∂ 2 Φ ∂y 2 = - z =y ∂ 2 Φ ∂y∂z , (12)
where Lag ψ (y, z) := Lag ψ (y) ∩ Lag ψ (z) and ds denotes (d -1)-dimensional Lebesgue measure. As a result, ∇ 2 Φ is the matrix of a Laplacian, on the graph dual to the Laguerre Diagram: for

v : Y → R v T (∇ 2 Φ)v = - y =z (∂ 2 yz Φ)(v(y) -v(z)) 2 . ( 13
)
This matrix is therefore negative semi-definite. In addition, Ker(∇ 2 Φ) is the space of constant functions on Y , provided µ(Lag ψ (y)) > 0 for all y ∈ Y , and one has the continuity of ρ and connectedness of {x ∈ X; ρ(x) > 0}. Our numerical implementation consists of maximizing the Kantorowitch functional using a damped Newton algorithm. The damping, a reduction in the step size, is chosen so as to enforce ν(Lag ψ (y)) > 0 at each iteration, hence also the invertibility of the Kantorowitch Hessian (13) up to constant functions. The performance of the method is evaluated by computing the optimal transport from the uniform density on the hollowed square X = [0, 3] 2 \ [1, 2] 2 , to a cartesian grid Y of Diracs approximating the the uniform density on the square [0, 1] 2 . Note that the non-convexity of the domain X would be a source of difficulty for PDE formulations of optimal transport (3).

Numerical experiments show that the computation time scales linearly in the cardinality of Y , see Figure 4, and that instances involving hundreds of thousands of Dirac masses can be solved in a few minutes on a modest 2.5Ghz Laptop using a single core. The semi-discrete approach to optimal transport has been extended to transport costs other than quadratic, in particular those subject to the Ma-Trudinger-Wang property [MTW05], and applied to inverse problems arising in geometric optics [START_REF] Machado | Intersection of paraboloids and application to Minkowski-type problems[END_REF].

Numerical resolution of Euler equations

Euler discovered in 1755 the equations, which bear his name, governing the motion of incompressible fluid flows without viscosity. His mathematical notation, extremely modern at the time, closely resembles today's form: on a domain X ⊆ R d , and an interval I of time

∂ t u + (u • ∇)u = -∇p, div u = 0, (14)
where p : I × X → R denotes pressure and u : I × X → R d is the velocity field, which is subject to a slipping condition u • n = 0 on the boundary ∂X. Despite their age, the mathematical analysis of Euler equations is still largely open, see [START_REF] Figalli | Variational models for the incompressible Euler equations[END_REF] for an overview.

Let us introduce the flow s : I × X → X of the velocity field u, which describes the motion of fluid particles in Lagrangian coordinates. It is defined for any time t ≥ 0 by an Ordinary Differential Equation (ODE), parametrized by the particle initial position x ∈ X s(0, x) = x, ∂ t s(t, x) = u(t, s(t, x)).

Depending on the occurence, we often denote s t (x) := s(t, x). A fundamental observation of Arnold [START_REF] Arnold | Sur la géométrie différentielle des groupes de Lie de dimension infinie et ses applications à l'hydrodynamique des fluides parfaits[END_REF] is that the Lagrangian coordinate system (s t ) t∈I formally obeys the equation of geodesics on the infinite dimensional space SDiff of diffeomorphisms of X with unit jacobian, equipped with the L 2 (X) metric, see §2.1. This leads to two natural problems for Euler equations:

• The Cauchy problem, forward in time: given the initial position s 0 and the initial velocity ṡ0 of the fluid particles, find their positions (s t ) t≥0 at all positive times. This amounts to computing the exponential map on the space SDiff. • The boundary value problem: given the initial s 0 and final s 1 positions of the fluid particles, find their intermediate states (s t ) 0≤t≤1 . This amounts to computing a minimizing geodesic on the space SDiff.

This manuscript is devoted to the boundary value problem. Regarding the Cauchy problem, we refer for instance to the recent works of Székelyhidi et al [START_REF] De | The Euler equations as a differential inclusion[END_REF], on which a course was given at the Journées EDP 2015. Let us mention that the tools of semi-discrete optimal transport can also be used for the Cauchy problem [dGWH + 15]. We give in §2.1 more detail on Arnold's formal interpretation of Euler equations as geodesics. Then we describe several relaxations in §2.2, 2.3, and the resulting discretizations including the one §2.4 proposed by Merigot and the author [MM15] and involving semi-discrete optimal transport.

Euler equations as a geodesic equation

In addition to incompressible fluid flows, Euler discovered in 1765 the equations governing the motion of a rigid body [Eul65]. Two centuries later, Arnold unified these two works [START_REF] Arnold | Sur la géométrie différentielle des groupes de Lie de dimension infinie et ses applications à l'hydrodynamique des fluides parfaits[END_REF], identifying Eulerian motions of a fluid and of a rigid object as the geodesics on the group of its configurations, equipped with a suitable left-invariant metric. For a rigid object this is the group SO 3 of spatial rotations with respect to some reference position.

The configurations of a fluid on a domain X are on the other hand described by maps s : X → X, where s(x) ∈ X is the present position of the fluid particle which reference position is x ∈ X. Smooth configurations are described by diffeomorphisms. Fluid incompressibility implies that the Jacobian det(∇s) is unit. In the following X ⊆ R d is a connected bounded domain of unit area, and we let

SDiff := {s ∈ C ∞ (X, X); det(∇S) = 1 on X}.
Arnold [START_REF] Arnold | Sur la géométrie différentielle des groupes de Lie de dimension infinie et ses applications à l'hydrodynamique des fluides parfaits[END_REF] embeds SDiff within the Hilbert space M := L 2 (X, R d ), which norm is denoted

• , and regards it as a submanifold. The energy and length of a path s ∈ H 1 ([0, 1], M) are respectively

E(s) := ˆ1 0 ṡt 2 dt, L(s) := ˆ1 0 ṡt dt, where m 2 := ˆX |m(x)| 2 dx. (15)
Here and below, ṡt denotes the time derivative of s t , hence ṡt 2 M = ´X | ṡt (x)| 2 dx is the kinetic energy of the fluid particles at time t. The geodesic distance between s * , s * ∈ SDiff is defined as

d(s * , s * ) 2 := inf{E(s); s ∈ H 1 ([0, 1], SDiff), s 0 = s * , s 1 = s * }. ( 16
)
By Cauchy-Schwartz's inequality, E(s) ≥ L(s) 2 , with equality if and only s is parametrized at constant speed. The minimal geodesic for ( 16), if it exists, thus also minimizes path length, and has constant speed, which corresponds to the fluid kinetic energy conservation. Following Arnold we derive below the equations of geodesics along SDiff, only formally since SDiff is not closed in L 2 (X, R d ). The tangent space to SDiff at the identity element Id, consists of vector fields v ∈ C ∞ (X, R d ) such that Id +tv has unit jacobian up to second order terms as

t → 0. det(Id +t∇v) = 1 + t Tr(∇v) + O(t 2 ) = 1 + t div v + O(t 2 ).
Hence

T Id SDiff = {v ∈ C ∞ (X, R d ); div v = 0}.
From the group law we obtain for any s ∈ SDiff

T s SDiff = s • T Id SDiff = {v • s; v ∈ C ∞ (X, R d ), div v = 0}.
Recall that a smooth path on an submanifold of a Riemannian manifold is a geodesic iff its acceleration is at each time normal to the submanifold tangent space. A geodesic path (s t ) 0≤t≤1 in SDiff thus satisfies for any divergence free vector field v and any t ∈ [0, 1]

0 = st , v • s t = st • s -1 t , v ,
where •, • denotes the L 2 scalar product For the second identity, we used the fact that s t preserves the Lebesgue measure on X, hence that right composition with s t is an isometry of M := L 2 (X, R d ). On a simply connected domain such as X, a vector field is orthogonal to all divergence free vector fields if and only it is a gradient field. Hence there exists a pressure field

p t ∈ C ∞ (X,
= ∂ t u(t, s(t, x)) + ∂ t u(t, s(t, x))∇u(t, s(t, x)).
Identifying these two expressions of the fluid particles acceleration ∂ 2 t s(t, x), and composing with s -1 t , we obtain Euler equations (12) as announced.

The main defect of this approach is that SDiff is not a submanifold of M := L 2 (X, R d ), in particular because it is not closed, hence geodesics need not exist. A possible remedy is to consider instead groups of diffeomorphisms with Sobolev regularity of index s > d/2 + 1, equipped with the associated Sobolev metric. The exponential map is then well defined and surjective, and minimizing geodesics exist between any pair of points [START_REF] Carverhill | Global solutions of the Navier-Stokes equation with strong viscosity[END_REF][START_REF] Bruveris | On Completeness of Groups of Diffeomorphisms[END_REF]. Of course, geodesic equations with respect to these stronger metrics differ from Euler's equations, and from any known physical model. We present in the next section an opposite approach, which relaxes the set of solutions, rather than strenghtening the metric.

The L 2 -relaxation, and flows as permutations

Diffeomorphisms of unit jacobian admit, as a superset, the collection S of maps preserving the Lebesgue measure on X

S := {s ∈ L 2 (X, R d ); s # Leb X = Leb X },
where Leb X denotes the Lebesgue measure on the domain X. See Definition 1.1 for the pushforward of a measure by a map. The set S is closed in L 2 (X, R d ), and is the completion of SDiff in dimension d ≥ 3. In dimension d ≤ 2 however, rigidity phenomena prevent the density of SDiff in S. For instance SDiff = {Id} if X = [0, 1], whereas S contains ‡the symmetric group of {1, • • • , N } for all N , see below. In dimension d ∈ {1, 2} the use of measure preserving maps can nevertheless be defended by arguing [START_REF] Brenier | Generalized solutions and hydrostatic approximation of the Euler equations[END_REF] that fluid flows on d-dimensional domains X generally arise by neglecting small dimensions in a three dimensional physical domain X × [0, ε] 3-d .

Permutations as measure preserving maps. Assume that the domain X can, for all N ≥ 1, be partitioned into N cubes, identical up to translation. For instance [0, 1] = ∪ N -1 k=0 [k/N, (k + 1)/N ]. Denote by X N the collection of their centers, equipped with the uniform probability measure. Then the collection Σ N of permutations of X N isometrically embeds into S, by permuting the partitioning cubes with their centers, and as N → ∞ it becomes dense in S.

A discretization of Euler's boundary value problem ( 16) based these permutations is proposed in [START_REF] Brenier | A combinatorial algorithm for the Euler equations of incompressible flows[END_REF]. Given a number T of time steps, and s 0 , s T ∈ Σ N approximating some boundary data s * , s * ∈ S, minimize over (s t ) T -1 t=1 the path energy (15) counterpart:

T N 0≤t<T x∈X N |s t+1 (x) -s t (x)| 2 . (18)
Figure 6 reproduces numerical experiments of [START_REF] Brenier | Generalized solutions and hydrostatic approximation of the Euler equations[END_REF] for the one dimensional case X = [0, 1], with boundary conditions s * = Id, s * = 1 -Id. The permutations s t associated to intermediate times 0 < t < T turn out to oscillate rapidly. They do not converge in S, but do so in an even more relaxed setting allowing particles to split and be distributed over the domain X at intermediate times, see the next subsection.

The linear relaxation, and the pressure field

Brenier [START_REF] Brenier | The least action principle and the related concept of generalized flows for incompressible perfect fluids[END_REF], introduced a second relaxation of Euler equations, which is linear, and for which the existence of a minimizer is guaranteed. It involves the space of all possible particles paths

Ω := X [0,1] ,
which is compact when equipped with Tychonoff's product topology. Probability measures µ ∈ Prob(Ω) are called generalized flows. They are said incompressible if the particle positions equidistribute at all times on the domain: ∀t ∈ [0, 1]

e t# µ = Leb X , (19)
where e t : Ω → X : ω → ω(t) denotes the evaluation map at time t. The boundary conditions are also modeled by a marginal constraint, coupling the the initial and end time: Finally, the action is the mean with respect to µ of the kinetic energy of the paths

(e 0 , e 1 ) # µ = (s * , s * ) # Leb X . (20)
E(µ) := ˆΩ ˆ1 0 | ω(t)| 2 dtdµ(ω). (21)
In this form, Euler's boundary value problem becomes the optimization of a lower-semi-continuous2 linear functionnal E(µ), subject to infinitely many marginal constraints(19), (20), which are linear and continuous. This optimization problem may have several solutions. Strikingly however, the dual linear problem has a unique3 solution, called the pressure field p : [0, 1] × X → R, which accounts for the incompressibility constraint. For any minimizer µ ∈ Prob(Ω) of the primal problem, µ-almost every trajectory obeys ω(t) = -∇p(t, ω(t)), an equation which makes sense thanks to the regularity p ∈ L 2 loc (]0, 1[, BV(X)) see [START_REF] Ambrosio | On the regularity of the pressure field of Brenier's weak solutions to incompressible Euler equations[END_REF].

Discretization as a multi-marginal, entropy penalized problem. Let N, T be positive integers, let X N be a sampling of N points equidistributed on the domain X, and let s : X N → X N be a permutation representing an incompressible map of interest. Benamou et al [BCC + 15] discretize the linear relaxation of Euler equations through the following multi-marginal optimization problem, which is strictly convex thanks to the added ε-entropy penalization: minimize over all generalized flows γ ∈ Prob(X T N )

x 1 ,••• ,x T ∈X N γ(x 1 , • • • , x T ) c(x 1 , • • • , x T ) + εγ(x 1 , • • • , x T )(ln γ(x 1 , • • • , x T ) -1), where c(x 1 , • • • , x T ) := n(|x 1 -x 2 | 2 + • • • + |x T -1 -x T | 2 + |x T -s(x 1 )| 2 ), (22)
subject to the marginal constraints

x 1 ,••• ,x T γ(x 1 , • • • , xi , • • • , x T ) = 1
for any 1 ≤ i ≤ N , and any fixed xi ∈ X N , discretizing (20). The large number of variables #(X T N ) = N T makes a direct approach hopeless. However, the dual problem only involves N T variables, instead of N T , associated to the incompressibility constraints and thus representing the pressure field: maximize over all p 1 , • • • , p T :

X N → R x∈X N 1≤t≤T p t (x) - x 1 ,••• ,x T ε exp((p 1 (x 1 ) + • • • + p T (x T ) -c(x 1 , • • • , x T )) /ε). (23)
This functional, which involves a sum over N T terms, can be evaluated using only T products of N ×N matrices, thanks to the additive structure of the cost (22) and the multiplicative properties of the exponential. Figure 6 reproduces numerical experiments of [BCC + 15]. The drawbacks of this approach are its numerical cost, as well as the blurring induced by the entropy penalization, which limits the visibility of the transition from classical solutions of ( 14) to non-deterministic generalized flows minimizing (21).

Polar factorization, and the distance to incompressible maps

We describe a discretization of Euler's boundary value problem, introduced by Merigot and the author in [MM15], which implements flow incompressibility by penalizing the distance to incompressible maps: for any m ∈

M := L 2 (X, R d ) d 2 (m, S) := inf s∈S m -s 2 . ( 24
)
The distance from m to S turns out to be the Wasserstein distance from the image measure of m to the Lebesgue measure. This surprisingly links Euler equations of incompressible fluid flows, to optimal transport which describes a pressure-less gaz [START_REF] Benamou | The Monge-Kantorovitch mass transfer and its computational fluid mechanics formulation[END_REF].

Theorem 2.1 (Brenier's polar factorization [START_REF] Brenier | Polar factorization and monotone rearrangement of vector-valued functions[END_REF]). For any m ∈ L 2 (X, R d ), one has

d(m, S) = W 2 (s # Leb X , Leb X ), (25)
where W 2 is the Wasserstein distance4 for the quadratic cost c(x, y) := |x -y| 2 . Furthermore, if m # Leb X is absolutely continuous with respect to the Lebesgue measure, then there exists a unique convex function ϕ : X → R, and a unique measure preserving map s ∈ S, such that

m = ∇ϕ • s. (26)
Brenier's polar factorization takes it name from the special instance where X is the euclidean unit ball and m : x → M x is a linear function. Let M = ΦS be the polar decomposition of M ∈ GL d R, into the product of a a positive definite matrix Φ ∈ S + d , and an orthogonal transformation S ∈ O n (R). Then ϕ : x → 1 2 x, Φx is convex, s : x → Sx preserves the Lebesgue measure on the ball, and one indeed has (26).

Let P N be a partition of X into N regions of area 1/N and diameter ≈ N -1/d , and let M N ⊆ L 2 (X, R d ) be the subset of maps which are piecewise constant on the partition P N . The push-forward of the Lebesgue measure by some m ∈ M N is a discrete probability measure equidistributed (if m is injective) on the image of m m

# Leb X = 1 N x∈Im(m) δ x . (27)
Evaluating the distance (25), from m ∈ M N to incompressible maps S, thus amounts to a semidiscrete optimal transport problem. We regard this as a numerical commodity in view of recent progress in the scalability, speed and reliability of this method, see §1.3. This functional depends on N (T +1) variables. It implements the constraints of incompressibility, and the boundary values, by penalization as λ 1. We show in [MM15], that minimizers of E(T, N T , λ T ) converge in a weak sense to minimizers of the linear relaxation of Euler's equations §2.3, for suitable sequences T, N T , λ T → ∞. The set S of measure preserving maps is non-convex, thus the distance m ∈ L 2 (X, R d ) → d 2 (m, S) is non convex, hence the optimized functional (28) is not convex either. Nevertheless a minimization using the LBFGS algorithm, and a multiscale (in time) initialization strategy, produced sensible results.

Figures 8 and9 from [MM15] illustrate, for the first time, a sharp result of Brenier [START_REF] Brenier | The least action principle and the related concept of generalized flows for incompressible perfect fluids[END_REF]: consider the linear relaxation of Euler equations on the time interval [0, t max ], and the unique associated pressure field p, see §2.3. If t 2 max ∇ 2 p < π Id, then the relaxation admits a unique minimizer, which is a classical solution to Euler equations. If this threshold is exceeded, then several minimizers may exist, and some or all of them may be non-classical generalized flows µ ∈ Prob(Ω). In that case the support of the flow µ has dimension > d, so that fluid particles follow non-deterministic paths, forming complex and turbulent structures. (Second to fifth row) Generalized fluid flows reconstructed by our algorithm, using the boundary conditions displayed in the first and last column. When t max < 1 we recover the classical flow, while for t max ≥ 1 the solution is not classical anymore and includes some mixing.

Figure 3 :

 3 Figure 3: Left: Domain X (blue boundary), population density ρ (grayscale), and bakeries Y (red dots). Center: Voronoi regions (red boundaries) attached to the bakeries. Right: Laguerre diagram, depending on both the bakeries positions Y and their prices ψ : Y → R.

Figure 4 :

 4 Figure 4: Left: Laguerre cells and their dual graph. Center: Numerical test case setup. Right: Computation time versus number of Dirac masses.

Figure 5 :

 5 Figure 5: The motion of inviscid incompressible fluids admits three formulations, either (I) Eulerian based on the local speed v : [0, 1] × X → R d , (II) Lagrangian based on diffeomorphisms s(t, •) which integrate the speed: ∂ t s(t, x) = v(t, s(t, x)), or (III) relaxed as a superposition of individual particles paths ω ∈ Ω, weighted by a measure µ.

Figure 6 :

 6 Figure 6: Relaxed Euler boundary value problem on X = [0, 1] with boundary conditions s * = Id, s * = -Id. The solution for this particular instance is unique and known [BFS09], it involves non-deterministic generalized flows see §2.3. Horizontal axis: initial particle position x ∈ [0, 1]. Vertical axis: particle position in [0, 1] at some specified time t. Left: discretization based on permutations [Bre08] see (18). Right: convex relaxation as multi-marginal problem [BCC + 15], regularized via entropy penalization, see (23).

Figure 7 :

 7 Figure7: Left: partition P N of the domain X into N regions. Right: illustration of the discretization (28). The maps (m t ) T t=0 belong to a linear subspace M N (black) of L 2 (X), whereas π S (m t ) denotes their projection onto the manifold S (blue).

Figure 8 :

 8 Figure 8: (First row) The beltrami flow in the unit square, shown at various timesteps, is a classical solution to Euler equations. The particles color depends on their initial position.(Second to fifth row) Generalized fluid flows reconstructed by our algorithm, using the boundary conditions displayed in the first and last column. When t max < 1 we recover the classical flow, while for t max ≥ 1 the solution is not classical anymore and includes some mixing.

Figure 9 :

 9 Figure 9: Particles trajectories in the case t max = 1.5, see Figure 8. Left: trajectories originating from inside this red circle seem strongly non-deterministic. Center: trajectories originating from this circle remain grouped. Right: particles around the domain center move little.

Except in the special case where these measures are equidistributed on an equal number of points.

Because path energy ω → ´1 0 ω(t) 2 dt is lower-semi continuous on Ω

Unique up to the trivial invariance by addition of a function depending only on time.

The square root of the optimal transport cost.

[Mer11] Q Merigot. A Multiscale Approach to Optimal Transport. Computer Graphics Forum, 30(5):1583-1592, August 2011.

[Mir15]

Jean-Marie Mirebeau. Adaptive, anisotropic and hierarchical cones of discrete convex functions. Numerische Mathematik, pages 1-47, 2015.

[MM15] Q Merigot and Jean-Marie Mirebeau. Minimal geodesics along volume preserving maps, through semi-discrete optimal transport. arXiv.org, May 2015.

[ A M Oberman and Yuanlong Ruan. An efficient linear programming method for Optimal Transportation. arXiv.org, September 2015.

[Sch15] Bernhard Schmitzer. A sparse algorithm for dense optimal transport. In Scale Space and Variational Methods in Computer Vision, pages 629-641. Springer International Publishing, Cham, May 2015.

[Shn94] A I Shnirelman. Generalized fluid flows, their approximation and applications. Geometric and Functional Analysis, 4(5):586-620, 1994.

[Vil08] C Villani. Optimal Transport: Old and New. 2008.


Q Shi

F Latourte

G Cailletaud

Micro-mechanical failure model of a ferritic-bainitic steel : consequence and interests of a dual-phase polycrystalline model

Keywords: Approche locale de la rupture, Rupture fragile, Micromécanique, Acier bainitique

    Micro-mechanical failure model of a ferritic-bainitic steel : consequence and interests of a dual-phase polycrystalline model Q. Shi, F. Latourte, Georges

Introduction

In the context of pressurized water nuclear reactor operation, the integrity of the reactor pressure vessel (RPV) has to be guaranteed during routine operation and when subjected to accidental loading conditions. For example, the loss of coolant accident is associated to a cold water injection that induces an increase of pressure together with a rapid decrease of temperature. In this context, several integrity assessments of the RPV steel have been proposed, as for example the Master Curve method [START_REF] Wallin | Structural integrity assessment aspects of the master curve methodology[END_REF][START_REF] Lee | Master curve characterization of the fracture toughness in unirradiated and irradiated RPV steels using full-and 1/3size pre-cracked charpy specimens[END_REF] that relies on fracture toughness data for irradiated and un-irradiated material. Because several microstructural changes during irradiation are key in understanding and predicting the plastic response and fracture properties modification with irradiation [START_REF] Bohmert | Correlation between irradiation-induced changes of microstructural parameters and mechanical properties of RPV steels[END_REF][START_REF] Kempf | Effect of lead factors on the embrittlement of RPV SA-508 cl 3 steel[END_REF][START_REF] Lu | The role of irradiation-induced phosphorus segregation in the ductile-to-brittle transition temperature in ferritic steels[END_REF][START_REF] Valo | Copper and phosphorus effect on residual embrittlement of irradiated model alloys and RPV steels after annealing[END_REF], European projects have been proposed to develop multi-scale modeling strategies for the RPV internals and the RPV steel with the objective of predicting the fracture toughness from microstructural initial composition, neutron flux and fluence, and other reference data [START_REF] Massoud | PERFECT : prediction of irradiation damage effects on reactor components: A summary[END_REF][START_REF] Leclercq | Perform 60 -prediction of the effects of radiation for reactor pressurevessel and in-core materials using multi-scale modelling -60 years foreseen plant lifetime[END_REF].

In the multi-scale modeling strategy, one key scale is the crystal plasticity whose different underlying physics, needs, implementations and applications have been reviewed in [START_REF] Roters | Overview of constitutive laws, kinematics, homogenization and multiscale methods in crystal plasticity finite-element modeling: Theory, experiments, applications[END_REF]. Crystal plasticity can be used to derive isotropic von Mises plasticity laws for the pressure vessel using homogenization procedures, based on the single-crystal tensile response obtained from discrete dislocation dynamics and lower scale simulations, predictive by nature. Recent advances in small-scale modelling of body-centered cubic (bcc) iron single crystals have been reported in [START_REF] Monnet | Atomic and dislocation dynamics simulations of plastic deformation in reactor pressure vessel steel[END_REF][START_REF] Queyreau | Orowan strengthening and forest hardening superposition examined by dislocation dynamics simulations[END_REF][START_REF] Monnet | Orowan strengthening at low temperatures in bcc materials studied by dislocation dynamics simulations[END_REF][START_REF] Monnet | Dislocation-dynamics based crystal plasticity law for the low-and high-temperature deformation regimes of bcc crystal[END_REF].

On the fracture mechanics side, most of the recently developed multi-scale models have been inspired by the local approach of Beremin [START_REF] Beremin | A local criterion for cleavage fracture of a nuclear pressure vessel steel[END_REF], in which the failure probability is obtained from a Weibull stress distribution identified from J2-plasticity finite element calculation results. In the case of polycrystalline aggregates, crystal plasticity is used to generate microstructural stress fields that are then compared with a local Griffith stress formulated by means of a surface energy, an average carbide density and a carbide size distribution function [START_REF] Mathieu | A micromechanical interpretation of the temperature dependence of beremin model parameters for french RPV steel[END_REF][START_REF] Libert | Temperature dependant polycrystal model application to bainitic steel behavior under tri-axial loading in the ductile-brittle transition[END_REF][START_REF] Vincent | Stress localization in bcc polycrystals and its implications on the probability of brittle fracture[END_REF]. The latter results have been obtained by considering idealized Voronoi microstructures, and therefore the validity of a micromechanical failure approach should be questionned in view of the complexity of the bainitic microstructure [START_REF] Altendorf | 3D reconstruction of a multiscale microstructure by anisotropic tessellation models[END_REF] and of the relationships between microstructural modeling and calculated stress fields [START_REF] Osipov | Fe modelling of bainitic steels using crystal plasticity[END_REF][START_REF] Nguyen | Micromechanical local approach to brittle failure in bainite high resolution polycrystals: A short presentation[END_REF].

In this context, the objective of this work is to introduce a novel RPV microstructural model taking advantage of recent physically based crystal plasticity constitutive equations, and to evaluate this model in term of brittle failure properties in the low temperature regime.

The paper is structured as follows. First, the RPV material is described. Then, a crystal plasticity model is tested for pure iron and identified for RPV steel according to two different modeling possibilities : single phase material or dual-phase material description. Finally,

Mn Ni Mo

Si C Cr Cu Al Co P S V Weight % 1.37 0.7 0.5 0.24 0.159 0.17 0.06 0.23 <0.01 0.005 0.008 <0.01 Table 1: Chemical composition of H1BQ12 material [START_REF] Diawara | Caracterisation quantitative de la microstructure de l'acier 16MND5 des cuves de REP[END_REF] these two models are compared in term of failure probability results.

Material description 1.Large scale heterogeneities due to manufacturing

The RPV is a large scale component manufactured from forged ferrules of 280mm thickness made from low carbon steel 16MND5 (ASTM A508-3) whose chemical composition is provided in Table 1. This composition has been obtained for a reference material considered in this study, designed as H1BQ12 material and sampled from ferrule B's knockout H1 of the pressure vessel Q12. The pressure vessel height is 12.6m, and its diameter is 4.4m.

At room temperature, the thermodynamic equilibrium of 16MND5 steel provides a balance of pearlite and ferrite. In order to obtain a bainitic microstructure composed of bainite and ferrite, a heat treatment is applied following the 3-step recipe:

• Two austenitizations at a temperature between 865 • C and 895 • C followed by a water quench. Each austenization lasts 4 hours and 40 minutes. After this heat treatment, the steel exhibits a tempered bainitic microstructure. Starting from an austenite grain, ferrite laths grow from the grain boundaries during the cooling. Since carbon is much less soluble in ferrite than in austenite, it will start to diffuse from supersaturated ferrite to austenite. This enriched austenitic phase will then lead to the development of cementite at the interface between ferrite and austenite. The precipitation of carbides will therefore diminish the carbon concentration of austenite which will be able to transform into ferrite. If the cooling is slow enough, carbide precipitates will only precipitate at lath boundaries and in the end we will obtain upper bainite : the ferritic phase will be free from cementite. However, if the cooling is quick, some carbon atoms contained into supersaturated ferrite might not have the time to diffuse to the interface with austenite and will precipitate inside the ferrite laths leading to lower bainite. These phenomena are described in [START_REF] Bhadeshia | Bainite in Steels -2nd Edition[END_REF].

The tempering treatments are used to decompose the austenite that may remain between the laths of ferrite and to transform carbides into cementite. The tempering will decrease the steel hardness an in turn its brittleness. The annealing reduces remaining internal stresses.

In 16MND5 steel, pro-eutectoïd ferrite can nucleate at high temperature from prior austenitic grain boundaries for slow cooling conditions by diffusion processes. These equiaxed grains exhibit low carbon content (less than 0.02 wt.%) resulting in a very fine carbide precipitation.

Let us insist on the key microstructural average dimensions for 16MND5 steel. The austenitic grain size is 50 µm, the baintite lath is about 10 µm long and 1 µm thick, the size of the cementite precipitates at the grain boundaries is 1 µm in length and 0.2 µm in thickness and the cementite precipitates within the ferrite laths are spheric with a diameter of 0.1 µm [START_REF] Naamane | Etude de la déformation plastique de la ferrite à basse température : simulation de dynamique des dislocations[END_REF]. Pro-eutectoïd ferrite grains are 10 µm in diameter. Given the production stages of 16MND5, and the sensitivity of the tempered bainitic microstructure to the process temperature and the cooling rate, the chemical composition, phase distribution and mechanical property are not uniform in the component thickness. As cooling rate is much slower in the core of the ferrule, more pro-eutectoid ferrite forms in this region, resulting in lower carbon content, lower tensile yield stress and strength. The variation in carbon content throught the thickness is shown in Figure 1 Table 2: Phase distribution at different sampling depths for H1BQ12 material [START_REF] Diawara | Caracterisation quantitative de la microstructure de l'acier 16MND5 des cuves de REP[END_REF] and H3BQ3) obtained for two different pressure vessels. The variations of yield stress and ultimate tensile stress at room temperature are shown in Figures 2 and3, respectively. A large campaign of EBSD analyses on H3BC18 steel has allowed to characterize the fraction of phases present in the microstructure for three different sampling depths s d : inner surface (s d ≈ 0 mm), quarter thickness (s d = 70 mm) and half thickness (s d = 140 mm) from the inner surface [START_REF] Diawara | Caracterisation quantitative de la microstructure de l'acier 16MND5 des cuves de REP[END_REF]. The main outcomes of this study has been summarized in Table 2. The chosen 16MND5 material, H1BQ12, has been used in the following to validate the simulations and identify the model parameters. This material has been used in several PhD works [START_REF] Renevey | Approches globale et locale de la rupture dans le domaine de la transition fragile-ductile d'un acier faiblement allié[END_REF][START_REF] Hausild | Transition ductile-fragile dans un acier faiblement allié[END_REF][START_REF] Mathieu | Analyse et modélisation micromécanique du comportement et de la rupture fragile de l'acier 16MND5 : prise en compte des hétérogénéités microstructurales[END_REF][START_REF] Rossol | Détermination de la ténacité d?un acier faiblement allié à partir de l"essai Charpy instrumenté[END_REF] and is therefore a good reference material candidate.

C content (wt%)

The H1BQ12 reference tensile specimens were taken at a depth of three quarters of the thickness (s d = 210mm). The tensile experiments were carried out at a strain rate of 4 × 10 -4 s -1 at temperatures ranging from -196 • C to 0 • C. The experiments on this material have been summarized in [START_REF] Renevey | Approches globale et locale de la rupture dans le domaine de la transition fragile-ductile d'un acier faiblement allié[END_REF], and the tensile curves are gathered in Figure 4. An increase of yield stress is observed when the temperature decreases, while the hardening is slightly decreasing.

plastic flow behaviour of pure iron 2.1 DD_CC constitutive law

The DD_CC (Dislocation Dynamics for body Centered Cubic crystals) relation is a law introduced in [START_REF] Monnet | Orowan strengthening at low temperatures in bcc materials studied by dislocation dynamics simulations[END_REF] and [START_REF] Monnet | Dislocation-dynamics based crystal plasticity law for the low-and high-temperature deformation regimes of bcc crystal[END_REF] which is based on physical properties of near pure iron through dislocation dynamics simulations, extending earlier physical models only considering a low temperature mechanism [START_REF] Rauch | Etude de l'écrouissage des métaux -Aspects microstructuraux et lois de comportement[END_REF]. The dislocation density on each slip system, ρ s is an essential internal variable of the DD_CC law. This law physically describes the plastic slip of the 12 {110} slip systems which are known as the most active among the 48 possible systems {110},{112},{123}. For each system the slip direction is noted m s and the normal to the slip plane is n s .

As usual, the single crystal law is written as a relationship between the applied shear stress τ s on a system s and the slip γ s , related to the local stress tensor σ and plastic strain tensor ε p by the Schmid tensor (or orientation tensor) µ:

τ s = σ : µ s , ε = Σ γs µ s , µ s = 1 2 (m s ⊗ n s + n s ⊗ m s ) (1)
The law has been implemented in the finite element code Code_Aster and aims at modeling two main deformation mechanisms at low and high temperature presented in the following.

Low temperature

At low temperature, a predominant screw dislocation network is observed, moving according to a double kink mechanism.

At low temperature the slip rate γs nuc can be expressed as a function of the effective shear stress, the average length of screw dislocations l s , the mobile dislocation density ρ s m , the absolute temperature T and the activation energy for double-kink nucleation under a given activation energy ∆G:

γs nuc = ρ s m bH • l s ρ, τ s ef f • exp   -∆G τ s ef f kT   τ s |τ s | (2)
where H is a constant frequency factor, b stands for the burgers vector magnitude and k is the Boltzmann constant. In the case of iron, direct analyses of experimental loading curves [START_REF] Naamane | Low temperature deformation in iron studied with dislocation dynamics simulations[END_REF] provided the following expression that was found to cover the entire stress range observed in the thermal regime:

∆G τ s ef f = ∆G 0     1 - τ s ef f + τ 0     (3)
where x + = 1 2 (|x| + x) stands for the positive part of x,and the index s relates to each slip system s.

High temperature

At high temperature the plastic deformation is associated to a drag mechanism that can be modeled by means of a Norton law:

γs drag = γ0 |τ s app | τ s c n • τ s app |τ s app | (4)
where γ0 is a constant, τ s c is the critical shear stress and n is the viscous exponent.

General case

In the general case, there is a competition between the two regimes and the total slip can be expressed as an harmonic average of γs drag and γs nuc , giving weight to the most limiting slip mechanism:

1 γs = 1 γs nuc + 1 γs drag (5)
To resolve the low temperature deformation given by equation ( 2) some additional relations define the screw segment length:

l s = max(λ s -2α s AT R s ; l c ) with l c = 500b (T /T 0 ) (6)
(λ s + D obs ) -1 = min ρ s f ; (D obs + 2R s )ρ s f with ρ s f = j =s ρ j , (7)
R s = λ s 2α s AT ; α s AT = j =s a sj AT ρ j ρ s f , ( 8
)
where l c is the critical length of a screw dislocation, λ is the average distance between two obstacles, D obs is the average obstacle diameter, a AT is the interaction matrix at high temperature, ρ s f is the forest dislocation density which is used to take into account the irradiation defects and finally T 0 is a reference temperature: 300 K. The interaction matrix of the 12 considered <111>{110} slip systems has been obtained from dislocations dynamics calculations [START_REF] Queyreau | Slip systems interactions in α-iron determined by dislocation dynamics simulations[END_REF] and is reduced such that only two independent parameters are considered: a AT coli for the collinear systems which interact strongly and a AT =coli for the non-collinear systems of weaker interaction.

Another important variable to resolve the equation ( 2) is the effective shear stress given by:

τ s ef f = τ s app -τ s c , with τ s c = τ s LR 2 + τ s LT 2 + τ F ( 9
)
τ s LR = µb a ss AT ρ s ; τ s LT = max 0 ; α s AT µb 1 λ s - 1 2α s AT R s + l c ( 10
)
where τ s app is the applied shear stress, τ s c is the critical shear stress, τ F is a term accounting for the network drag, τ s LR is a self hardening term and τ LT is related to the line tension. The dislocation density evolution depends on the dislocations mean free path and requires the average diameter of the bainitic grains and two multiplication parameters called K self and K f .

ρs = 1 D lath + √ a ss ρ s K self + α s λ s ρ s f K f -y s ρ s γs b ( 11
)
a sj = a sj AT 1 - τ s ef f τ 0 2 ; 1 y s = 1 y AT + 2πτ s ef f µb ( 12
)
K self and K f are related to each other by the following equation:

K f = K self β , β = min 3 ; 1 + 2 • T T 0 (13)

Tensile behaviour calculation technique

This model requires 15 different parameters but most of them can be determined using physical measurements like the average obstacle length, initial dislocation density, average lath dimension. Table 3 presents the values of these parameters identified in [START_REF] Monnet | Dislocation-dynamics based crystal plasticity law for the low-and high-temperature deformation regimes of bcc crystal[END_REF] on pure iron at temperatures ranging from -196

• C to 0 • C. b (nm) H (s -1 ) ∆G 0 (eV) τ 0 (MPa) ρ m (mm -2 ) ρ ini (mm -2 ) γ0 (s -1 ) n 0.
D obs (nm) y AT (nm) τ F (MPa) D lath (mm) (a sj AT ) coli (a sj AT ) =coli K self 10 2 0 0.01 0.7 0.1 100
Table 3: Parameters used in DD_CC law for pure iron [START_REF] Monnet | Dislocation-dynamics based crystal plasticity law for the low-and high-temperature deformation regimes of bcc crystal[END_REF] With the reference parameters, the strain-stress curve has been calculated first using a selfconsistent Berveiller-Zaoui homogenization technique [START_REF] Berveiller | An extension of the self-consistent scheme to plasticallyflowing polycrystals[END_REF] and compared with experimental results.

When considering a self consistent model, each grain of the polycrystal is successively regarded as an inclusion within the matrix of the homogeneous equivalent medium (HEM). Each grain interacts therefore with the equivalent medium of the whole volume considered. The Berveiller-Zaoui rule assumes that elasticity is uniform and isotropic, that the grains are spherical and that the loading is monotonic. The localization equation writes:

σ g = Σ + 2αµ (1 -β) E p -ε p g , ( 14
)
where g is a given grain, σ g is the average stress tensor on g, Σ is the macroscopic stress tensor, µ is the elastic shear modulus, α is a parameter that can be approximated by equation [START_REF] Mathieu | A micromechanical interpretation of the temperature dependence of beremin model parameters for french RPV steel[END_REF], β is a parameter described in Kröner's model [START_REF] Kröner | Zur plastischen verformung des vielkristalls[END_REF] and given by equation ( 17), E p is the macroscopic plastic strain and ε p g is the average plastic strain on g.

α -1 ≈ 1 + 3 2 µE p Σ -1 , with: (15)
(E p ) 2 = 3 2 E p : E p , Σ 2 = 3 2 S : S , S = Σ - 1 3 Tr(Σ) ( 16
)
β = 2 (4 -5ν) 15 (1 -ν) ≈ 0.5 for ν = 0.3 (17)
From the approximation in equation ( 17), equation ( 14) now writes:

σ g = Σ + αµ E p -ε p g , (18)
with the following homogenization rule:

E p = g x g εp g , (19)
where x g is the volumetric fraction of grain g. The equations ( 18) and ( 19) are the equations solved by Code_Aster that was used to calculate the homogenized tensile response. An aggregate of 100 grains, of equal volume, and giving an isotropic texture, was considered. For pure iron, the complete tensile curve is only available at room temperature and a good agreement is observed, except for the very early hardening stage (ε < 0.015) since the DD_CC law does not yet physically describe the complex transition from elastic elastic regime to microplasticity and plasticity. By sake of synthesis, this tensile response is not reported here. Some validations also not reported herein consisted in comparing the tensile behaviour obtained by homogenization with the one obtained by a finite element calculation of a 98 grain Voronoi polycrystal as already developed in [START_REF] Barbe | Intergranular and intragranular behavior of polycrystalline aggregates.part 2: Results[END_REF], and the gap was found to be less than 5 MPa for 3% of macroscopic strain, for various orientation samplings and various tensile directions applied to the homogenized model.

In the following, the homogenized model response will be plotted against the experimental tensile yield stress from -196 • Cto 0 • C.

Parameter sensitivity study

The following sensitivity study has been conducted for the main adjustable parameters. The parameters were modified one by one, leaving the other parameters fixed. This procedure allows to identify the mechanical role of each parameters and to provide some complementary validation of the reference value since the same experimental data as used in the original identification [START_REF] Monnet | Dislocation-dynamics based crystal plasticity law for the low-and high-temperature deformation regimes of bcc crystal[END_REF], obtained from [START_REF] Cottu | Etude en traction de l'adoucissement du fer de haute pureté par le carbone à basse temperature[END_REF] is shown in the cross plots.

Because the conventional yield stress σ Y 0.2 is obtained during the early stage of deformation, not well captured by the model, the sensitivity analysis is conducted by considering From the results gathered in Figure 2.3, the following conclusions can be drawn:

ρ ini =2.E8mm -2 ρ ini =1.E8mm -2 ρ ini =5.E7mm -2 ρ ini =2.E7mm -2 ρ ini =1.2E7mm -2 exp (a) (b) (c) (d)
• Figure 2.3a show that reducing the lath size D lat significantly increases the yield stress Σ 0.02 . This increase is rather uniform over the temperature range,

• Figure 2.3b show that increasing the dislocation density ρ ini significantly increases the yield stress, more specifically in the high temperature regime (T > -100 • C). The highest dislocation densities correspond to a work hardened material.

• Figure 2.3c and 2.3d show the effect of average obstacle size D obs for low and high initial dislocation density respectively. The low value of ρ ini is 7.2 • 10 6 mm -2 , and the high value is 10 8 mm -2 . While no effect is obtained for low initial dislocation density, a strong yield stress increase is observed for high initial dislocation density. As in Figure 2.3b this effect is more pronounced in the high temperature regime.

From the results of Figure 2.3, the yield stress increase as a function of the lath size as been plotted at room temperature against the square root inverse of D lat . For values above 5 (i.e. for D lat < 5µm) the evolution is linear as shown by the green solid fitting line. This fit corresponds to the following Hall-Petch relationship: In addition, the obtained slope A lies between the experimentally obtained values for interstitial free steel (0.14MPa • m 1/2 ) and 0.03wt% carbon steel (0.81MPa • m 1/2 ) [START_REF] Lim | Simulation of polycrystal deformation with grain and grain boundary effects[END_REF].

σ Y = σ 0 + A • D -1/2 lat , σ 0 = 85.4 MPa, A = 0.25MPa • m 1/2 (20)
Equation 11 is showing that K self is influent on the dislocation multiplication rate. This is confirmed by figure 7

K self =30 K self =60 K self =90 K self =120 K self =150 exp (
-Σ 0.02 (MPa) K self =30 K self =60 K self =90 K self =120 K self =150 exp (b)

Adaptation to more complex materials

The model has been validated for the so called pure iron, which is a large grain polycristalline material with low carbon content. To now adapt the model to a more complex material of bainitic microstructure and of higher carbon content, therefore of higher yield stress, the principle is as follows:

• increase the dislocation density ρ ini and eventually provide measured values for 16MND5. Note that this density is sensitive to the manufacturing (hot/cold work, heat treatment) and that measures rely on transmission electron microscopy analyses on thin foils and therefore are difficult to achieve and need to be repeated enough to be statistically representative,

• increase the obstacle size D obs to the maximal obstacle size present in the microstructure, obstacles being carbides and dislocation junctions,

• decrease the lath size D lat to the minimal mean free path for dislocations, e.g. grain dimension or lath thickness.

• eventually increase the alloy friction parameter τ F to integrate both the solid solution hardening and the precipitation strengthening.

Starting from these general principles, the work conducted on 16MND5 material will be discussed in the following.

plastic flow behaviour of RPV steel

A specific strategy has been introduced in this work to model the RPV steel. The final objective of this work being to use polycrystalline stress fields to compute a failure criterion, an accurate micro-structural description is desired. As already introduced, depending on the sampling depth several phases coexist, the softer being pro-eutectoid ferrite, varying in volume proportion from 0% to 25%. It is then chosen to model the RPV steel as a dual phase model composed by ferrite (soft) and bainite (hard). To start with a simple model, the carbide clusters and martensite are considered as similar to the bainite. Using this approach, varying the ferrite proportion for example allows to capture the variations in yield stress and tensile strength with the sampling depth as reported experimentally in Figures 2 and3.

To illustrate this model, a realization example is shown in Figure 3. A 162 grain aggregate was generated first by a Voronoi tesselation and meshed using the open source Salome software embedding Netgen meshing procedure. A phase label is then attributed to the different grains to reach the desired ferrite fraction, by randomly picking ferrite grains. The ferrite fraction for this work has been set at 25%, representative of the ferrule half thickness.

In the following, a smaller aggregate of 98 grains was used for the finite element calculations. Some informations on this mesh are provided in table 4. For the tensile loadings of this aggregate, the boundary conditions are imposed as average displacement conditions on the upper and lower faces (without planar condition) and the four lateral faces are let totally traction free. A list of parameters specific to each phase was identified by conducting tensile simulation on a dual-phase aggregate and optimizing the parameters by a Levenberg-Marquardt algorithm. Several parameters introduce the mechanical properties contrast between the two phases: the initial dislocation density, the alloy friction, the average lath size, the obstacle size, while the other parameters are kept identical (see Table 5).

Volume

b (nm)

H (s -1 ) ∆G 0 (eV) A second optimization was conducted by considering the material as single-phase. Although less physically detailed, this optimization provide a reference material of same effective response in tension, allowing to compare stress field heterogeneities between the two models without changing the average mechanical response. This second list of parameters is provided in Table 6.

τ 0 (MPa) ρ m (mm -2 ) γ0 (s -
Figure 9 demonstrates the good agreement between the two models, single-phase and dual-phase, with the experimental tensile data at -150 • C, a relevant temperature to conduct the fracture study. The voronoi aggregate mesh used to compute the tensile curves are also shown on the same Figure. Average stress-strain curves per grain in the tensile direction are reported Figure 10 and show that stresses in bainitic grains are higher than in ferritic grains, as expected, with a large scatter in stress levels as soon as yielding occurs. At 10% macroscopic strain, the average deformation of ferritic phase is 11% higher than that of bainite.

failure behaviour of RPV steel

At low temperature a possibility of brittle failure by cleavage exists for 16MND5 steel [START_REF] Wallin | Statistical model for carbide induced brittle fracture in steel[END_REF][START_REF] Tanguy | Local approach to fracture based prediction of the delta-T56J and deltatK1c100 shifts due to irradiation for an A508 pressure vessel steel[END_REF]. Modeling the brittle failure is of industrial relevance and being able to predict its probability for irradiated material has been a key objective of European projects like Perfect and then Perform60 [START_REF] Massoud | PERFECT : prediction of irradiation damage effects on reactor components: A summary[END_REF][START_REF] Leclercq | Perform 60 -prediction of the effects of radiation for reactor pressurevessel and in-core materials using multi-scale modelling -60 years foreseen plant lifetime[END_REF] Energy absorbed during fracture for 16MND5 steel shows that the ductile to brittle transition takes place between -150 • C and 0 • C [START_REF] Tanguy | Modélisation de l'essai Charpy par l'approche locale de la rupture[END_REF]. Using the same single-phase and dual-phase 98 grain aggregate as introduced earlier to study the behaviour, additional crystal plasticity calculations were conducted at -196 • C and -150 • C since these two temperatures belong to the brittle failure domain. Loading conditions were triaxial with a chosen triaxiality of 2.5 (cf. Eq. 21) to represent the average stress field in front of a compact tension specimen crack tip. Practically, the loading was prescribed by imposing average loads to the 3 sets of adjacent faces by means of discrete linking elements and integral relationships, using the following relationship between average reactions in each direction:

αΣ 1 = Σ 2 = Σ 3 , α < 1 , χ = 1 + 2α 3 -3α (21)
The triaxiality was also computed after the calculation by averaging the triaxiality directly computed from the principal stresses, which validated the loading conditions. Since the conditions are imposed as average conditions, displacement fluctuations are allowed on each face and artificial stress concentrations are prevented. It is then possible to use a reasonably small representative volume as the elements near the boundary can be included in the post treated volume. The volume of the aggregate is set to (50µm) 3 .

As shown by Figure 11, the macroscopic stresses in the three principal directions slightly differ between the single-phase and the dual phase calculations, while no mismatch was found in uniaxial tension. This difference can be attributed to the increased stress heterogeneities induced by dual-phase model and more pronounced for triaxial loadings, resulting in lower yield stress compared to the single phase model.

Chosen local approach to failure model

The model used to predict failure is based on an earlier work in which a Beremin-type post treatment [START_REF] Beremin | A local criterion for cleavage fracture of a nuclear pressure vessel steel[END_REF] was applied to the stress and strain fields computed for a polycrystalline aggregate [START_REF] Mathieu | Analyse et modélisation micromécanique du comportement et de la rupture fragile de l'acier 16MND5 : prise en compte des hétérogénéités microstructurales[END_REF]. This failure model is based on the following 6 hypotheses:

• H1: The carbides and other second phase particles are the initiating locations of fracture. The effect of larger particles (MnS) is not taken into consideration, since they mostly affect the ductile failure properties [START_REF] Tanguy | Modélisation de l'essai Charpy par l'approche locale de la rupture[END_REF].

• H2: The plasticity is considered as mandatory (together with a critical stress) in the initiation process. This will lead to consider a threshold of plastic deformation in the initiation criterion.

• H3: The carbide population is composed of clusters in the middle of bainitic laths and of larger isolated carbides. All these carbides are potential candidates to provoke cleavage. The defect sizes follow a Frechet distribution [START_REF] Ortner | Characterisation of Euro-A reference steel for application of a model of brittle fracture[END_REF].

• H4: Once a crack initiates at a defect, its propagation in the neighbouring ferritic matrix is considered to be driven by a Griffith-style criterion. The ferrite cleavage principally takes place in the {001} family. The corresponding Griffith stress is given by:

σ Grif = πEγ p 2(1 -ν 2 )a , ( 22
)
where E corresponds to Young's modulus, a is the crack size and γ p , surface energy in J/m 2 term, is derived from a former optimization for 16MND5 [START_REF] Mathieu | Analyse et modélisation micromécanique du comportement et de la rupture fragile de l'acier 16MND5 : prise en compte des hétérogénéités microstructurales[END_REF]:

γ p = 2.15 + 1.53 • e 0.009 • T (23)
• H5: The largest defects (of dimension greater than r max c ) are initiated in the beginning of plasticity deformation and are blunted instead of being propagated. As a result, there exists a threshold in carbide size, beyond which the cleavage mechanism will not be activated. In the following this threshold is set to one micron.

• H6: The overcome of micro-structural barriers is not a critical mechanism, as long as the stress needed for the propagation of a defect is reached. The weakest chain hypothesis is applied here : the propagation of a micro-defect in the ferritic matrix leads to the fracture of the considered volume.

A key data to apply such micromechanical failure model is about the carbide distribution in material. The study of steel Euro Material A by Ortner et al. [START_REF] Ortner | Characterisation of Euro-A reference steel for application of a model of brittle fracture[END_REF] not only gives the size distribution of carbides, but also the carbide number per volume unit N:

p(r ≥ r 0 ) = f (r 0 ) = 1 -exp(-( r 0 3.6 • 10 -8 ) -2.7 ), N = (7.6 + 2.4) • 10 17 m -3 (24
)

n ave = N v, σ 0 = N v 2 • 10 11 (25)
A given carbide number obtained from a truncated normal distribution of average n ave and standard variation σ 0 [START_REF] Macdougall | Simulating Computer Systems[END_REF] is distributed to each gauss point associated to an elementary the volume v. In this manner some gauss points do not get carbides, as shown in Figure 13(a), and some others get multiple carbides. In the latter case only the largest carbide is kept for the calculation of critical stress σ Grif . Those gauss points without carbides are allocated with very high σ Grif . By this method a distribution of critical cleavage stress is obtained, that is a function of r max c and γ p . Most local approaches calculate the critical stress based on the first principal stress but another stress measure is chosen in this model. As the crystal orientation for each gauss point is known, the projection of stress tensor on the cleavage planes can be easily performed,

σ max c = M ax{σ p } (26)
This effective maximal cleavage stress is compared with the Griffith stress and the equivalent plastic deformation is compared with the critical value ε mini p . According to the model hypothesis H6, the fracture occurs as soon as the following criterion is verified:

σ max c ≥ σ Grif , ε equ p ≥ ε mini p (27)
In the finite element calculation, the simulation is performed step by step. Inside time steps all the variables are linearly interpolated.

The post treatment described above is repeated N s times, for each realization the spatial arrangement of carbides is re-computed. In order to obtain the failure probability distribution function, the whole results are sorted in increasing order of fracture time. A rank i=1,2... N s is allocated to each result and the corresponding fracture probability P R is given as :

P R = i -0.5 N s ( 28
)
Equation 24 is employed to simulate the distribution of carbides in both ferritic and bainitic phases. It should be noted that this is an approximation and that additional calculations can be carried out once the realistic distributions are obtained in each phase and eventually within the phases.

Results and discussion

The first calculation of cumulated failure probability has been carried out directly on the Voronoï aggregates. A result summary is provided in table 7 and the crack initiating locations (elements) are shown in figure 14, as well as the grain boundaries. The simulated failure curve at -150 • C and the separated responses of the two phases are shown in figure 15.

A conclusion can be drawn here that given the mechanism and hypothesis explicated above, while the crack initiation takes place only in ferrite at -196 • C, it occurs in both phases at -150 • C, especially at grain boundaries and model interfaces. This result is in with anticipation, as in comparison to -150 • C, the σ Grif is smaller at -196 • C, due to the smaller value of γ p . In addition, softer ferritic grains are prone to higher strain levels thus are more likely to pass the threshold of plastic deformation compared to bainitic grains. For each calculation, the last criterion to be reached (Griffith stress or plastic strain) is denoted as critial. For each temperature, this critical criterion is reported in The initiating elements for the 200 carbide samplings are shown in Figure 14, together with the ferrite-bainite interfaces. Initiation sites are often close to these interfaces, also correlated with stress concentrations.

Figure 15 shows, for the dual phase calculation at -150 • C, the repartition of failure initiation site between the two phases. About 80% of initiation events occur in bainitic phase, which is consistent with the Griffith stress limiting failure mechanism as reported in Table 7 Figure 16 shows that while the principal stresses at 50% failure possibility for single-phase and dual-phase models are equal, the failure probability distribution significantly differs in the low probability regime for the two models. The failure curve of dual-phase model spreads broader than the single-phase one, indicating a lower Beremin exponent. This is rather expected since increasing the microstructural and mechanical heterogeneties likely increases failure probability when considering a weakest link scenario. Naturally, this conclusion could be different depending on the considered carbide spatial distribution.

Figure 17 shows the cloud of carbide size and principal stress as crack initiates. It can be concluded that while relatively smaller (0.6µm to 0.8µm) carbides can lead to failure at high

Conclusion

In conclusion, a preliminary step of this work consisted in identifying plastic flow models of the RPV steels, for two distinct micromechanical models : one single phase, and one dual phase. This identification, especially in the dual-phase case, could be improved by considering more information about phase-to-phase heterogeneities like initial dislocation densities. Other improvement possibility would be to take advantage from in-situ SEM experiments together with digital image correlation to identify the mechanical contrast between the two phases.

Introducing the dual phase model was motivated by the phase distribution variation through the thickness. Coupling accurate microstructural description with plastic flow model and failure model would address the need of fracture toughness prediction at different depths of the component. Although these large scale heterogeneities were not yet predicted in this work, some preliminary assessments have indicated the interest of a dual phase model. Indeed, the dual phase model is associated with a lower macroscopic stress for lower failure probabilities.

This work also intends to illustrate the variety of computational micromechancals tools available that might help to understand or identify some root causes of failure, experimentally difficult or impossible to obtain. Naturally, the relevance of the obtained results strongly depends of the modeling hypothesis and the accuracy of mechanical fields computed at the microstructural scale. This in turn motivates the development of additional and complementary validation procedures for the different elementary bricks used in this approach : metallurgical quantitative information, plastic flow behaviour, failure behaviour.

• A 7 . 5

 75 hours tempering at 630 • C-645 • C • A final annealing lasting 8 hours at a temperature of 610 • C

Figure 1 :Figure 2 :

 12 Figure1: Carbon content in H1BQ12 steel versus the sampling position[START_REF] Bethmont | Evolution de la ténacité dans l'épaisseur des viroles de cuves des réacteurs à eau préssurisée[END_REF]

Figure 3 :

 3 Figure 3: Ultimate tensile stress of H1BQ12 steel versus the sampling position [24] Phase 1/2 depth 1/4 depth inner diameter s d = 140mm s d = 70mm s d = 0 -15mm Ferrite % 25 ± 3 16 ± 3 x Carbide clusters % 25 ± 2 25 ± 1 15 ± 2 Bainite % 50 54 77 Martensite % 0 5 ± 4 18 ± 2

1. 2 Figure 4 :

 24 Figure4: Stress-strain curves of 16MND5 samples of 3/4 depth, H1BQ12 material, reproduced from[START_REF] Renevey | Approches globale et locale de la rupture dans le domaine de la transition fragile-ductile d'un acier faiblement allié[END_REF]

Figure 5 :

 5 Figure 5: Yield stress sensitivity analysis results for some key model parameters

Figure 6 :

 6 Figure 6: Linear fitting of Σ 0.02 over D lat

  (b) showing the effect of K self on the hardening of pure iron, presented by the difference between Σ 0.1 and Σ 0.02 . Decreasing K self increases the hardening. As shown by Figure7(a), K self has a negligible influence on yield stress Both Figures7(a) and 7(b) prove the validity of the reference value of K self = 100 for pure iron.

Figure 7 :

 7 Figure 7: Influence of parameter K self .

Figure 8 :

 8 Figure 8: Realization example of a 162 grain polycrystal, Voronoi tesselation (a) and exploded mesh showing the distribution of phases with 75% bainite content (b)

Figure 9 :

 9 Figure 9: Macroscopic stress-strain curves at -150 • C compared with experimental values (a). The meshed aggregate used to compute the tensile curves is shown with colored grains (a) and colored phases for the 75% bainite fraction case (b).

Figure 10 :

 10 Figure 10: stress-strain curves for each grain for a model with 75% bainite with parameters of table 5

Figure 11 :

 11 Figure 11: Average stress-strain curves in the three tensile components obtained during the triaxial tensile calculation at -150 • C.

Figure 12 :Figure 13 :

 1213 Figure 12: Distribution of carbides

Figure 14 :

 14 Figure 14: Distribution of crack initiating elements (in blue) at -150 • Cfor the dual-phase model. The ferrite-bainite interface mesh is displayed in magenta.

Figure 15 :

 15 Figure 15: Failure curve calculated at -150 • C

Figure 16 :

 16 Figure 16: Failure curves of single-phase and dual-phase models at -150 • C

Figure 17 :

 17 Figure 17: Carbide size at crack initiation against principal stress

Table 4 :

 4

			N(elements)	N(grains)
	Whole model	100%	26718	98
	Ferrite	24.96% 7092 (26.5%)	24
	Bainite	75.04% 19626 (73.5%)	74

Basic information for dual-phase Voronoï model with 75% bainite

Table 5 :

 5 Parameters at -150 • C used in DD_CC law for dual-phase model

	b (nm)	H (s -1 )	∆G 0 (eV)	τ 0 (MPa) ρ m (mm -2 ) γ0 (s -1 )	n	y AT (nm)
	0.248	1 • 10 11	0.84	363	1.2 • 10 7	10 -6	20	2
	Phase D obs (nm) ρ ini (mm -2 ) τ F (MPa) D lath (mm) (a sj AT ) coli (a sj AT ) =coli	K self
	Ferrite	10	2.4 • 10 7	8	0.088	0.7	0.1	30
	Bainite	20	1.3 • 10 8	85.4	0.005	0.7	0.1	30

Table 6 :

 6 Parameters at -150 • C used in DD_CC law for single-phase model

	1 )	n	y AT (nm)

obs (nm) ρ ini (mm -2 ) τ F (MPa) D lath (mm)

Table 7 .

 7 At the lowest of -196 • C, the plastic deformation is more critical than the cleavage stress, and the trend is opposite at -150 • C.

		Failure initiation	Critical : σ Grif	Critical : ε mini
	Whole model	200	200	81	120	119	80
	Ferrite	200	39	81	39	119	0
	Bainite	0	161	0	81	0	80

p Temperature -196 • C -150 • C -196 • C -150 • C -196 • C -150 • C

Table 7 :

 7 Failure calculation information at -150 • C and -196 • C for dual-phase Voronoï model