Marc Bonnet
email: bonnet@athena.polytechnique.fr

Regularized direct and indirect symmetric variational BIE formulations for three-dimensional elasticity. 1

Regularized direct and indirect symmetric variational

BIE formulations for three-dimensional elasticity.

1 Introduction and motivation.

Although variational Galerkin-type boundary integral equation (GBIE) formulations have received a lot of attention in the past, mostly in the fielf of applied mathematics (see e.g. Bendali [START_REF] Bendali | Problème aux limites extérieur et intérieur pour le système de Maxwell en régime harmonique[END_REF], Cortey-Dymont [START_REF] Ph | Simulation numérique de problèmes de diffraction d'ondes par une fissure[END_REF], Nedelec [START_REF] Nedelec | Integral equations with non integrable kernels[END_REF], Wendland [START_REF] Wendland | Mathematical properties and asymptotic error estimates for elliptic boundary element methods[END_REF]), developments and applications in the field of solid mechanics (see Sirtori et al. [START_REF] Sirtori | A Galerkin symmetric boundary element method in elasticity: formulation and implementation[END_REF], Parreira and Guiggiani [START_REF] Parreira | On the implementation of the Galerkin approach in the boundary element method[END_REF], Becache [START_REF] Becache | Résolution par une méthode d'équations intégrales d'un problème de diffraction d'ondes élastiques transitoires par une fissure[END_REF], Nishimura and Kobayashi [START_REF] Nishimura | A boundary integral equation method for an inverse problem related to crack detection[END_REF], Kane and Balakrishna [START_REF] Kane | Symmetric Galerkin boundary formulations employing curved elements[END_REF]) have been few compared to the very popular collocation boundary integral equation (CBIE) methods.

Variational approaches do have attractive conceptual features compared to usual collocation methods: known convergence results are more extensive [START_REF] Nedelec | Integral equations with non integrable kernels[END_REF] and actually better than those of CBEM using the same discretization. They lead to symmetric stiffness-like BEM matrices, so that the solution CPU time, which dominates for large DOF numbers, is asymptotically half that necessary for solving the full nonsymmetric system produced by collocation methods. Indeed the numerical examples in Balakrishna et al. [START_REF] Balakrishna | Efficient analytical integration of symmetric Galerkin boundary integrals over curved elements: elasticity formulation[END_REF] show clearly that, for a given discretization, GBEA becomes computationally more efficient than CBIE for sufficiently large numbers of DOFs. GBIE formulations can be viewed as weighted residual formulations of CBIEs, and thus usually involve double surface integrals (see however Bui [START_REF] Bui | On the variational boundary integral equations in elastodynamics with the use of conjugate functions[END_REF], where a symmetric BIE formulation with a single surface integration is discussed). In order to get a symmetric GBIE formulation, weighted residual forms of displacement and traction CBIEs are used respectively on the Dirichlet and Neumann parts of the boundary. Although the collocation traction BIE is known to require a C 1,α continuity for the boundary displacement u at the collocation point, which causes some technical difficulties especially when using curved elements for 3-D situations, its Galerkin counterpart requires only a C 0,α continuity for u. The usual C 0 conformal shape functions can then be used without difficulty in a GBEA context.

On the other hand, the implementation of GBIE formulations for general engineering problems is somewhat more involved than that of displacement CBIE, which explains perhaps its relative lack of popularity. Also, investigation of GBIE formulations for mixed boundary value problems began only recently. Another important issue is the handling of singular integrals: the double boundary integrals that appear in Galerkin-type BIE formulations are convergent but some of the successive simple integrals are not. One possibility is a proper direct handling of the strongly singular and hypersingular simple surface integrals that appear, using either a numerical method (Guiggiani et al. [START_REF] Guiggiani | A general algorithm for the numerical solution of hypersingular boundary integral equations[END_REF]) or analytical integration of the singularity (Sirtori et al. [START_REF] Sirtori | A Galerkin symmetric boundary element method in elasticity: formulation and implementation[END_REF], Balakrishna et al. [START_REF] Balakrishna | Efficient analytical integration of symmetric Galerkin boundary integrals over curved elements: elasticity formulation[END_REF]; the latter used symbolic computation). Another possibility is to perform, as is done here, an analytical regularization before any discretization.

In the present paper, we investigate the statement of a GBIE formulation for the mixed boundary value problem of linear elastostatics, for either interior or exterior problems, obtained from the stationarity conditions associated to the minimization principle of an augmented potential energy. This variational BIE formulation appears to be precisely that obtained from a weighted residual viewpoint. Also, the present formulation is regularized in the spirit of previous works on indirect regularization of displacement CBIEs [START_REF] Bonnet | Méthode des équations intégrales régularisées en élastodynamique tridimensionnelle[END_REF], [START_REF] Bonnet | Regularized Boundary Integral Equations for Three-dimensional Bounded or Unbounded Elastic Bodies Containing Curved Cracks of Arbitrary Shape Under Dynamic Loading[END_REF], [START_REF] Bonnet | Regularization of the Displacement and Traction BIE for 3D Elastodynamics using indirect methods[END_REF] and of traction BIEs ( [START_REF] Nishimura | A boundary integral equation method for an inverse problem related to crack detection[END_REF]). The resulting variational BIE formulation is expressed in terms of weakly singular inner surface integrals followed by outer regular integrals. Also, it is a direct formulation, whose unknowns are the physical displacements and tractions on the boundary. As a by-product, the indirect variational BIE formulation for mixed problems, whose unknowns are fictitious densities, is also established from the direct one.

Finally, some implementation features are discussed.

Elastic problem, potential energy and stationarity conditions

Let Ω denote a three-dimensional elastic body, either bounded or unbounded; n denotes the unit normal vector outwards from Ω. The elastic displacement u(y) (y ∈ Ω) is governed by the homogeneous Navier equation:

C ijab u a,bj = 0 (1)
where C denotes the Hooke tensor of elastic constants. The comma indicates partial differentiation; Einstein summation convention is used throughout. Besides, mixed boundary conditions are considered:

   u = u D on S u T n (u) = t D on S T (2)
where S u , S T = ∂Ω -S u are complementary portions of ∂Ω and

T n (u) = [C : ε(u)].n
is the elastic traction vector associated with u. If Ω is unbounded, suitable decay conditions at infinity must be added to [START_REF] Bendali | Problème aux limites extérieur et intérieur pour le système de Maxwell en régime harmonique[END_REF].

The solution u to the elastic problem (1-2) minimizes the following augmented potential energy functional E(v):

E(v) = Ω ε(v) : C : ε(v) dV - S T t D .v dS - Su t.(v -u D ) dS (3)
The last integral term in (3) has the effect that u is an unconstrained minimizer of E(v):

no kinematical constraints are imposed on v. Thus u can be sought such that the first variation of E vanishes at v = u, i.e., letting v = u + δu, ve have: The next important step is the actual construction of test functions δu that satisfy [START_REF] Becache | Résolution par une méthode d'équations intégrales d'un problème de diffraction d'ondes élastiques transitoires par une fissure[END_REF].

δE(u).δu = Ω ε(u) : C : ε(δu) dV - S T t D .
In order to do so, recall that any such δu admits the integral representation (x interior to Ω):

δu k (x) = - ∂Ω δu i (y)n j (y)Σ k ij (x, y) dS y + ∂Ω δt i (y)U k i (x, y) dS y ( 6
)
where δt is the traction vector T n (δu). Similarly, any δu + which solves (1) in the exterior domain Ω + = R 3 -Ω satisfies at any x interior to Ω the complementary (with reference to Ω + ) integral representation:

0 = - ∂Ω δu + i (y)n j (y)Σ k ij (x, y) dS y + ∂Ω δt + i (y)U k i (x, y) dS y (7)
where δt + stands for T n (δu + ). The normal n outwards from Ω is used in both [START_REF] Bonnet | Regularization of the Displacement and Traction BIE for 3D Elastodynamics using indirect methods[END_REF] and [START_REF] Bui | On the variational boundary integral equations in elastodynamics with the use of conjugate functions[END_REF]. The symbols U k i (x, y), Σ k ij (x, y) denote the i-and ij-components of the elastic displacement and stress fields created at y by a unit point force applied at x along the k-direction (elastostatic fundamental solution). The full space (Kelvin) solution, a half space solution with a free-surface condition (Mindlin), or any other fundamental solution defined on a subset of R 3 that includes Ω may be used for this purpose. Any of them has the symmetry property:

U a k (x, y) = U k a (y, x) (8)
which in turn implies:

C ijab ∂ ∂x b U k a (x, y) = Σ a ij (y, x) C ijab ∂ ∂x b Σ a kl (x, y) = C k ab ∂ ∂y b Σ a ij (y, x) (9)
Next, ( 7) is subtracted from (6), giving:

δu k (x) = ∂Ω ũi (y)n j (y)Σ k ij (x, y) dS y - ∂Ω ti (y)U k i (x, y) dS y
using the displacement and traction jumps accross ∂Ω: ũ = δu + -δu and t = δt + -δt.

It is indeed legitimate to restrict the above formula by imposing δu + = δu, i.e. ũ = 0 on S u and δt + = -δt, i.e. t = 0 on S T , which gives:

δu k (x) = S T ũi (y)n j (y)Σ k ij (x, y) dS y - Su ti (y)U k i (x, y) dS y (10) = δu T k (x) + δu U k (x) (11)
Such test functions δu satisfy (1) and reflect the boundary condition structure of the initial mixed elastostatic problem under study. Note that, since ũ must be continuous all over ∂Ω, one has:

ũ | ∂ ST = 0 (12)
The stress tensor δσ = σ(δu) is then given from ( 10) by the representation:

δσ ij (x) = C klab S T ũk (y)n (y) ∂ ∂y b Σ a ij (y, x) dS y - Su tk (y)Σ k ij (y, x) dS y (13)
where use has been made of ( 8), ( 9).

Derivation of the direct symmetric bie formulation

Basically the symmetric BIE formulation comes from the substitution of ( 10), ( 13) into the stationarity equation ( 5). This leads to a formulation in terms of double surface integrals. However, integrability problems, related to the presence of kernels behaving like r -1 , r -2 , r -3 where r = |y -x| (especially the last two), require a limiting process of some kind.

Introduction of an auxiliary surface

Let S be a closed, regular surface exterior to Ω (hence exterior to ∂Ω for interior problems, but interior to ∂Ω for exterior problems), defined by means of a one-to-one mapping F onto ∂Ω:

y ∈ ∂Ω → z = F (y) ∈ S
which is left unspecified. This surface S is substituted to ∂Ω in ( 10), [START_REF] Krishnasamy | Hypersingular boundary integral equations: their occurrence, interpretation, regularization and computation[END_REF], following an approach used by Sirtori et al. [START_REF] Sirtori | A Galerkin symmetric boundary element method in elasticity: formulation and implementation[END_REF]. The idea is to perform some analytical manipulations, for regularization purposes, on the double surface integrals over ∂Ω × S (in which case the elastic kernels are nonsingular) and then consider the limiting process S → ∂Ω.

The test functions δu, δσ are thus given by: (with Su = F (S u ), ST = F (S T )). Recall that eqns. ( 14), ( 15) apply for Ω either bounded or unbounded, n being the unit normal outwards from Ω.

δu k (x) = ST ũi (y)n j (y)Σ k ij (x, y) dS y - Su ti (y)U k i (x, y) dS y ( 14
)
δσ ij (x) = C klab ST ũk (y)n (y) ∂ ∂y b Σ a ij (y, x) dS y - Su tk (y)Σ k ij (y, x) dS y (15)
Upon substitution of ( 14), ( 15), the stationarity condition ( 5) for E(u) can be split into two independent parts:

   ∀ũ δE T = 0 ∀ t δE U = 0 (16)
where δE T and δE U collect the terms containing ũ and t respectively that arise from the substitution of ( 14), ( 15) into (5).

Evaluation of δE U

From eqns ( 5), ( 14), [START_REF] Nedelec | Integral equations with non integrable kernels[END_REF], one has:

δE U = I 1 + I 2 (17)
I 1 = Su t k (x) Su ti (z)U k i (x, z) dS z dS x + S T t D k (x) Su ti (z)U k i (x, z) dS z dS x (18)
I 2 = - Su u D i (x)n j (x) Su tk (z)Σ k ij (z, x) dS z dS x - S T u i (x)n j (x) Su tk (z)Σ k ij (z, x) dS z dS x (19)
First, the term I 1 is made of a (potentially) weakly singular surface integral followed by a regular one and thus needs no regularization; it is left unchanged.

Then, I 2 , eqn. [START_REF] Polizzotto | An energy approach to the boundary element method. Part I: elastic solids[END_REF], is rewritten as:

I 2 = - Su tk (z) Su [u D i (x) -u D i (y)]n j (x)Σ k ij (z, x) dS x dS z - Su tk (z) S T [u i (x) -u D i (y)]n j (x)Σ k ij (z, x) dS x dS z - Su tk (z)u D i (y) ∂Ω n j (x)Σ k ij (z, x) dS x dS z ( 20
)
where y ∈ S u is the point such that F (y) = z ∈ Su (hence u(y) = u D (y)). Integrations over ∂Ω and S have been interchanged, this operation being valid because of the current nonsingular character of the integrals.

Besides, one notes that the fundamental stress tensor satisfies, by virtue of equilibrium, the following identity (the source point z being exterior to Ω):

∂Ω n j (x)Σ k ij (z, x) dS x = κδ ik (21)
where κ = 0 (Ω bounded, n exterior to ∂Ω) or κ = 1 (R 3 -Ω bounded, n interior to ∂Ω).

Thus, the last integral in [START_REF] Rizzo | A boundary integral equation method for radiation and scattering of elastic waves in three-dimensions[END_REF] becomes:

Su tk (z)u D i (y) ∂Ω n j (x)Σ k ij (z, x) dS x dS z = κ Su tk (z)u D k (y) dS z (22)
These manipulations result in a regularization of the initially strongly singular integral w.r.t. x in [START_REF] Polizzotto | An energy approach to the boundary element method. Part I: elastic solids[END_REF]. Indeed one notes that, since

U (z, x) ∼ |z -x| -1 Σ(z, x) ∼ |z -x| -2
the expressions ( 18), [START_REF] Rizzo | A boundary integral equation method for radiation and scattering of elastic waves in three-dimensions[END_REF] involve weakly singular integrals w.r.t. x followed by nonsingular integrals w.r.t. z, provided the regularity requirement u ∈ C 0,α (∂Ω), i.e.:

∃C > 0, ∃α ∈]0, 1] |u(y) -u(x)| ≤ C |y -x| α (23)
is met. Then, the limiting expression of I 2 , I 1 for S → ∂Ω is obtained by a mere substi-tution of ( Su , z) by (S u , y) in eqns. [START_REF] Parreira | On the implementation of the Galerkin approach in the boundary element method[END_REF], [START_REF] Rizzo | A boundary integral equation method for radiation and scattering of elastic waves in three-dimensions[END_REF]. The result is thus:

I 1 = Su t k (x) Su ti (y)U k i (x, y) dS y dS x + S T t D k (x) Su ti (y)U k i (x, y) dS y dS x (24)
I 2 = - Su tk (y) Su [u D i (x) -u D i (y)]n j (x)Σ k ij (y, x) dS x dS y - Su tk (y) S T [u i (x) -u D i (y)]n j (x)Σ k ij (y, x) dS x dS y -κ Su tk (z)u D k (y) dS y (25)
For eqn. (25) above, use has been made of the substitution of ( 21) into [START_REF] Polizzotto | An energy approach to the boundary element method. Part I: elastic solids[END_REF].

Evaluation of δE T

From eqns ( 5), ( 14), ( 15), one has:

δE T = J 1 + J 2 (26)
J 1 = - Su t k (x) ST ũi (z)n j (z)Σ k ij (x, z) dS z dS x - S T t D k (x) ST ũi (z)n j (z)Σ k ij (x, z) dS z dS x (27) J 2 = Su u D i (x)n j (x)C klab ST ũk (z)n (z) ∂ ∂z b Σ a ij (z, x) dS z dS x + S T u i (x)n j (x)C klab ST ũk (z)n (z) ∂ ∂z b Σ a ij (z, x) dS z dS x (28)
First, J 1 is subjected to a treatment similar to that made to I 2 . For this purpose, eqn.

(27) is rewritten, upon addition and subtraction of ũ(x) to ũ(z), as:

J 1 = - Su t k (x) ST ũi (z)n j (z)Σ k ij (x, z) dS z dS x - S T t D k (x) ST [ũ i (z) -ũi (x)]n j (z)Σ k ij (x, z) dS z dS x - S T t D k (x)ũ i (x) ST n j (z)Σ k ij (x, z) dS z dS x ( 29
)
where the fact that ũ = 0 for x ∈ S u has been used. Besides, the fundamental stress tensor satisfies the following identity, similar to [START_REF] Sirtori | A Galerkin symmetric boundary element method in elasticity: formulation and implementation[END_REF] but with the source point x now being interior to Ω:

S n j (z)Σ k ij (x, z) dS z = (κ -1)δ ik ( 30
)
where again κ = 0 (Ω bounded, n exterior to ∂Ω) or κ = 1 (R 3 -Ω bounded, n interior to ∂Ω). The last integral in (29) then becomes:

S T t D k (x)ũ i (x) Su n j (z)Σ k ij (x, z) dS z dS x = S T t D k (x)ũ i (x) (κ -1)δ ik - Su n j (z)Σ k ij (x, z) dS z dS x = (κ -1) S T t D k (x)ũ k (x) dS x - S T t D k (x)ũ i (x) Su n j (z)Σ k ij (x, z) dS z dS x (31)
Substitution of (31) into (29) allows for an expression of J 1 which is made, in the limiting case S → ∂Ω, of weakly singular integrals over S T followed by nonsingular integrals over S u , S T , as follows:

J 1 = - Su t k (x) S T ũi (y)n j (y)Σ k ij (x, y) dS y dS x - S T t D k (x) S T [ũ i (y) -ũi (x)]n j (y)Σ k ij (x, y) dS y dS x + (1 -κ) S T t D k (x)ũ k (x) dS x + S T t k (x)ũ i (x) Su n j (y)Σ k ij (x, y) dS y dS x (32)
The weakly singular character of the inner integral in the above equation relies upon ũ being C 0,α -continuous over ∂Ω, eqn. (23). The property (12) of ũ, together with (23), implies:

∃C > 0, ∃α ∈]0, 1] |u(y)| ≤ C |y -x| x ∈ S u , y ∈ S T
and thus also contributes to the weakly singular character of (32):

Next, let us consider the second integral J 2 in (26), which involves the hypersingular

kernel C klab ∂ ∂z b Σ a ij (z, x)
. Following an argument given by Nishimura and Kobayashi [START_REF] Nishimura | A regularized boundary integral equation method for elastodynamic crack problems[END_REF] and using the symmetry property [START_REF] Guiggiani | A general algorithm for the numerical solution of hypersingular boundary integral equations[END_REF], one notes that, for x = y:

∂ ∂x j C klab ∂ ∂x b Σ a ij (z, x) = C klab ∂ ∂x b ∂ ∂x j Σ a ij (z, x) = 0 ∂ ∂z l C klab ∂ ∂z b Σ a ij (z, x) = C ijab ∂ ∂x b ∂ ∂z l Σ a kl (x, z) = 0
This implies the existence of a fourth-order tensor A such that:

C klab ∂ ∂z b Σ a ij (z, x) = e iep e jf q e kgr e lhs ∂ ∂x e ∂ ∂x f ∂ ∂z g ∂ ∂z h A pqrs (z, x) (z = x) (33)
where e ijk are the components of the permutation tensor. For instance, the tensor A(z, x) associated with the three-dimensional isotropic Kelvin fundamental solution is given (Nedelec [START_REF] Nedelec | Integral equations with non integrable kernels[END_REF]) by:

A pqrs (z, x) = µ 8π δ pr δ qs + δ ps δ qr + 2ν 1 -ν δ pq δ rs |z -x| (34)
The decomposition (33) allows for the use of the following variant of Stokes' formula:

S ge abc n a f ,b dS = - S f e abc n a g ,b dS
where S is any regular surface; moreover either S is closed or f = 0 on ∂S. Here it is applied twice (once w.r.t. z and once w.r.t. x) to J 2 , eqn. (28); this gives:

J 2 = Su (Ru D ) iq (x) ST (Rũ) ks (z)B ikqs (z, x) dS z dS x + S T (Ru) iq (x) ST (Rũ) ks (z)B ikqs (z, x) dS z dS x (35)
using the notations

(Ru) iq = e jf q n j u i,f (36)

Statement of the direct symmetric BIE formulation

The previous results are now put together. From ( 16), the definitions ( 17), (26) and the results (24), ( 25), ( 32), (38), the direct symmetric BIE formulation for the mixed elastic boundary value problem reads:

Find u ∈ V u , t ∈ V T    B uu (u, ũ) + B tu (t, ũ) = L u (ũ) B ut (u, t) + B tt (t, t) = L t ( t) ∀ũ ∈ V 0 u , ∀ t ∈ V T (39)
using the bilinear forms

B uu (u, ũ) = S T (Ru) iq (x) S T
B iqks (y, x)(Rũ) ks (y) dS y dS x (40)

B tu (t, ũ) = - Su S T t k (x)ũ i (y)n j (y)Σ k ij (x, y) dS y dS x (41) B ut (u, t) = - Su S T tk (y)[u i (x) -u D i (y)]n j (x)Σ k ij (y, x) dS x dS y (42) B tt (u, t) = Su Su t k (x) ti (y)U k i (x, y) dS y dS x (43)
and the linear forms

L u (ũ) = (κ -1) S T t D k (x)ũ k (x) dS x + S T t D k (x) S T [ũ i (y) -ũi (x)]n j (y)Σ k ij (x, y) dS y dS x - S T t D k (x)ũ i (x) Su n j (y)Σ k ij (x, y) dS y dS x - S T (Ru D ) iq (x)
Su B iqks (y, x)(Rũ) ks (y) dS y dS x (44)

L t ( t) = κ Su tk (x)u D k (x) dS x + Su tk (y) Su [u D i (x) -u D i (y)]n j (x)Σ k ij (y, x) dS x dS y - Su S T t D k (x) ti (y)U k i (x, y) dS y dS x (45)
and where the spaces of admissible functions V u , V 0 u , V T are defined as:

V u = u | u ∈ C 0,α (S T ) and u = u D on ∂S T V 0 u = ũ | ũ ∈ C 0,α (S T ) and ũ = 0 on ∂S T V T = t | t piecewise continuous on S T
Note that, contrary to (u, ũ), no restriction is made on the values of t, t along ∂S T .

Eqns. (39) to (45) constitute the main result of this paper. Being the stationarity condition (5) applied to an appropriate family of test functions δu, they define a variational BIE formulation for the mixed elastostatic boundary-value problem. In addition, the formulation is symmetric by virtue of the relations

∀(u, ũ) ∀(t, t)          B uu (u, ũ) = B uu (ũ, u) B tt (t, t) = B tt ( t, t) B ut (u, t) = B tu ( t, u) (46)
which are visible consequences of formulae ( 40), ( 41), ( 42), (43). Those symmetry properties may be viewed as consequences of the symmetry properties ( 8), [START_REF] Guiggiani | A general algorithm for the numerical solution of hypersingular boundary integral equations[END_REF], which are true for any elastostatic fundamental solution.

Moreover, the present variational symmetric BIE formulation is direct, i.e. expressed in terms of the physical elastic boundary unknowns u | S T , t | Su . The prescribed displacements u D and tractions t D appear in the linear forms L u and L t .

Also, the formulation (39 to 45) is valid whatever fundamental solution is used. However, the kernel B(y, x) depends on the actual fundamental solution used and is currently known analytically only in a few cases like the Kelvin solution (see eqn. (34) for the isotropic case or Becache [START_REF] Becache | Résolution par une méthode d'équations intégrales d'un problème de diffraction d'ondes élastiques transitoires par une fissure[END_REF] for the anisotropic case).

The present formulation is established for elastic domains either bounded or unbounded. They slightly differ from each other by the presence of the coefficient κ in the expressions (44), (45) for the linear forms L u , L t .

Finally, the regularization of the J 2 term, subsection 4.3, was first proposed by Nedelec [START_REF] Nedelec | Integral equations with non integrable kernels[END_REF], for Neumann problems. It is valid provided u, ũ are C 0,α -continuous over ∂Ω and S T respectively. On the other hand, the weighted residual formulation of the traction CBIE requires a C 1,α continuity for u at the collocation point no matter which integration technique is used. This is of practical importance since C 0,α interpolations are easier to handle and more flexible, especially for 3-D situations.

Relation to weighted residual statements

The displacement and traction exterior representation formulas are respectively:

0 = ∂Ω u i (x)n j (x)Σ k ij (z, x) dS x - ∂Ω t i (y)U k i (z, x) dS x (47) 0 = n j (z) C ijab ∂Ω u k (x)n (x) ∂ ∂z b Σ a k (z, x) dS x - ∂Ω t k (y)Σ k ij (z, x) dS x ( 48
)
where z is any point exterior to Ω and n(z) any unit vector attached to z. One then readily sees that equations (17-24-25) are obtained by taking the inner product of eqn.

(47) by a fictitious traction ti (z), then integrating the result over ST . Also, (26-32-38) come from the inner product of eqn. (48) by a fictitious displacement ũi (z), integrated over Su , with n the unit normal to ST pointing away from Ω. The surface S = Su ∪ ST is defined as in subsection 4.1 above, and the symmetry properties ( 8), ( 9) are used.

As a result, the variational symmetric BIE formulation (39 to 45) is identical to that obtained from a weighted residual statement of displacement and traction collocation BIEs and a regularization approach following the lines of subsections 4.2, 4.3. It is also similar in nature to the weighted residual statement developed in [START_REF] Kane | Symmetric Galerkin boundary formulations employing curved elements[END_REF], from which it differs mainly by the kind of regularization used for the originally hypersingular kernel.

Formulations for Dirichlet amd Neumann problems

The formulation (39 to 45) is applicable to the special cases of Dirichlet (S u = ∂Ω, S T = ∅)

or Neumann (S T = ∂Ω, S u = ∅) boundary conditions. The resulting variational BIEs are as follows, for the Dirichlet elastic problem

Find t ∈ V T B tt (t, t) = L t ( t) ∀ t ∈ V T (49)
with

B tt (u, t) = ∂Ω ∂Ω t k (x) ti (y)U k i (x, y) dS y dS x L t ( t) = κ ∂Ω tk (x)u D k (x) dS x + ∂Ω tk (y) ∂Ω [u D i (x) -u D i (y)]n j (x)Σ k ij (y, x) dS x dS y
and the Neumann elastic problem:

Find u ∈ V u B uu (u, ũ) = L u (ũ) ∀ũ ∈ V 0 u (50)
with

B uu (u, ũ) = ∂Ω (Ru) iq (x) ∂Ω B iqks (x, y)(Rũ) ks (y) dS y dS x L u (ũ) = (κ -1) ∂Ω t D k (x)ũ k (x) dS x + ∂Ω t D k (x) ∂Ω [ũ i (y) -ũi (x)]n j (y)Σ k ij (x, y) dS y dS x

Plane elasticity

The result (39) can be specialized to plane elasticity. In this case, formulae (41) to (45) remain valid, using the plane strain/stress counterpart of the 3D fundamental solution (the range of all indices is {1, 2}).

Moreover, one readily sees that in eqn. (33) indices i, j, k, l and e, f, g, h should also range in {1, 2}, which implies p = q = r = s = 3. Thus:

C klab ∂ ∂z b Σ a ij (z, x) = e 3ep e 3f q e 3gr e 3hs ∂ ∂x e ∂ ∂x f ∂ ∂z g ∂ ∂z h A 3333 (z, x) = e ep e f q e gr e hs ∂ ∂x e ∂ ∂x f ∂ ∂z g ∂ ∂z h A(z, x) (z = x) (51)
for either plane strain or plane stress (with e 12 = -e 21 -1, e 11 = e 22 = 0). Eqn. (40)

then reduces to:

B uu (u, ũ) = S T S T du i ds (x) dũ k ds (y)e ie e kg ∂ 2 ∂x e ∂y g A(x, y) ds y ds x (52)
and a similar expression should be used in the relevant term of L u , eqn. (44).

For the two-dimensional isotropic plane strain Kelvin solution, one has (Bonnemay [START_REF] Bonnemay | Equations intégrales pour l'élasticité plane[END_REF]):

A(z, x) = µ 4π(1 -ν) r 2 ln r r = |z -x| (53)
5 Indirect variational BIE formulation

The direct formulation (39 to 45), being valid for both bounded or unbounded domains, can be used to derive an indirect formulation as well. Consider a bounded regular domain Ω together with its complement Ω = R 3 -Ω and denote by u I , t I and u E , t E respectively the solutions to the coupled interior and exterior elastostatic problems with equal prescribed displacements and opposite prescribed tractions:

u I = u E = u D on S u t I = -t E = t D on S T
Let also n and n denote the opposite normals to ∂Ω, respectively outwards from Ω and Ω.

Then, the stationarity condition ( 16) is considered for the coupled problem. The test functions δu should be equal for both the interior and exterior problem, thus ( 14) and [START_REF] Nedelec | Integral equations with non integrable kernels[END_REF] should ultimately be taken with equal ũ, t and the same unit normal, say n(y). This extends to eqns. ( 24), ( 25), (32), (38). On the other hand, the unit normal n(x) in the latter equations is the one that appears in T n (δu) in [START_REF] Bonnet | Regularized Boundary Integral Equations for Three-dimensional Bounded or Unbounded Elastic Bodies Containing Curved Cracks of Arbitrary Shape Under Dynamic Loading[END_REF] and thus should change sign.

Therefore, the variational formulation (39) becomes, for the coupled problem:

Find u E , u I ∈ V u , t E , t I ∈ V T    -B uu (u E , ũ) + B uu (u I , ũ) + B tu (t E , ũ) + B tu (t I , ũ) = -L u (ũ) -B ut (u E , t) + B ut (u I , t) + B tt (t E , t) + B tt (t I , t) = L t ( t) ∀ũ ∈ V 0 u , ∀ t ∈ V T
where B uu , etc. are given by ( 40), ( 41), ( 42), (43) with the normal n outwards from Ω, while from ( 44), (45) many terms cancel out when appropriately combined, so that:

L u (ũ) = S T t D i (x)ũ i (x) dS x L t ( t) = Su u D i (x) ti (x) dS x
Let ∆u = u Eu I and ∆t = -t Et I denote respectively the displacement and traction jumps accross ∂Ω (note that t I and t E are defined in terms of the opposite normals n and n, and also that ∆u should be continuous over ∂Ω and satisfy [START_REF] Balakrishna | Efficient analytical integration of symmetric Galerkin boundary integrals over curved elements: elasticity formulation[END_REF]). Thus, due to the linearity of B uu , etc. w.r.t. each argument, the variational formulation for the coupled problem takes the following, clearly symmetric, form:

Find ∆u ∈ V 0 u , ∆t ∈ V T    B uu (∆u, ũ) + B tu (∆t, ũ) = L u (ũ) B ut (∆u, t) + B tt (∆t, t) = -L t ( t) ∀ũ ∈ V 0 u , ∀ t ∈ V T (54)
One notices that the complexity of the right-hand sides is dramatically reduced between the direct formulation (39) and the indirect one (54), which means substantial savings in terms of computational time. However, it should be emphasized that, in order to recover the physical quantities u on S u and t on S u , one still has to compute the appropriate integral representations, the limiting expression of:

u k (z) = S T ũi (y)n j (y)Σ k ij (z, y) dS y - Su ti (y)U k i (z, y) dS y σ ij (z ) = C klab S T ũk (y)n (y) ∂ ∂y b Σ a ij (y, z ) dS y - Su tk (y)Σ k ij (y, z ) dS y
for z → x ∈ S T , z → x ∈ S u , both leading to strongly singular integrals. Note that their computation thus relies on ∆u, ∆t being both C 0,α at the point x of interest, while the assumption made on t for the direct formulation (39) is weaker. Also, there may be some cases for which the boundary elastic quantity themselves are only of secondary interest (e.g. in some fracture or optimization problems), giving a good reason for using the indirect formulation (54).

Another comment of interest is that, in the indirect case, the value E at equilibrium (i.e. for ∆u, ∆t solution to ( 54)) of the elastic potential energy is given (e.g. from (3))

by:

E = 1 2 [L t (∆t) -L u (∆u)]
Thus the Formulation (54) applied to the special case of the Dirichlet elastic problem gives:

Find ∆t ∈ V T , ∀ t ∈ V T ∂Ω ∂Ω ψ k (x) ti (y)U k i (x, y) dS y dS x = ∂Ω tk (x)u D k (x) dS x
while one has for the Neumann elastic problem:

Find ∆u ∈ V 0 u , ∀ t ∈ V 0 u ∂Ω (Rφ) iq (x) ∂Ω B iqks (x, y)(Rũ) ks (y) dS y dS x = ∂Ω ũk (x)t D k (x) dS x
The two above formulations may be found e.g. in Nedelec [START_REF] Nedelec | Integral equations with non integrable kernels[END_REF].

Numerical implementation

The symmetry properties (46) allow a reduction of computational effort for the linear system solution step due to the symmetric character of the governing matrix; this has been extensively discussed elsewhere. We wish to focus the discussion on the treatment of the singular surface integrals in (40) to (45), i.e. repeated integrals over the same element E, for the simplest case of isoparametric conformal BEM interpolation. The surface ∂Ω is divided into boundary elements, which are mapped on a reference element E 0 , (usually

either the square ξ = (ξ 1 , ξ 2 ) ∈ [-1, 1] 2 or the triangle 0 ≤ ξ 1 + ξ 2 ≤ 1)
. Let the discussion be restricted to the consideration of a single element E. The location of a point y on E is expressed in terms of N shape functions N K and N geometrical nodes y K (K = 1 . . . n) located on the boundary of E:

y = N K=1 N K (ξ)y K (ξ ∈ E 0 ) (55)
Then, for (α, β) = 1, 2, the natural basis (a α ), metric tensor (g αβ ) and unit normal n on E are given by:

a α (ξ) = N K=1 N K ,α (ξ)y K g αβ (ξ) = a α (ξ) • a β (ξ) J(ξ)n(ξ) = a 1 ∧ a 2 J 2 (ξ) = (g 11 g 22 -g 2 12 )(ξ) (ξ ∈ E 0 ) (56)
The shape functions N K are also used to interpolate u, t, ũ, t. Moreover, it can be shown that:

∀ξ ∈ E 0 e abc n b (y) ∂f ∂y c (y)J(ξ)e a = ∂f ∂ξ 1 a 2 - ∂f ∂ξ 2 a 1 (y = y(ξ)) (57)
Let η = (η 1 , η 2 ) denote the antecedent of x on E 0 . Following a common practice in BEM (see e.g. Rizzo et al. [START_REF] Rizzo | A boundary integral equation method for radiation and scattering of elastic waves in three-dimensions[END_REF]), set ξ 1 = η 1 + ρ cos α, ξ 2 = η 2 + ρ sin α. Then:

dS y = J(ξ)dξ 1 dξ 2 = J(ξ)ρdρdα (58)
Define N q (ρ, α; η), r(ρ, α; η), Û k i (ρ, α; η), Σk is (ρ, α; η), Biqks (ρ, α; η) as follows:

N K (ξ) -N K (η) = ρ N K (ρ, α; η) r(x, y) = K K=1 N K (ξ)y K = ρr(ρ, α; η) U k i (x, y) = 1 ρ Û k i (ρ, α; η) Σ k is (x, y) = 1 ρ 2
Σk is (ρ, α; η) [4], allows their numerical evaluation to be performed with standard product Gaussian quadrature formulas.

The double surface integrals over E × E are computable by using standard Gaussian quadrature (points η i , weights w i ) for the outer one:

I sing e,e = E 0 E 0 f (η, ξ) dη dξ ≈ n i=1 w i E 0 f (η i , ξ) dξ
The remaining, weakly singular, integral is evaluated using formulas (60), (61), (62).

Figure 1 :

 1 Figure 1: Boundary ∂Ω and auxiliary surface S

B iqks (x, y) = 1 ρ

 1 Biqks (ρ, α; η) (59) where r(ρ, α; η) = 0 for sufficiently small ρ and Û (ρ, α; η), Σ(ρ, α; η), Biqks (ρ, α; η) are thus regular at ρ = 0. Hence typical (inner) singular surface integrals in (40) to(45) become:E n s (y)Σ k is (x, y) (u i (y) -u i (x)) dS y is (ρ, α; η) M K (ρ, α; η)J(ξ)dρdα (60) E B iqks (x, y)e lhs n l (y)u k,h (y) dS y , α; η) (a β .e s )M K ,α -(a α .e s )M K ,β dρdα (61) E U k i (x, y)t i (y) dS y (ρ, α; η)M K (ξ)J(ξ)dρdα(62)Equations (60), (61), (62) take full advantage of the regularization. A further coordinate change (ρ, α) → (v 1 , v 2 ), in order to convert them to integrals over the square [-1, 1] 2

7

  Concluding comments.In this paper, symmetric variational regularized BIE formulations are established, for the mixed elastostatic boundary-value problem. A direct version (in terms of unknown boundary displacements and tractions) is first presented. Then an indirect version (in terms of unknown fictitious densities) is in turn established. Both are expressed using at most weakly singular integrals followed by regular integrals. These formulations provide a basis for the numerical solution of general 3D situations. Conventional boundary element interpolations of any degree for the geometry as well as the unknowns can be implemented in a straightforward manner, using conventional techniques for numerical integration together with a now well-established singularity cancellation method at the shape function level.Moreover, the variational BIE formulations are shown to express the stationarity of an augmented potential energy functional, thus being truly a variational BIE formulation.The present approach can be developed to other situations, such as potential problems or elastodynamic problems, as well. Moreover, similar variational BIE formulations are already known in the literature, notably for Dirichlet or Neumann problems in various contexts, including crack problems (see references quoted herein).

  δu dS -

				t.δu dS = 0	(4)
				Su
	with no kinematical constraint on δu. Indeed, from
	ε(u) : C : ε(δu) dV =	T n (u).δu dS -	div [C : ε(u)].δu dV
	Ω	∂Ω		Ω
	it is readily shown that u solves (1-2) and that the Lagrange multiplier t equals the
	traction vector T n (u) on S u .
	Now let us suppose that the minimization of E(v) is attempted only for those v which
	satisfy the local equilibrium equation (1). Since this is also satisfied by the solution u,
	one is then led to restrict (4) to the trial functions δu = v -u which themselves satisfy
	(1). Upon integration by parts of the domain integral, the stationarity condition (4) thus
	takes the form
	δE(u).δu =	u D .T n (δu) dS +	u.T n (δu) dS
	Su		S T
	-	t D .δu dS -	t.δu dS = 0	(5)
		S T	Su
	which constitutes the starting point for the derivation of a direct variational BIE formu-
	lation.


Meriem Ben

Haj Slama

Vincent Taupin

Nabila Maloufi

Kaustubh Venkatraman

Anthony D Rollett

Ricardo A Lebensohn

Stephane Berbenni

Benoît Beausir

Antoine Guitton
email: antoine.guitton@univ-lorraine.fr

Stéphane Berbenni

Electron Channeling Contrast Imaging characterization and crystal plasticity modelling of dislocation activity in Ti21S BCC material

Keywords:

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

Introduction

In metals, elementary mechanisms such as dislocation glide, operate at the atomic and electronic scale; typically, the Ångström and the picosecond. They can strongly influence the macroscopic strength at the laboratory scale, e.g. the centimeter and the second. For example, abinitio simulations of dislocation cores in Body Centered Cubic (BCC) metals highlighted the occurrence of non-Schmid effects, which could explain asymmetry between tension and compression at the macroscopic scale [START_REF] Ismail-Beigi | Ab initio study of screw dislocations in Mo and Ta: a new picture of plasticity in bcc transition metals[END_REF][START_REF] Rodney | Ab initio modeling of dislocation core properties in metals and semiconductors[END_REF][START_REF] Lia | Dislocation-core symmetry and slip planes in tungsten alloys: Ab initio calculations and microcantilever bending experiments[END_REF] .

To understand and ultimately control the mechanical behavior of metallic materials with regards to their processed microstructures, one has to face this impressively complex multiscale problem that can span more than 10 orders of magnitude in time and space. This is why advanced theoretical and experimental tools are constantly being developed, aiming at bridging relevant length scales and at understanding the mechanisms involved.

In support of these aims, experimental methods include several imaging (High-Resolution Transmission HR-TEM, Scanning Tunneling STM, Atomic Force Microscopy AFM) and diffraction techniques (High-Energy tion densities inside the grains and their interactions with grain boundaries. Modeling tools include Discrete Dislocation Dynamics (DDD), Crystal Plasticity (CP) models based on Finite Elements or Fast Fourier Transform algorithms (CPFEM, CPFFT), or mean-field approaches (e.g., self-consistent models such as the ViscoPlastic Self-Consistent (VPSC), and Elasto-Visco-Plastic Self-Consistent (EVPSC) schemes) [START_REF] Bertin | Frontiers in the simulation of dislocations[END_REF] . Recent studies have been carried out that aim to couple models to experiments such as CPFEM (or CPFFT) with EBSD [START_REF] Hansen | An investigation of geometrically necessary dislocations and back stress in large grained tantalum via EBSD and CPFEM[END_REF][START_REF] Zhang | Korsunsky , Datasets for multi-scale diffraction analysis (synchrotron XRD and EBSD) of twinning-detwinning during tensile-compressive deformation of AZ31B magnesium alloy samples[END_REF][START_REF] Pi | Modeling uniaxial tensile deformation of polycrystalline Al using CPFEM[END_REF] , or CPFEM with High-Resolution Digital Image Correlation HRDIC [START_REF] Di Gioacchino | An experimental study of the polycrystalline plasticity of austenitic stainless steel[END_REF] . Such coupled approaches can provide information on textures, Geometrically Necessary Dislocation (GND) density patterning, plastic strain localization and transfer across grain boundaries, etc. Among the results, such coupled modeling-experimental methods show that it is possible to link the microstructural evolution to stress-strain behavior.

However, the majority of models revealed that each configuration usually has its own behavior, especially concerning dislocation interactions, such that a huge number of simulations are required to characterize all possible mechanisms. Perhaps more importantly, the analysis of individual "model " configurations, by means of molecular dynamics for example [ 13 , 14 ] although very instructive, usually do not or cannot consider the non-local influence of the surrounding microstructure on the local mechanisms investigated. For instance, the activation of a particular dislocation system in a grain or the interaction of a dislocation with a grain boundary are obviously strongly dependent on the elastic fields induced by nearby grains, grain boundaries, triple lines and other possible defects such as solutes or precipitates.

A few studies have tried to make an overview of constitutive laws and theories for crystal plasticity modelling [ 15 , 16 ], but to our knowledge, none of them provided quantitative and statistical comparisons between these several modelling methods that could shed new light on the inherent mechanisms. In this paper, a combined experimental/theoretical approach is proposed to fill this gap.

Concerning our experimental approach, the strategy reported in this paper is based on Scanning Electron Microscope (SEM) measurements. SEM is flexible essentially because of the large sample surface area that can be imaged, as compared to TEM. Moreover, in SEM, the nondestructive Electron Channeling Contrast Imaging ECCI imaging mode provides TEM-like diffraction contrast imaging of near-surface defects with a good accuracy [START_REF] Kriaa | Fundamental and experimental aspects of diffraction for characterizing dislocations by electron channeling contrast imaging in scanning electron microscope[END_REF][START_REF] Mansour | Accurate electron channeling contrast analysis of a low angle sub-grain boundary[END_REF][START_REF] Guitton | A dislocation-scale characterization of the evolution of deformation microstructures around nanoindentation imprints in a TiAl Alloy[END_REF][START_REF] Mansour | Accurate electron channeling contrast analysis of dislocations in fine grained bulk materials[END_REF][START_REF] Guyon | Sub-micron resolution selected area electron channeling patterns[END_REF] . In terms of length scales, centimetric bulk materials can be characterized, enabling statistical studies. A more powerful technique, Accurate Electron Channeling Contrast Imaging (A-ECCI) developed some years ago [START_REF] Guyon | Sub-micron resolution selected area electron channeling patterns[END_REF] , has an angular accuracy better than 0.1°with a spatial resolution finer than 500 nm [START_REF] Sachs | Zur Ableitung einer Fliessbedingung[END_REF] . By coupling A-ECCI to an in situ tensile test inside the SEM [START_REF] Ben | In situ macroscopic tensile testing in SEM and electron channeling contrast imaging: pencil glide evidenced in a bulk 𝛽-Ti21S polycrystal[END_REF] , several areas on a bulk sample can be accessed with a single tool while under mechanical load, thus covering the scale from sample to dislocations. The availability of multiple regions of interest is helpful in enumerating relevant statistical parameters from the mesoscopic and microscopic descriptions of microstructures, which in turn provides large enough datasets for data mining.

For numerical calculations, we consider three mean-field models and a full-field CPFFT model. The first mean field model uses the Sachs assumption, i.e., that the local stress state is equal to the macroscopic one [START_REF] Reuss | Berechnung der Fliessgrenze von Mischkristallen auf Grund der Plastizitätsbedingung für Einkristalle[END_REF] . This assumption is commonly made in experimental characterization of dislocation activity, where uniform stress is applicable in contrast to the uniform strain assumed in the Taylor model [START_REF] Taylor | Plastic strains in metals[END_REF] , which is more successful for texture prediction at large strain. Hence, our hypothesis is that the uniform stress model is more reliable than others, at low plastic strains, for a material with many possible slip systems. In addition, we use two popular models, i.e. , the VPSC model [ 26 , 27 ], a well-known model for prediction of texture development in polycrystals, and a more recent affine-type EVPSC model based on the translated field method by Mareau and Berbenni [START_REF] Mareau | An affine formulation for the self-consistent modeling of elasto-viscoplastic heterogeneous materials based on the translated field method[END_REF] . The affine-type EVPSC model was recently applied to study the effect of anisotropic elasticity stiffness components (C 11 , C 12 , C 44 ) on the macroscopic behavior of polycrystalline 𝛽-Ti alloys, especially the transition between elastic and plastic states, and the development of incompatibility stresses ( i.e., inter-granular internal stresses) within the polycrystal [START_REF] Lhadi | Micromechanical modeling of the effect of elastic and plastic anisotropies on the mechanical behavior of beta-Ti alloys[END_REF] . The fourth full-field model (CPFFT) is an elasto-viscoplastic formulation with anisotropic elasticity/plasticity [ 30 , 31 ]. The different model predictions are compared with experimental data in terms of slip system activities at the onset of plasticity. For the majority of metals, slip activation at yield follows the well-known Schmid law [START_REF] Christian | Some surprising features of the plastic deformation of Body-Centered Cubic metals and alloys[END_REF][START_REF] Taylor | Thermally-activated deformation of BCC metals and alloys[END_REF][START_REF] Duesbery | Plastic anisotropy in BCC transition metals[END_REF] and the expected results are that: ( i ) the slip system with the highest Resolved Shear Stress (RSS) is activated (before others) and ( ii ) the RSS of the activated slip system is necessarily equal to or higher than its associated Critical Resolved Shear stress (CRSS).

Marked by an atypical and still not well understood plastic behavior, BCC metals are reported to clearly violate this law [START_REF] Duesbery | Plastic anisotropy in BCC transition metals[END_REF] . Theoretical calculations showed that the CRSS of the three BCC slip families depend on the crystal orientation with respect to the loading direction [START_REF] Orlans-Joliet | Plane strain compression of silicon-iron single crystals[END_REF] . More recent work proposed a modified parameter-free Schmid law to quantify the deviations from the Schmid law based on the deviations of the screw dislocation trajectory away from a straight path between equilibrium configurations [START_REF] Dezerald | Plastic anisotropy and dislocation trajectory in BCC metals[END_REF] . This deviation from the Schmid law is metal-dependent [START_REF] Aono | Orientation dependence of slip in niobium single crystals at 4.2 and 77 K[END_REF] . Furthermore, the number of activated slip systems in BCC metals can be as high as 48, and "pencil glide" (see §.3.1) is frequently reported in experiments. For those reasons, we choose a BCC metal as a challenging benchmark for our proposed coupled experimental / theoretical strategy, i.e. , 𝛽-Ti21S. Note further that elastic anisotropy is important in this material [START_REF] Wang | An atomistically-informed dislocation dynamics model for the plastic anisotropy and tension-compression asymmetry of BCC metals[END_REF] , making the analysis more complex, even at the onset of plasticity.

In what follows, we present the development of an experimentaltheoretical approach that allows comparing several modeling methods while contrasting their results with experiment.

The experimental work that we describe below is driven by the following objectives: ( i ) perform observations on a bulk specimen, ( ii ) provide an in situ observation of dislocations and slip lines during deformation, ( iii ) characterize a large polycrystalline area where hundreds of plastic events (dislocations and slip lines) can be captured, ( iv ) conduct an automated statistical analysis of the data collected.

Modeling with several crystal plasticity theories enables ( i ), documentation and quantification of the long-range influence of internal stresses due to neighboring grains and defects on the local mechanisms observed and ( ii ), the use of statistical analysis to improve mechanical models.

Comparison with simulations is assisted by an automated analysis of the observed microstructures so as to provide a sufficiently rich set of plastic events, in view of quantifying statistically the errors made by the models. It contributes as well to bridging the gap between elementary mechanisms at the dislocation scale and the overall behavior of the polycrystalline microstructure.

Experiment and modeling

Material and experimental setup

As reported in our previous work [START_REF] Ben | In situ macroscopic tensile testing in SEM and electron channeling contrast imaging: pencil glide evidenced in a bulk 𝛽-Ti21S polycrystal[END_REF] detailing the experimental procedure:

Bulk tensile sample ( Fig. 1 a) was cut from a 1.78 mm thickness rolled sheet of 𝛽-Ti21S alloy, produced by Titanium Metals Corporation / Toronto, Canada. The chemical composition of the material is Ti-15.9Mo-2.7Nb-2.9Al-0.2Fe-0.2Si (wt%). As a heat treatment, a soaking at 843 °C for 14 min was first applied before air cooling. The obtained microstructure is reported to be fully 𝛽 phase [START_REF] Weiss | Mechanical properties and microstructural study of homogeneous and heterogeneous laser welds in 𝛼, 𝛽, and 𝛼 + 𝛽 titanium alloys[END_REF] . Tensile test was performed at room temperature with a DEBEN (Suffolk, UK) machine of a maximum load cell of 1 kN. The sample was deformed under an imposed strain rate of 3.3 × 10 -4 s -1 . Note that the unknown machine stiffness was not removed from the stressstrain curve. Therefore, the slope of the linear part ( Fig. 1 d) does not correspond to the Young modulus of the material. The stress field was confirmed to be uniaxial and homogeneous along the operational length of the sample. The sample was mechanically polished with 1 μm diamond paste, followed by chemo-mechanical polishing with colloidal silica suspension. Finally, 2 h ion-polishing in a PECS II (GATAN, USA) machine was applied with a 3 keV beam. Then, detailed in situ characterizations of the microstructure during deformation, were performed by A-ECCI for enabling numerous observations of defects on several and large areas of interest (Dislocations and Slip-lines), in a Zeiss Auriga Scanning Electron Microscope (SEM, Oberkochen, Germany) operating at 10 kV.

Since the BSE yield strongly depends on the orientation of the crystal relative to the incident electron beam ( i.e., optical axis of the SEM), the orientation was accurately determined through High-Resolution Selected Area Channeling Patterns HR-SACP. The best defect contrast is obtained when the incident beam is oriented near a pseudo-Kikuchi band edge where the BSE yield is minimal. This corresponds to a channeling condition associated with the darkest contrast of the grain [START_REF] Morin | Scanning electron microscope observation of single defects in solid crystalline materials[END_REF] . All observations were performed at constant strain.

Prior to this step, i.e ., before deformation, EBSD mapping provided grain orientations to an accuracy of about 2°in a Zeiss Supra 40 SEM (Oberkochen, Germany) operating at 20 kV.

Automated Hough transform based detection of slip line crystallography

In order to establish an automated procedure for slip line crystallography, a trace analysis has been combined with an image analysis algorithm for slip line detection. Several experimental approaches have been widely applied for decades to study the deformation of polycrystals with a statistical significance. Yet, comparing our procedure to the existing methods was not the current work's purpose.

In the present study, the automation process is desirable not only for extracting useful statistics associated with the slip process but eventually for employing micrograph processing procedures for appropriate model selection. The ECC images were converted into a binary format with appropriate thresholding using Image-J [START_REF] Schneider | NIH image to imagej: 25 years of image analysis[END_REF] , to highlight the slip lines on a given region of interest. The thresholding procedure was conducted manually to highlight the slip lines in a given region of interest and to eliminate all other sources of contrast (noise). The thresholding procedure was carried out such that maximum number of lines were contrasted in the image with minimum noise. The binary image was used in a MATLAB [START_REF] Benker | MATLAB. In: Ingenieurmathematik mit Computeralgebra-Systemen[END_REF] program for automated slip line detection. The edges corresponding to each line and their corresponding pixels (x,y) are detected using a 'Canny edge detection operator' available in MATLAB. When viewed in Hough parameter space, points which are collinear in the Cartesian image space become readily apparent as they yield curves with a common intersection point. The standard Hough transform function available in MATLAB is used on the binary image to highlight the detected slip lines. Each detected line is highlighted in green with red and blue end points and tagged by its pixel co-ordinates and a line number ( Fig. 2 ). The angle ' 𝛼 1 ' for each line (computed with respect to the loading direction) along with its pixel co-ordinates and label number are stored as an output, based on analysis of the Hough transform.

With the EBSD map of the investigated zone all the theoretical slip plane traces corresponding to each plane family ({110}, {112}, {123}) were determined using the ATEX software [START_REF] Beausir | Université de Lorraine -metz, ATEX© -analysis tools for electron and X-ray diffraction[END_REF] . Their angles ( 𝛼) with respect to the loading direction were computed at each pixel location on the EBSD microstructure and stored separately as an output file. For all the lines that were detected, using the Hough transform algorithm, the angular difference (| 𝛼 1 -𝛼 2 |) between the experimental slip line angle and the theoretical angle of a given slip plane family is used to determine the most consistent slip plane family. A correlation coefficient was set up taking into account the ambiguity in determination.

The correlation coefficient for each detected line, based on the angular difference can be defined as follows:

⃗ R = ⎛ ⎜ ⎜ ⎝ r 1 r 2 r 3 ⎞ ⎟ ⎟ ⎠ = ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ cos [ min | | α 1 -α 2 | |{ 110 } ] cos [ min | | α 1 -α 2 | |{ 112 } ] cos [ min | | α 1 -α 2 | |{ 123 } ] ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ With |α1 -α2 | ∈ [ 0 , π 2 ]
. For some of the detected lines, there can be an ambiguity in determining the most consistent slip plane family {110}, {112} or {123}. To account for such ambiguities, a modified correlation vector is defined as follows:

⎛ ⎜ ⎜ ⎝ R 1 R 2 R 3 ⎞ ⎟ ⎟ ⎠ = | | | | | | | ⎛ ⎜ ⎜ ⎝ r 1 r 2 r 3 ⎞ ⎟ ⎟ ⎠ - ⎛ ⎜ ⎜ ⎝ max ( r 1 , r 2 , r 3 ) max ( r 1 , r 2 , r 3 ) max ( r 1 , r 2 , r 3 ) ⎞ ⎟ ⎟ ⎠ | | | | | | |
This will return the best matching 𝑅 𝑖 ( 𝑖 = 1 , 2 , 3 ) as zero and others to a non-zero value. A threshold can now be defined on the R i 's, with R i < 0.001, corresponding to an angular difference of less than 3°. The value corresponding to R i < 'threshold' is set to zero, thus classifying that particular line, 'i' as belonging to either two of the three possible slip systems.

According to the given convention, a detected slip line can belong to one of the six possible cases:

Based on this, an algorithm is formulated for counting the number of slip lines belonging to each slip system family and each detected slip line is assigned one of the six possible vectors:

The total number of slip traces corresponding to each slip plane family and the probability of observing a slip line corresponding to a given slip plane family P[N(R i )] is defined as:

⎛ ⎜ ⎜ ⎝ N ( R 1 ) N ( R 2 ) N ( R 3 ) ⎞ ⎟ ⎟ ⎠ = k=M ∑ k=1 ⎛ ⎜ ⎜ ⎝ n ( R 1 ) n ( R 2 ) n ( R 3 ) ⎞ ⎟ ⎟ ⎠ k
for M slip lines detected in a given region of interest.

P [ N ( R i )] for i =1 , 2 , 3 = N ( R 𝑖 ) N ( R 1 ) + N ( R 2 ) + N ( R 3 )
Probabilities are shown in pie charts ( Fig. 5 ). The uncertainties are estimated using Cochran statistics [START_REF] William | Experimental Designs[END_REF] . This is generally used for calculating the minimum sampling size ( 𝑛 𝛼 ) for a representative statistic:

n α = Z 2 C p α q α m 2 α
Where 𝑍 𝐶 is the selected critical value of desired confidence level, 𝐶, obtained from the 𝑍 𝐶 -score tables; 𝑝 𝛼 is the proportion of the glide plane family 𝛼 according to the random distribution ( 𝑝 { 123 } = 0.5; 𝑝 { 110 } = 𝑝 { 112 } = 0.25); 𝑞 𝛼 = 1 -𝑝 𝛼 and 𝑚 𝛼 is the desired margin of error for the glide plane family 𝛼.

In our study, 𝑁= 123 plastic events were tracked. Therefore, the sampling will be relevant statistically if 𝑁 ≥ 𝑛 𝛼 .

The Therefore, this estimation of uncertainties has to be understood as follows: one is sure at 95 % that the true percentage of { 123 } is localized between [ 48 % -9 %; 48 % + 9 % ] = [ 39 %; 57 % ] , where 48 % is the probability taken from Fig. 5 .

Brief description of the polycrystal plasticity models

As far as the simulation part is concerned, we performed calculations according to several numerical models to predict plastic activity (activated slip systems), with direct input from EBSD map of the same deformed zone.

For comparison purposes, the paper focused on the four following polycrystal plasticity models using different assumptions, which are here summarized:

1 The mean-field model using the Sachs hypothesis (uniform stress assumption) [START_REF] Reuss | Berechnung der Fliessgrenze von Mischkristallen auf Grund der Plastizitätsbedingung für Einkristalle[END_REF] . The local 3D Cauchy stress tensor is simply set equal to the macroscopic tensile stress (only one non-zero component for uniaxial tension). With this model, each point in the microstructures can be considered independently of all others and elastic/plastic interactions between grains are neglected. More importantly, only the effect of each grain crystallographic orientation through the Schmid tensor has an effect on the projection of the stress state onto each slip systems. Note that, as aforementioned, the Sachs hypothesis is commonly used in microscopy when studying dislocation activity. 2 The Visco-Plastic Self-Consistent (VPSC) model with a "tangent linearization " of the viscoplastic flow rule [START_REF] Molinari | A self-consistent approach of the large deformation polycrystal viscoplasticity[END_REF][START_REF] Lebensohn | A self-consistent approach for the simulation of plastic deformation and texture development of polycrystals: application to Zirconium alloys[END_REF][START_REF] Molinari | Tuning a self-consistent viscoplastic model by finite element results -I. Modeling[END_REF] . It has been used extensively in the literature. It predicts well in particular the evolution of crystallographic textures in polycrystals at large strains. Although this model would normally be compared with experimental data in the plastic domain (ignoring the elastic part of deformation), we applied it at the onset of plasticity. 3 The Elasto-Visco-Plastic Self-Consistent (EVPSC) mean field homogenization scheme solved by the "translated field " (TF) method with an "affine linearization " of the viscoplastic flow rule [START_REF] Mareau | An affine formulation for the self-consistent modeling of elasto-viscoplastic heterogeneous materials based on the translated field method[END_REF] . Here, the constitutive equations are the same as in the recent contribution of Lhadi et al. [START_REF] Lhadi | Micromechanical modeling of the effect of elastic and plastic anisotropies on the mechanical behavior of beta-Ti alloys[END_REF] , where it was applied to a fully 𝛽-Ti alloy. Here, we used the EVPSC model for the present 𝛽-Ti21S alloy, considering spherical grains with crystallographic orientations extracted from EBSD data, and compared its predictions with the other two mean-field models. In particular, elastic anisotropy is included in the EVPSC model, such that internal stresses due to both elastic and plastic anisotropies are considered. Therefore, it is expected to improve predictions with respect to the two previous mean field models. 4 Full-field polycrystal plasticity FFT (CPFFT) [START_REF] Ra Lebensohn | An elastoviscoplastic formulation based on fast Fourier transforms for the prediction of micromechanical fields in polycrystalline materials[END_REF] . This elastoviscoplastic model should give the best predictions because it includes elastic anisotropy and it solves locally the balance of Cauchy stress tensor at each grid point in the loaded microstructure. Although all models are applied to the microstructure digitized from the EBSD data, this model is the best at capturing the long-range elastic strains and thus stresses due to grain interactions in the microstructure. The simulations were performed using the EBSD maps as direct input. Accordingly, they corresponded to 2D simulations with no buffer zone and strictly periodic boundary conditions. Therefore, the obtained microstructure was constituted of columnar grains ( i.e. , elongated along the normal direction).

Comparison procedure

Using EBSD data of the same zone of interest (the zone is shown in Fig. 3 and contains 110 grains), polycrystalline plasticity calculations were done as follows:

-For EVPSC and CPFFT models: the strain hardening parameters were adjusted to fit the experimental tensile curve and used the same elastic parameters. Local stress tensors were extracted for each grain point and were used as MATLAB input data to finally map the RSS and the associated slip systems. The mapping procedure is described below. -For Sachs model, the local stress tensor was taken for each grain point equal to the macroscopic stress associated with uniaxial tensile test (only 𝜎 11 ≠0). -For VPSC model, all the steps from simulation to mapping were integrated in an in-house code, ATEX [START_REF] Beausir | Université de Lorraine -metz, ATEX© -analysis tools for electron and X-ray diffraction[END_REF] .

The local stress tensors were obtained in the sample coordinate system. With ATEX and MATLAB codes and as illustrated in Fig. 3 , they were rotated in the crystal coordinate system then multiplied by the Schmid tensor to get the RSS associated with slip systems. Then, using MATLAB and ATEX codes, 1 st maximum RSS, 2 nd maximum RSS… until n th maximum RSS maps could be generated, as well as the corresponding slip system maps. We chose to consider the n maximum RSS values that are less than 10 % difference from the 1 st maximum RSS value. This percentage corresponds to the threshold value above which all models give 100 % correct predictions for all active slip planes. In other terms, if we consider the n th maximum RSS different from the 1 st maximum RSS by more than 10 % (of the latter value), the four models can predict all the experimental slip planes. But below 10 %, the models start to be unable to predict some of the slip events.

To better illustrate this criterion, an "ambiguity map" was generated with each maximum RSS map (starting from the 2 nd ). The pixels associated to an RSS value higher than 10 % of the 1 st maximum RSS value are colored in black, the rest of pixels are colored in white.

For the studied material 𝛽-Ti21S there are 48 possible slip systems belonging to three slip plane families: {110} (12 slip systems), {112} (12 slip systems) and {123} (24 slip systems) [START_REF] Lewis | Slip systems and initiation of plasticity in a body-centered-cubic titanium alloy[END_REF] . This motivated showing, for each maximum RSS map, three slip system maps corresponding to each of the three plane families. Fig. 3. Numerical calculation path from EBSD map to slip system map. (Technically speaking, "R " is a matrix constructed from the Eulers that performs transformation of axes from sample to crystal, not rotation).

Results

Slip trace analysis

Using ECC micrographs obtained during the in situ tensile test, a slip trace analysis was performed. The analysis was done at the onset of plasticity (curve in Fig. 4 a) to determine the elementary mechanisms involved in the first steps of the deformation. After yielding and a few percent strain, the resulting microstructure is slightly deformed. Such a microstructure is investigated using A-ECCI to maximize the contrast from slip lines and dislocations [START_REF] Ben | In situ macroscopic tensile testing in SEM and electron channeling contrast imaging: pencil glide evidenced in a bulk 𝛽-Ti21S polycrystal[END_REF] .

The studied area from the tensile sample was chosen initially to cover slightly more than 100 grains, to allow a statistical study and enough data for micromechanical simulations using a statistically relevant region based on this sufficiently large number of grains. For detailed experimental observations and characterizations, ten particular regions of interest (each with a few grains) were selected. The examples reported in this communication correspond to a representative microstructure of the whole sample, allowing to explain the procedure and to discuss the main results.

Fig. 4 b shows the example region, loaded to a maximum stress of 810 MPa, that contains seven grains of interest and 55 plastic events (dislocations and slip lines). The slip trace analysis associates a slip plane type, {110}, {112} or {123}, to each event observed. This is done using the local orientation from EBSD data, based on which one can superimpose all possible slip systems on the slip plane trace. More precisely, the angles between the slip/dislocation lines and the loading direction ( i.e., with the horizontal) were first measured. Then, using the ATEX software [START_REF] Beausir | Université de Lorraine -metz, ATEX© -analysis tools for electron and X-ray diffraction[END_REF] and local orientation from the EBSD map of the investigated region as inputs, all the theoretical slip plane traces corresponding to each plane family were determined, as well as their angles to the loading direction. The smallest angle between a theoretical trace and the experimental one selects the most consistent slip plane family.

According to much previous literature, BCC plastic activity is governed by a pencil glide mechanism that does not involve a well-defined glide plane [START_REF] Taylor | The distortion of iron crystals[END_REF] . That is because in BCC crystals, {110}, {112} and {123} planes have comparable atomic packing densities and all contain < 111 > directions. Therefore, < 111 > -dislocations can slip into each of these glide planes [START_REF] Seeger | Anomalous slip -a feature of high-purity body-centred cubic metals[END_REF] . It turns out that dislocations can easily change their slip planes, i.e ., cross-slip during deformation. The dislocations move from one slip plane to another forming a stair-like shape, which results in a curved appearance in the slip trace [ 23 , 52 ].

Fig. 4 b shows the example for one curved line in grain G1 [1;1'] exhibiting pencil glide. The theoretical {110} ( Fig. 4 c), {112} and {123} plane traces are given by the ATEX software.

For instance, line [l] forms an angle closer to 𝛾°with respect to the horizontal. Calculating the differences with the real values of angles gives the closest expected slip plane family (smallest difference). Taking into consideration the angular resolution given by EBSD (2°) and the uncertainty of measuring the experimental angle (~1°), the slip plane family is determined at 3°of precision.

To automate the above slip trace analysis, the algorithm described in Section 2.2 was run on a total of 123 detected slip events and the probability distribution for the same is shown in Fig. 5 . The results show that the majority of observed slip traces belong to the {123} family. Besides, one can note that the statistical distribution between the three slip plane families {100},112} and {123} (almost 25 %/25 %/50 %), roughly reflects their respective multiplicity (12/12/24).

Statistical comparison of polycrystal plasticity models with experimental data

In order to evaluate the capability of each polycrystal plasticity model in predicting the experimentally observed slip lines, we followed an approach based on Resolved Shear Stress (RSS) analysis. We tested the predictions of the four models introduced in Section 2.3 , with expected increasing fraction of successful predictions: Sachs [START_REF] Reuss | Berechnung der Fliessgrenze von Mischkristallen auf Grund der Plastizitätsbedingung für Einkristalle[END_REF] , VPSC [START_REF] Molinari | A self-consistent approach of the large deformation polycrystal viscoplasticity[END_REF] , EVPSC [START_REF] Mareau | An affine formulation for the self-consistent modeling of elasto-viscoplastic heterogeneous materials based on the translated field method[END_REF] , and CPFFT [START_REF] Molinari | Tuning a self-consistent viscoplastic model by finite element results -I. Modeling[END_REF] .

Our comparative method allowed first, to statistically quantify the difference in the deduced parameters between the four models. Fig. 6 b shows an example of quantification of the maximum RSS (MRSS) obtained with the Sachs and the CPFFT models (the model hypotheses are detailed in Section 2.3 ). The local differences are not negligible; they are the highest especially near grain boundaries at which the physical properties are discontinuous, and they range up to 100 MPa.

Then, a comparison of the RSS values obtained by each model at the beginning of the plastic domain was performed in order to identify ( i ) which potentially active slip system corresponds to the maximum RSS and ( ii ) if the latter is higher than the CRSS i.e., if the Schmid law is locally respected. To do so, after determining the experimentally observed slip plane families for the studied 55 plastic lines, a mapping of the RSS associated with each slip system (48 slip systems for BCC crystals) was performed. Then, we compared the activated slip systems predicted by each model (according to the mapping procedure described in Section 2.4 ) with the experimental results. Fig. 6 c presents an example of comparison between experiments and CPFFT predictions. The first map (showing the center grain in green) is obtained by considering the 1 st MRSS and the second map (showing the middle grain in orange) is obtained by considering the 2 nd MRSS. Plane traces for lines 1 and 1' were experimentally determined to belong to {110} and {123} respectively. CPFFT calculations gave a 1 st MRSS value corresponding to {110} slip plane and a 2 nd MRSS value corresponding to {123} slip plane. The 2 nd MRSS value is retained since it was not very different from the 1 st MRSS (less than 10 %). Thus, we can reasonably say that the model predicts the appropriate slip plane family activation for both slip events.

The results obtained by applying the comparison procedure for the four models, show that the CPFFT and affine EVPSC give -as expectedbetter prediction than the Sachs (uniform stress) and the VPSC meanfield models. Yet, the difference in the well-predicted event percentage was quite small (86 % up to 90 % i.e., error percentage decreased from 14 % to 10 %); the statistics are summarized in Table 1 .

These four different models give different incompatibility stresses with different mechanical interactions between grains. Among the four models, the Sachs model does not include internal stresses: only the macroscopic stress is considered and the incompatibility stresses are zeros. In the VPSC model, only plastic incompatibilities are included: it does not describe elastic incompatibilities owing to material's anisotropic elasticity. Both the CPFFT and the EVPSC models include elastic and plastic incompatibilities in their formulations. Therefore, they give the best predictions among the different models about plastic activities in the elastic-plastic transition.

Comparison between MRSS and CRSS was made for the two models: EVPSC and CPFFT. For the first model, they were almost identical for all three slip plane families. For the second model, the CRSS was higher for {112} and {123} families (close values for both), than for the {110}. For both models, CRSS did not vary much compared to the initial value which is as expected because essentially no plastic slip was accumulated. Moreover, all MRSS corresponding to the activated slip planes were higher than the associated CRSS. From a statistical point of view, and for small strains, the Schmid law seems to be satisfied for this BCC material.

The same statistics as the ones reported in Table 1 were obtained for the analysis of shear rates and also when considering MRSS/CRSS values. Consequently, all interpretations and discussions made on RSS analysis can be drawn for shear rates and MRSS/CRSS fraction analysis too.

Discussion

Due to the SEM techniques used (2D EBSD and ECCI), all the reported results correspond to surfaces.

For relevant comparison, all the experimental observations as well as the numerical calculations were performed at the mesoscale.

For the studied material, 𝛽-Ti21S alloy, {123} and {112} slip-planes form 76 % of the experimentally observed traces, against only 24 % for {110}. This result is unexpected. The literature reports mainly {110} slip-planes for experiments with single crystals. Some of our {123} observations could be related to the surrounding microstructure in a polycrystal. In this context, Castany et al. [START_REF] Castany | In situ TEM study of dislocation slip in a metastable 𝛽 titanium alloy[END_REF] reported dislocation glide on {123} during an in situ TEM study of polycrystalline Ti alloy. In other cases, the observation of {112} and {123} can as well be explained by the fact that what we observed were "apparent " slip-planes, since the observations were made at meso-scale. These "apparent " {112} and {123} slip-planes may, at a finer scale, be sequences of short {110} segments. In fact, an atomic-scale investigation by Douat et al. [ 54 , 23 ], showed that microscopic {112} slip traces in Nb BCC structure are composed of periodic sequences of {110} traces detected by STM. In situ TEM investigations have also highlighted that dislocations were in some cases cross-slipping between different {110} close in average to a (123) plane [START_REF] Mompiou | Conventional vs harmonic-structured 𝛽-Ti-25Nb-25Zr alloys: a comparative study of deformation mechanisms[END_REF] . Even in the latter case, the experimental and theoretical results can still be compared since they were obtained both at the mesoscale.

For small strains, the Schmid law appears to be statistically satisfied. However, the Schmid law was not applicable locally for several cases, and especially not near grain boundaries, but in average it was applicable in a high fraction of cases. Yet, the tension/compression asymmetry of the yield stress is a good indicator of the presence of non-Schmid effects; this aspect, however, was not investigated here.

Concerning the critical resolved shear stress (CRSS) values, the EVPSC and CPFFT models gave different results; for the first model, CRSS was almost identical for all three slip plane families, while for the second, the CRSS was higher for {112} and {123} families. In the literature, contradictory results have been reported concerning CRSS values in BCC structures. {112} slip itself was reported to have a higher CRSS in the anti-twinning direction than in the twinning direction [ 56 , 57 ]. Compared with {110} slip, {112} slip was reported to exhibit higher CRSS in the anti-twinning direction but similar one for slip in the twinning direction [START_REF] Franciosi | Investigation of slip system activity in iron at room temperature by SEM and AFM in-situ tensile and compression tests of iron single crystals[END_REF] . Other authors stated that the CRSS of {112} slip plane family in general is lower than that of the {110} one [START_REF] Orlans-Joliet | Plane strain compression of silicon-iron single crystals[END_REF] . Here, it raises the question as to which key parameters control the CRSS values of each slip family in BCC crystals.

Although the different models compared above are based on different assumptions and take different microstructural parameters into account (e.g. grains are considered spherical in EVPSC calculations, whereas grains are elongated in the normal direction in CPFFT calculations), they provided good predictions at the onset of plasticity. Nevertheless, for more deformed states, lower efficiency of these models is expected, which can be confirmed (or not) by applying the approach described in this paper. Moreover, to be more rigorous, a reliable model must not only consider the interactions of the crystals with their neighboring grains [START_REF] Abdolvand | Strong grain neighbour effects in polycrystals[END_REF] , but also must include additional hypothesis for material behavior particularly for grain boundary-adjacent areas (different from those used for the rest of microstructure i.e., inside the grains). The latter is mandatory for the study of materials with very small grain sizes, and more complex microstructures. The procedure reported here can be applied to test how the activation of slip systems depends on the values of elastic parameters (C 11 , C 12 and C 44 ) measured and proposed in the literature for 𝛽-Ti21S or near-𝛽 Ti alloys [START_REF] Lhadi | Micromechanical modeling of the effect of elastic and plastic anisotropies on the mechanical behavior of beta-Ti alloys[END_REF] , and on CRSS values that can also be better fitted under higher plastic strains.

Conclusion

In this paper, we performed a multiscale and statistical analysis of plastic activity in BCC 𝛽-Ti21S, particularly at the onset of plasticity, coupling: in situ dislocation-scale characterization of plastic deformation mechanisms during macroscopic tensile test on bulk specimen, and numerical predictions of four polycrystal plasticity models based on mean field models: Sachs, VPSC, EVPSC and full-field simulations: CPFFT.

The coupling of in situ tensile testing with ECCI microstructural characterization represents a successful combination to avoid the limitations related to non-reliable statistics. It has the advantage of allowing analysis of many zones of interest and covers a maximum area of study in one single testing with a good accuracy. The relevance and reliability of our experimental results allowed a novel comparison with four different polycrystal plasticity models. This innovative procedure opens the way to validate and calibrate mesoscale models with large and statistically relevant datasets on more complex materials.

The study reported here can also be applied to in-depth statistical understanding of dislocation transfer mechanisms across grain bound-aries. Moreover, the described method provides input to data mining processes, and can be easily adopted by the material science community since the observation of slip lines is accessible in all SEMs equipped with Back-Scattered Electron (BSE) detectors. This approach can be coupled to DIC for a deeper investigation of slip step height, for instance.

Fig. 1 .

 1 Fig. 1. Experimental setups (a) tensile specimen form and dimensions (b) specimen mounted in the mini-tensile machine (c) mounting inside the microscope (d) applied tensile curve; reported A-ECCI observations were performed at the red cross point of the curve (after 15 min relaxation).

Fig. 2 .

 2 Fig. 2. Binary image used for the slip line detection algorithm.

Fig. 4 .

 4 Fig. 4. Determination of the experimental plane trace family. (a) Stress/strain tensile curve with loading/unloading parts (b) ECCI micrograph of the analyzed zone showing the slip trace segments [1] and [1'] from grain G1. (c) Theoretical slip plane traces of {110} associated to G1 crystallographic orientation.

Fig. 5 .

 5 Fig. 5. Pie chart showing the probability distribution of possible slip plane families using the slip line detection algorithm for the total of 123 slip events.

Fig. 6 .

 6 Fig. 6. Example of comparison between experimental observations and modeling. (a) ECC image of the large area of interest (b) the map Δ𝑀𝑅𝑆𝑆 = 𝑀 𝑅𝑆 𝑆 𝐹 𝐹 𝑇 -𝐸𝑉 𝑃 -𝑀 𝑅𝑆 𝑆 𝑆𝑎𝑐ℎ𝑠 of the same zone (c) ECC image of grain G1 and slip plane maps according to the associated color code: in green grains where {110} slip families were activated, in fuchsia grains where {112} slip families were activated and in orange grains where {123} slip families were activated.

  table below shows the best compromise for the Cochran parameters considering our sampling sizes:

	Parameters	Values
	𝑁	123
	𝐶	95 %
		{ 123 }	{ 110 }; { 112 }
	𝑝 𝛼	0.5	0.25
	𝑚 𝛼	9 %	8 %
	𝑛 𝛼	96	72

Table 1

 1 Statistics of predicted slip events in agreement with experimental determination.

		Error percentage of
		predictions Vs
	Model	experimental results
	Sachs (uniform stress)	14 %
	VPSC (mean field)	14 %
	EVPSC (mean field)	10 %
	CPFFT (full field)	12 %

Acknowledgments

This work was supported by the French State through the program "Investment in the future " operated by the National Research Agency (ANR) and referenced by ANR-11-LABEX-0008-01 (LabExDAMAS) and was also supported by the Région Grand Est. The authors would like to thank Dr. Nathalie GEY, Dr. Hana KRIAA and Dr. Jean-Jacques FUN-DENBERGER for fruitful discussions, and, Dr. Julien GUYON and Dr. Jean-Sébastien LECOMTE for technical support. Anthony D. Rollett is grateful for support from the LabEx DAMAS and the Department of Energy, Basic Energy Sciences, Mechanical Behavior and Radiation Effects Program led by Dr. John VETRANO, grant number DE-SC0019096.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Author Contributions

Meriem Ben Haj Slama has contributed to the conception, design of the work, acquisition, analysis, interpretation of data, and has drafted the work. Vincent Taupin, has contributed to the conception, interpretation of data, and has substantively revised the paper. Nabila Maloufi, has contributed to the conception, design of the work, and has substantively revised the paper. Kaustubh Venkatraman, has contributed to the analysis of data and the drafting of the work. Anthony D. Rollett, Ricardo A. Lebensohn, Stéphane Berbenni and Benoît Beausir have contributed to the analysis of data and the creation of new software used in the work. Antoine Guitton, has contributed to the conception, design of the work, interpretation of data, and has substantively revised the paper.

Supplementary materials

Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.mtla.2020.100996 .