domain: phys spi

text:
A Vattré

V Chiaruttini

Singularity-free theory and adaptive finite element computations of arbitrarily-shaped dislocation loop dynamics in 3D heterogeneous material structures

Keywords: Anisotropic elasticity theory, adaptive finite element method, heterogeneity, dislocation loops, dislocation shrinkage and propagation, short-range core reaction, dislocation-precipitate interaction, polycrystalline micropillar

The long-standing problem of arbitrarily-shaped discrete dislocation loops in three-dimensional heterogeneous material structures is addressed by introducing novel singularity-free elastic field solutions as well as developing adaptive finite element computations for dislocation dynamics simulations. The first framework uses the Stroh formalism in combination with the biperiodic Fourier-transform and dual variable and position techniques to determine the finite-valued Peach-Koehler force acting on curved dislocation loops. On the other hand, the second versatile mixed-element method proposes to capture the driving forces through dissipative energy considerations with domain integrals by means of the virtual extension principle of the surfacial discontinuities. Excellent agreement between theoretical and numerical analyses is illustrated from simple circular shear dislocation loops to prismatic dislocations with complicated simply-connected contours in linear homogeneous isotropic solids and anisotropic elastic multimaterials, which also serves as improved benchmarks for dealing with more realistic boundary-value problems with evolving dislocations. Examples of sophisticated dislocation applications include the short-range core reaction between intersecting dislocation loops in interaction with a spherical cavity, as well as the Orowan dislocation-precipitate bypass mechanism in a compressed micropillar of polycrystalline copper. The latter multiscale investigation spans three orders of magnitude in size scale, and is thus enabled by computationally efficient and robust adaptive mesh generation procedures for explicit dislocation propagation, interaction, and coalescence in three-dimensional materials and material structures.

Introduction

In crystalline materials, most mechanical properties, such as ductility, hardening and strengthening mechanisms, and fracture, are governed by the arrangement and the stress-driven motion of individual dislocations and their collective dynamics [START_REF] Kubin | Dislocations, mesoscale simulations and plastic flow[END_REF][START_REF] Cai | Imperfections in Crystalline Solids[END_REF], Bulatov et al., 2006). The fundamentals of dislocations are mainly studied top-down and bottom-up to describe plastic deformation processes as well as dislocation-mediated relaxation mechanisms in three-dimensional material structures. As topological defects that arise from translational lattice incompatibility in real solids, the closed dislocations are discrete Volterra dislocation loops [START_REF] Volterra | Sur l'équilibre des corps élastiques multiplement connexes[END_REF], for which the short-and long-range interactions with solute atoms and vacancies, free surfaces, grain boundaries and heterophase interfaces, cracks, inclusions and voids are characterized by continuous elastic distortion fields. The prediction of the internal dislocation-induced stresses within a complex full-scale environment is also ranged over a large portion of solid-state physics and continuum mechanics, which steadily require further in-depth treatment of various theoretical and numerical developments.

The present contribution is part of the long-standing problems, namely the continuum representation of discrete dislocation loops interacting with various strain-induced microstructural features in three-dimensional multimaterial and polycrystalline samples subjected to complex loading conditions. This offers i) novel singularity-free solutions for arbitrarily-shaped dislocation loops within the frame of the anisotropic elasticity theory in heterogeneous systems, and ii) a finite element methodology using adaptive remeshing techniques to treat more complicated configurations with large-scale computations whose the previous solutions are not yet determined, and certainly never will be for the most realistic industrial situations. The former therefore serves as a basis for comparison and validation for the numerical finite-element field solutions as well as the configurational forces for arbitrarily-shaped dislocation loops, while the latter proposes to handle computational dislocations dynamics in more advanced microstructures and sophisticated structures.

In recent decades, modern high-performance structures made of interface-dominated materials are used in a wide range of engineering applications from the aerospace and aviation industry with thermal barrier coatings on nickel-base superalloy substrates for turbine blades to micro-electrics with insulating films on semiconductor devices, through nuclear energy systems with radiation resistant metallic multilayers. At the microscopic scale, the presence of dislocation loops is also of paramount significance as elementary sources of internal stresses, while their motion and interaction are influenced by the elastic modulus mismatch and crystallographic misorientation across the internal interfaces. When the average radius of the dislocation loops is on the order of the representative length of the internal heterogeneities, the discrete distribution of dislocations cannot be described by constitutive relations. At this scale of observation, the complicated and important problem of dislocation interactions by the magnitude of the Burgers vectors.

Without the need to postulate the additive decomposition of the total solution fields into elastic and inelastic parts, the extended finite element method (Belytschko andBlack, 1999, Möes et al., 1999) has been developed for dislocation loops [START_REF] Belytschko | On XFEM applications to dislocations and interfaces[END_REF], Gracie et al., 2007, 2008[START_REF] Oswald | An extended finite element method for dislocations in complex geometries: thin films and nanotubes[END_REF][START_REF] Huang | An en efficient 2D discrete dislocation dynamics-XFEM coupling framework and its application to crystal plasticity[END_REF]. This computational method explicitly describes the dislocation loops as interior discontinuities with tangential enrichments [START_REF] Belytschko | Arbitrary discontinuities in finite elements[END_REF]. The corresponding slip is therefore prescribed on internal surfaces in terms of specific displacement boundary conditions, similar to the classical theory of dislocations. The non-singular elastic field solutions are also accurately captured with finite energy measures due to regularization functions at the dislocation cores, thus without assuming the discrete defects as platelet misfitting inclusions. The driving forces acting the dislocation loops are computed by a contour integral of the Eshelby energy-momentum tensor [START_REF] Eshelby | The elastic energy-momentum tensor[END_REF], which corresponds to the well-known J -integral [START_REF] Rice | Conserved integrals and energetic forces[END_REF][START_REF] Weertman | Dislocation Based Fracture Mechanics[END_REF][START_REF] Kirchner | The force on an elastic singularity in a non-homogeneous medium[END_REF][START_REF] Lubarda | Dislocation Burgers vector and the Peach-Koehler force: a review[END_REF]. The configurational force is the Peach-Koehler force in linear elasticity of dislocations [START_REF] Peach | The forces exerted on dislocations and the stress fields produced by them[END_REF]. Without remeshing techniques, the applicability and the capability of the extended finite element method with relatively simple topological configurations have been illustrated to interaction problems of dislocation loops with free surfaces as well as internal heterophase interfaces [START_REF] Belytschko | On XFEM applications to dislocations and interfaces[END_REF], Gracie et al., 2007, 2008[START_REF] Oswald | An extended finite element method for dislocations in complex geometries: thin films and nanotubes[END_REF]. However, the J -integral requires very fine meshes to resolve the configurational forces for dislocation propagation, which yields a notable limitation in numerical simulations using structured meshes of bilinear elements [START_REF] Gracie | On a new extended finite element method for dislocations: core enrichment and nonlinear formulation[END_REF], so that the application calculations are commonly limited to two-dimensional problems.

The remarkable versatility of the finite element method for heterogeneous multiscale and multiphysics problems motivates the present computationally efficient framework. In contrast to the previous extended finite element method, the present alternative approach is based on the finite element formulation using advanced and flexible adaptive remeshing strategies for fine mesh resolution close to arbitrarily-curved dislocation topologies and for dislocation loop dynamics in realistic materials and structures. The numerical framework is capable of computing the stress fields and the driving forces acting on dislocation loops with a high degree of accuracy and with significant precision, as compared to the novel theoretical elastic solutions for dislocation loops in the context of anisotropic hetero-elasticity theory in three dimensions. In what follows, a sign convention is adopted to define the positive and negative magnitude of the forces. The signed magnitude of the driving forces in the Peach-Koehler expression as well as of the configurational formulation is positive to describe the expansion of the dislocation loops, while the magnitude is reported as negative for dislocation shrinkage.

The content of the paper is organized as follows. Section 2 continues beyond previous contributions in the research area by treated dislocation loops in ideal hetereostructures formed by elastically anisotropic crystals using the Stroh formalism combining with biperiodic Fourier-transform and dual variable and position techniques. A single arbitrarily-shaped dislocation loop is embedded in multilayered systems made of dissimilar, linear and anisotropic plates with perfectly bounded internal interfaces. The standard stress singularity from the classical continuum theory is removed by spreading the Burgers vectors isotropically over a finite width on the glide planes. Extensions of this theoretical work by use of a monolithic finite element method with adaptive remeshing techniques are described in section 3. In this computational context, the fundamental Peach-Koehler force that governs the motion of the discrete dislocation loops is therefore captured by dissipative energy considerations. Section 4 illustrates an excellent agreement between the theoretical and numerical solutions in terms of non-singular stresses and driving forces. Thus, the finite element approach is used to explore more sophisticated microstructural situations, offering more realistic boundary-value description of dislocation dynamics in heterogenous materials and advanced multimaterial structures in presence of external and internal interfaces with arbitrary orientations.

2. Singularity-free theory of arbitrarily-shaped dislocation loops in multilayered structures Figure (1a) shows the ideal three-dimensional multilayered materials in Cartesian coordinates (x, y, z) = (x 1 , x 2 , x 3 ), for theoretical considerations. The complete system is composed of an arbitrary number of N -bonded orthotropic, dissimilar, rectangular, and linearly anisotropic elastic solids with perfectly connected internal interfaces. A singularity-free and arbitrarilyshaped dislocation loop is described by the Burgers vector over a horizontal planar cut surface S within a single layer, as displayed by the gray area in Fig. (1a). The individual and finite thickness is h λ = z λ -z λ-1 for the λ th layer, with λ = 1, . . . , N , while the origin O is located at one of the four corners on the bottom surface and all plates are defined in the positive z-region. The lower and upper interfaces of layer λ are also defined as z λ-1 and z λ , respectively, and the in-plane x x 1 -and y x 2 -directions are aligned along the horizontal edges of the plate boundaries. It follows that z 0 = 0 and z N = H at the bottom and top traction-free surfaces, respectively, with H the total thickness of the multilayers.

In the following, the general expressions for the displacements and tractions for each homogeneous plate are derived, and then propagated over the entire multilayered structures in presence of perfectly bonded interfaces. The prescribed displacement jump at the dislocation loops is consistently defined in the Fourier-transformed space, while the regularization treatment at the dislocation cores is able to determine the non-singular Peach-Koehler force for general curvilinear geometries, without introducing artificial cutoff distances to avoid the standard singularities.

Displacement and stress field solutions for each anisotropic elastic layer

In the absence of body forces and thermal effects, the classical stress equilibrium equation of anisotropic linear solids is fulfilled at any point x = [x 1 , x 2 , x 3 ] t of the multilayered structures, as follows where a comma stands for differentiation, with repeated indices denoting summation convention ranging from 1 to 3, unless stipulated otherwise. In eq. ( 1), the three-dimensional stress field σ kj (x) is related to the displacement field u k (x) by

σ λ kj,k (x) = 0 , (1)
h 1 h λ = z λ -z λ -1 h N z = z 0 = 0 z = z N = H L a y e
σ λ kj (x) = c λ kjml u λ m,l (x) , (2)
where the anisotropic elastic constants c λ kjml are independently defined with respect to the local crystallographic orientation of each layer λ, and are fully symmetric, i.e., c kjml = c jkml = c kjlm = c mlkj . Here and in the following, the superscript λ is omitted to avoid notational complexity, unless needed. As proposed in various dislocation-related problems for intrinsic and extrinsic defects [START_REF] Wu | Interface traction stress of 3D dislocation loop in anisotropic bimaterial[END_REF], Vattré, 2017a,b, Gao and Larson, 2015[START_REF] Vattré | Semicoherent heterophase interfaces with core-spreading dislocation structures in magneto-electro-elastic multilayers under external surface loads[END_REF], a biperiodic Fourier series formulation of the displacement vector is written as a linear superposition of harmonic plane-wave functions by

u k (x 1 , x 2 , x 3 = z) = Re ∑ η 1 ∑ η 2 e -i2π(η 1 x 1 +η 2 x 2 ) ũk (η 1 , η 2 , z) = Re ∑ η α e -i2πη α x α ũk (η α , z) , (3)
where i 2 = -1, and Re stands for the real part of a complex quantity, while α takes the value 1 or 2. Due to the mathematical completeness of Fourier series, the general field solutions are analyzed using plane waves e -i2πη α x α in eq. ( 3), where the unknown complex Fourier coefficients ũk (η α , z) are defined in the reciprocal space for each Fourier mode (η 1 , η 2 ) parallel to the perfectly bounded interface planes. Substituting eq. (3) to eq. ( 2) and then to eq. ( 1), three second-order differential equations for each layer are obtained in the Fourier-transformed domain, i.e.

4π 2 c jαkβ η α η β ũk (η α , z) + i2π (c jαk3 + c j3kα ) η α ũk,3 (η α , z) -c j3k3 ũk,33 (η α , z) = 0 , (4)
with β = {1, 2}. Following the algebraic manipulation from [START_REF] Pan | Static response of a transversely isotropic and layered half-space to general surface loads[END_REF], eq. ( 4) can be converted and cast into the following linear system of first-order differential equations (Vattré andPan, 2019, Liu et al., 2018), i.e.,

d dz ũ (η α , z) t (η α , z) = i2πη T -1 R t T -1 -4π 2 η 2 -Q + R T -1 R t i2πη RT -t ũ (η α , z) t (η α , z) , (5)
leading to six convenient relations between the expansion coefficients of the Fourier-transformed displacement and traction vectors, for which the latter is defined by tk (η α , z) = σkj (η α , z) n j = σk3 (η α , z), with unit normal n j x 3 . Equation ( 5) is met for each homogeneous layer, individually, within which the involved stiffness matrices Q, R, and T are given by

[Q] kj = c lkjs m l m s , [R] kj = c lkjs m l n s , and, [T] kj = c lkjs n l n s , (6)
while the components of the wave vector η are given by

η = [η 1 , η 2 , 0] t = η m = η 2 1 + η 2 2 m , (7)
where the corresponding magnitude η and direction m of the wavenumber vector are defined in the oblique plane basis spanned by (m(θ), n), i.e.,

m(θ) = m = [cos θ, sin θ, 0] t , and, n = [0, 0, 1] t , (8)
as commonly sought in the form of surface waves propagating in the direction m parallel to the internal and external surfaces, with m • n = 0 and |m| = |n| = 1. The general expression of the primary elastic solutions of eq. ( 5) is expressed as follows

ũk (η α , z) = e -i2πpηz a k tk (η α , z) = -i2πη e -i2πpηz b k , (9)
where p, a k = [a 1 , a 2 , a 3 ] t , and b k = [b 1 , b 2 , b 3 ] t are complex-valued unknowns. Substituting eq. ( 9) into eq. ( 5), these unknowns are determined by solving the eigenvalues and the corresponding complex eigenvectors of the Stroh eigenproblem [START_REF] Ting | Anisotropic elasticity: theory and applications[END_REF][START_REF] Stroh | Dislocations and Cracks in Anisotropic Elasticity[END_REF][START_REF] Stroh | Steady state problems in anisotropic elasticity[END_REF], i.e.,

-T -1 R t T -1 -Q + R T -1 R t -RT -t a b = p a b , (10)
recast into a classical linear eigensystem. The Stroh eigenvalues of eq. ( 10) as well as the eigenvectors are conveniently arranged such that Im p k > 0, and p k+3 = pk , since these complex quantities appear in complex conjugate pairs [START_REF] Eshelby | Anisotropic elasticity with applications to dislocation theory[END_REF].

Im stands for the imaginary part of a complex quantity, while the overbar denotes the complex conjugate. Superposing the six eigensolutions, the general expressions of the displacements and tractions in the Fourier-transformed domain are expressed in terms of the Stroh formalism, as follows

-i2πη ũ (η α , z) t (η α , z) = A Ā B B e -i2πp † η(z-z λ ) 0 3×3 0 3×3 e -i2π p † η(z-z λ-1 ) K 1 K 2 , ( 11
)
which are defined in any given layer λ bonded by the lower and upper interfaces at z = z λ-1 and z = z λ , respectively. In eq. ( 11), A and B are 3 × 3 eigenvector matrices defined by

A = a 1 , a 2 , a 3 B = b 1 , b 2 , b 3 = R t A + TA e i2πp † η(z-z λ ) , (12)
within which the z-dependent diagonal and exponential matrix is represented by

e i2πp † η(z-z λ ) = diag e i2πp 1 η(z-z λ ) , e i2πp 2 η(z-z λ ) , e i2πp 3 η(z-z λ ) , (13)
while K 1 and K 2 are two 3 × 1 complex and constant unknown matrices to be determined by the specific boundary-value conditions related to the presence of arbitrarily-shaped dislocation loops embedded in elastic multilayered structures.

Exact expression of the primary field solutions between two adjacent layers

The recent dual variable and position technique [START_REF] Pan | Green's functions for geophysics: a review[END_REF][START_REF] Vattré | Semicoherent heterophase interfaces with core-spreading dislocation structures in magneto-electro-elastic multilayers under external surface loads[END_REF][START_REF] Liu | General surface loading over layered transversely isotropic pavements with imperfect interfaces[END_REF]) is adopted in the present dislocation theory, in lieu of the conventional transfer matrix method for laminated systems [START_REF] Pan | Static response of a transversely isotropic and layered half-space to general surface loads[END_REF][START_REF] Thomson | Transmission of elastic waves through a stratified medium[END_REF][START_REF] Haskell | The dispersion of surface waves on multilayered media[END_REF][START_REF] Singh | Static deformation of a multilayered half-space by internal sources[END_REF] that might arise numerical instability issues, such as spurious oscillations for large amplitudes of frequencies or half-wave numbers [START_REF] Pan | Green's functions for geophysics: a review[END_REF], Vattré and Pan, 2021, Vattré et al., 2021). Substituting z by z λ-1 into the linear system in eq. ( 11), the general solutions in terms of displacements and tractions read

-i2πη ũ (η α , z λ-1 ) t (η α , z λ-1 ) = A Ā B B e i2πp † ηh λ 0 3×3 0 3×3 I 3×3 K 1 K 2 , (14)
while similarly substituting z = z λ into eq. ( 11), analogous relations are obtained, i.e.,

-i2πη ũ (η α , z λ ) t (η α , z λ ) = A Ā B B I 3×3 0 3×3 0 3×3 e -i2π p † ηh λ K 1 K 2 , ( 15
)
such that both unknown complex vectors K 1 and K 2 in eqs. ( 14) and ( 15) can be removed to meet cross field relations (Vattré andPan, 2019, Liu et al., 2018) at z = z λ-1 and z = z λ . Thus, it follows

-i2πη ũ (η α , z λ-1 ) t (η α , z λ ) = S λ 6×6 -i2πη ũ (η α , z λ ) t (η α , z λ-1 ) = S λ 11 S λ 12 S λ 21 S λ 22 -i2πη ũ (η α , z λ ) t (η α , z λ-1 ) , ( 16
)
where the relevant submatrices of the six-dimensional matrix S λ 6×6 are defined by

S λ 6×6 = S λ 11 S λ 12 S λ 21 S λ 22 = A e i2π p † η h λ Ā B B e -i2π p † η h λ A Ā e -i2π p † η h λ B e i2π p † η h λ B -1 , (17)
for which the inverse of the 2 × 2 block matrices over the field of complex quantities in the right-hand side can be explicitly written using the Banachiewicz-Schur form, since the matrix S λ 6×6 contains real-valued elements [START_REF] Baksalary | Generalized inverses of partitioned matrices in Banachiewicz-Schur form[END_REF]. On the other hand, the propagation relations of the expansion coefficient solutions at both interfaces z λ and z λ+1 for the adjacent layer λ + 1 are specified by

-i2πη ũ (η α , z λ ) t (η α , z λ+1 ) = S λ+1 6×6 -i2πη ũ (η α , z λ+1 ) t (η α , z λ ) = S λ+1 11 S λ+1 12 S λ+1 21 S λ+1 22 -i2πη ũ (η α , z λ+1 ) t (η α , z λ ) , (18)
similarly constructed as eq. ( 16), including the associated built-in elements in terms of the Stroh eigensolutions and thicknesses. Furthermore, assuming that the interface at z λ between the two adjacent layers is perfectly bonded, the continuous Fouriertransformed displacement and traction vectors at

z = z λ read ũ (η α , z = z λ ) + -= ũ (η α , z λ+ ) -ũ (η α , z λ-) = 0 3×1 t (η α , z = z λ ) + -= t (η α , z λ+ ) -t (η α , z λ-) = 0 3×1 , (19)
which lead to important recursive relations between interfaces z λ-1 and z λ+1 by combining eqs. ( 16) and (18) together, i.e.,

-i2πη ũ (η α , z λ-1 ) t (η α , z λ+1 ) = S λ:λ+1 6×6 -i2πη ũ (η α , z λ+1 ) t (η α , z λ-1 ) = S λ:λ+1 11 S λ:λ+1 12 S λ:λ+1 21 S λ:λ+1 22 -i2πη ũ (η α , z λ+1 ) t (η α , z λ-1 ) , ( 20
)
where the superscript λ:λ+1 means the resulting propagation matrix from layer λ to layer λ+1, and the layer-to-layer submatrices S λ:λ+1 γδ of the six-dimensional block matrix S λ:λ+1 6×6 are expressed as

              
I 3×3 -S λ 21 S λ+1 12 -1 S λ 22 , (21)
as analogously proposed in recent multilayered problems under multiple surface loadings [START_REF] Liu | General surface loading over layered transversely isotropic pavements with imperfect interfaces[END_REF] and specific internal conditions for semicoherent interfaces [START_REF] Vattré | Semicoherent heterophase interfaces with core-spreading dislocation structures in magneto-electro-elastic multilayers under external surface loads[END_REF].

Recursive relations in presence of a dislocation loop embedded in multilayered structures

In contrast with the perfectly bonded interfacial conditions expressed in eq. ( 19), the horizontal and finite dislocation loop is characterized by considering specific boundary conditions in terms of continuous tractions and displacement discontinuity with respect to the given Burgers vector b

k = [b 1 , b 2 , b 3 ] t .
Assuming that the layer λ is the source layer that encompasses the closed dislocation loop at z = z s , the standard condition in the three-dimensional physical domain reads

u k (x 1 , x 2 , x 3 = z s+ ) -u k (x 1 , x 2 , x 3 = z s-) = b k χ S (x 1 , x 2 ) t k (x 1 , x 2 , x 3 = z s+ ) -t k (x 1 , x 2 , x 3 = z s-) = 0 , (22)
where the indicator function χ S (x 1 , x 2 ) of the set S yields χ S (x 1 , x 2 ) = 1 if (x 1 , x 2 ) ∈ S, and χ S (x 1 , x 2 ) = 0, otherwise. The corresponding boundary conditions in the Fourier-transformed domain are written as follows

ũk (η α , z = z s+ ) -ũk (η α , z = z s-) = ũp k (η α , b k ) tk (η α , z = z s+ ) -tk (η α , z = z s-) = 0 , (23)
where ũp k (η α , b k ) is the complex-valued function related to the prescribed displacement discontinuity of arbitrary Burger vector in the reciprocal space.

In the context of dislocation loops embedded in multilayered structures, the global propagation submatrices in eq. ( 21) from the bottom surface at z = z 0 = 0 to the top surface to z = z N = H can be partitioned to explicitly determine the displacement and traction solutions at any z-level including the conditions for the dislocation loop and the perfectly bonded interfaces. Thus, the solutions in the Fourier-transformed domain at z f in layer λ can be obtained from eq. ( 20) by propagating first the transformed displacement and traction vectors from the bottom surface z = 0 to the lower side of the dislocation loop at z = z s-, as follows

-i2πη ũ (η α , z = 0) t (η α , z = z s-) = S 1:λ 11 S 1:λ 12 S 1:λ 21 S 1:λ 22 -i2πη ũ (η α , z = z s-) t (η α , z = 0) , (24)
and then, propagating the solution from the upper side at z s+ of the dislocation loop to the top surface z = H to similarly obtain the subsequent relations

-i2πη ũ (η α , z = z s+ ) t (η α , z = H) = S λ:N 11 S λ:N 12 S λ:N 21 S λ:N 22 -i2πη ũ (η α , z = H) t (η α , z = z s+ ) = -i2πη ũ (η α , z = z s-) -i2πη ũp (η α , b k ) t (η α , z = H) , (25)
for which the dislocation boundary condition in terms of displacement in eq. ( 23) is incorporated in the right-hand side. Both eqs. ( 24) and ( 25) can be conveniently recast into the following linear system, i.e.,

    0 3×3 -I 3×3 S λ:N 11 S λ:N 12 0 3×3 0 3×3 S λ:N 21 S λ:N 22 -I 3×3 S 1:λ 11 0 3×3 0 3×3 0 3×3 S 1:λ 21 0 3×3 -I 3×3         -i2πη ũ (η α , z = 0) -i2πη ũ (η α , z = z s-) -i2πη ũ (η α , z = H) t (η α , z = z s-)     =     -i2πη ũp (η α , b k ) t (η α , z = H) -S 1:λ 12 t (η α , z = 0) -S 1:λ 22 t (η α , z = 0)     , (26)
where the submatrices S 1:λ γδ and S λ+1:N γδ are defined in eq. ( 21). The system in eq. ( 26) is therefore solved for each Fourier mode (η 1 , η 2 ) to determine the displacement fields at the three involved locations, i.e., ũ (η α , z = 0), ũ (η α , z s-) and ũ (η α , z = H), as well as the internal traction, i.e., t (η α , z = z s ) at z = z s-= z s , with respect to the given boundary conditions at the external surfaces, i.e., t (η α , z = 0), t (η α , z = H) as well as the planar dislocated interface, i.e., ũp (η α , b k ) at z = z s .

Importantly, for any field point z f below the location of the dislocation loop embedded in the layer λ f , i.e., 0 < z f < z s-, the recursive relations in eq. ( 20) are conveniently split as follows

           -i2πη ũ (η α , z = 0) t (η α , z = z f ) = S 1:λ f 11 S 1:λ f 12 S 1:λ f 21 S 1:λ f 22 -i2πη ũ (η α , z = z f ) t (η α , z = 0) -i2πη ũ (η α , z = z f ) t (η α , z = z s ) = S λ f :s 11 S λ f :s 12 S λ f :s 21 S λ f :s 22 -i2πη ũ (η α , z = z s-) t (η α , z = z f ) , (27)
which can also be recast into the following linear system of equations, i.e.,

     S 1:λ f 11 0 3×3 -I 3×3 0 3×3 S 1:λ f 21 -I 3×3 0 3×3 0 3×3 I 3×3 S λ f :s 12 0 3×3 0 3×3 0 3×3 S λ f :s 22 0 3×3 -I 3×3          -i2πη ũ (z = z f ) t (z = z f ) -i2πη ũ (z = 0) t (z = z s )     =      -S 1:λ f 12 t (z = 0) -S 1:λ f 22 t (z = 0) -2πη S λ f :s 11 ũ (η α , z = z s-) -2πη S λ f :s 21 ũ (η α , z = z s-)      , ( 28
)
and be solved to determine the displacement and traction fields at z = z f , using the known displacement field solutions ũ (η α , z s-) at z = z s-= z s in the right-hand side, previously obtained from the resolution of eq. ( 26). On the other hand, for any field solution at point z f such that z s+ < z f < H in the layer λ f , a similar system as eq. ( 28) leads to

     S s:λ f 11 0 3×3 0 3×3 S s:λ f 12 S s:λ f 21 -I 3×3 0 3×3 S s:λ f 22 -I 3×3 S λ f :N 12 S λ f :N 11 0 3×3 0 3×3 S λ f :N 22 S λ f :N 21 0 3×3          -i2πη ũ (η α , z = z f ) t (η α , z = z f ) -i2πη ũ (η α , z = H) t (η α , z = z s )     =     -i2πη ũ (η α , z = z s-) -i2πη ũp (η α , b k ) 0 3×1 0 3×1 t (η α , z = H)     , (29)
within which the known displacement field solutions ũ (η α , z s-) is analogously introduced. For the present particular case of free-standing multilayered structures, the traction-free boundary conditions at the bottom and top surfaces, i.e., t (η α , z = 0) = t (η α , z = H) = 0, are imposed in eqs. ( 26), (28), and (29), without loss of generality, while the singularity-free displacement discontinuity ũp (η α , b k ) is formulated in the next section.

Prescribed displacement discontinuity for core-spreading dislocation loops

The displacement discontinuity represented by the Burgers vectors of the closed dislocation loops in the Fourier reciprocal space is analytically derived for two shapes, namely the rectangular and the elliptical Volterra dislocation loops in sections 2.4.1 and 2.4.2, respectively. The closed-form solutions are obtained using an isotropic core-spreading procedure to expand the compact planar dislocation cores over a finite area, characterizing the Burgers vector distribution within the physical dislocation cores. The general arbitrarily-shaped dislocation loops with simply-connected contours, which can be characterized by parametric equations, are treated using a univariate Gauss-Legendre cubature based on the Green integration formula, as detailed in section 2.4.3.

Rectangular dislocation loops

For a given dislocation loop with rectangular shapes R that is composed of four straight segments located at x 3 = z = z s , x 1 = ±R 1 , and x 2 = ±R 2 , with the Burgers vector b k and the center at [x 0 , y 0 , z s ] t , the displacement discontinuity in the Fouriertransformed domain is given by the following surface integral, i.e.

ũp k (η α , b k ) = ˆR b * k (x 1 , x 2 ) e i2πη 1 x 1 e i2πη 2 x 2 dx 1 dx 2 = ˆR1 +x 0 -R 1 +x 0 ˆR2 +y 0 -R 2 +y 0 b * k (x 1 , x 2 ) e i2πη 1 x 1 e i2πη 2 x 2 dx 1 dx 2 , ( 30
)
where b * k is the core-spreading Burgers vector associated with the compact vector b k . Such compact dislocation cores can therefore be spread out by convoluting the discontinuity displacement conditions with specific two-dimensional weighted function on the dislocation plane [START_REF] Vattré | Semicoherent heterophase interfaces with core-spreading dislocation structures in magneto-electro-elastic multilayers under external surface loads[END_REF]. The isotropic Gaussian function is used to introduce a continuous distribution of the Burgers vectors as defined by

(x 1 , x 2 ) = 1 πr 2 0 e -(x 2 1 +x 2 2 )/r 2 0 , with : ˆ∞ -∞ ˆ∞ -∞ (x 1 , x 2 ) dx 1 dx 2 = 1 , (31)
with r 0 ≥ 0 the core-spreading dislocation parameter, so that convolution term in eq. ( 30) with respect to (x 1 , x 2 ) leads to

ũp k (η α , b k ) = b k ˆR1 +x 0 -R 1 +x 0 ˆR2 +y 0 -R 2 +y 0 e i2πη 1 x 1 e i2πη 2 x 2 * (x 1 , x 2 ) dx 1 dx 2 , ( 32
)
where * is the convolution product. By use of standard tables of integrals, the integral over the rectangle can be explicitly determined by the following analytical closed form, i.e.

ũp k (η α , b k ) = 4R 1 R 2 b k sinc (2πη 1 R 1 ) sinc (2πη 2 R 2 ) e i2πη 1 x 0 -π 2 η 2 1 r 2 0 e i2πη 2 y 0 -π 2 η 2 2 r 2 0 , (33)
which is also substituted into eqs. ( 26) and ( 29), as stated therefore as internal boundary conditions. For physical reasons, a multivariate anisotropic Gaussian distribution, consequently including more core-spreading parameters, can be convoluted in the same straightforward manner, without subsequent mathematical difficulties, as recently proposed by Lazar, Po, and co-authors [START_REF] Po | A non-singular theory of dislocations in anisotropic crystals[END_REF][START_REF] Lazar | Three-dimensional nonlocal anisotropic elasticity: a generalized continuum theory of Ångström-mechanics[END_REF][START_REF] Lazar | The non-singular Green tensor of Mindlin's anisotropic gradient elasticity with separable weak non-locality[END_REF]. The calibration of these additional length-dependent parameters, however, requires ad-hoc simulations from the first principles or atomistics, which is beyond the scope of the present work.

Elliptical dislocation loops

Following the line integration derivation of dislocation fields for parametric dislocation curves using polar coordinates by [START_REF] Han | Fields produced by three-dimensional dislocation loops in anisotropic magneto-electro-elastic materials[END_REF], the two-dimensional Fourier-transformed displacement discontinuity is similarly formulated for a horizontallyoriented elliptical loop E with major and minor semi-axes a 1 and a 2 centered about [x 0 , y 0 , z s ] t , as follows

ũp k (η α , b k ) = b k ˆE e i2πη 1 x 1 -π 2 η 2 1 r 2 0 e i2πη 2 x 2 -π 2 η 2 2 r 2 0 dx 1 dx 2 = 2πR E η E J 1 (η E R E ) b k e i2πη 1 x 0 -π 2 η 2 1 r 2 0 e i2πη 2 y 0 -π 2 η 2 2 r 2 0 , (34)
where J 1 is the Bessel function of the first kind and order 1, while

η E = 2π R E (η 1 a 1 ) 2 + (η 2 a 2 ) 2 , and R E = √ a 1 a 2 .

Arbitrarily-shaped dislocation loops with simply-connected contours

For the general integration domain D over arbitrarily-shaped dislocation loops, closed simply-connected contours with convex and nonconvex geometrical regions are characterized by simply-connected polygons and described by parametric curves in closed-form or represented using B-splines and non-uniform rational B-splines. In practice, the surface integral over D can be achieved using efficient and suitable cubature formulas for planar regions, as commonly derived for polynomials over regular polygons [START_REF] Nooijen | Symmetric numerical integration formulas for regular polygons[END_REF] and complex polygons [START_REF] Sommariva | Product Gauss cubature over polygons based on Green's integration formula[END_REF]. Based on the Green integration theorem, i.e., Stokes's theorem in two dimensions, Gauss-related cubatures over polygons are exact for bivariate polynomials of degree at most 2M D using ∼ M 2 D nodes, reducing the integral over the polygons to a sum of edge integrations. Thus, the boundary ∂D is considered as a Jordan curve and is also described counterclockwise by the sequence of vertices

V i = {x i 1 , x i 2 }, with i = 1, . . . , N D , such that ∂D = N D i=1 [V i , V i+1 ] = N D i=1 ν i , with : V 1 = V N D , (35)
so that N D is the total number of edges on the boundary of a simply-connected polygonal domain. By virtue of the Green formula in the context of cubature rules, the surface integral of the prescribed displacement discontinuity ũp (η α , b k ) in the reciprocal space is transformed into a one-dimensional problem, as follows

ũp k (η α , b k ) = b k ˆD g (x 1 , x 2 ) dx 1 dx 2 = b k L ∑ ν i ∈∂D ˆ[V l ,V l+1 ] G (x 1 , x 2 ) dx 2 ≈ b k ∑ ν i ∈∂D M D +1 ∑ j=1 M D ∑ m=1 ijm g x i 1 (τ j , τ m ) , x i 2 (τ j ) , (36)
where

G(x 1 , x 2 ) is any fixed x 1 -primitive of g(x 1 , x 2 ), i.e., G (x 1 , x 2 ) = ˆx1 g (w 1 , x 2 ) dw 1 , (37)
where

x i 1 is a fixed convex combination of the x i 1 coordinates. The second discrete sum over j in the right-hand side of eq. ( 36) is performed for M D + 1 points along a base-edge, whose intersection is the connected with the domain D, while the sum over m is associated with M D -oriented orthogonal lines to the base-edge [START_REF] Sommariva | Product Gauss cubature over polygons based on Green's integration formula[END_REF]. A minimum decomposition of a given simple polygon into convex subpolygons is usually carried out for nonconvex and multiply-connected polygons [START_REF] Gerdjikov | Decomposing a simple polygon into pseudo-triangles and convex polygons[END_REF].

When the integrand g(x 1 , x 2 ) is a bivariate polynomial of degree at most 2M D -1, the quadrature in the right-hand side of eq. ( 36) is exact with M D nodes and gives a polynomial of degree at most 2M D . According to the present dislocation-based problem, the integrand is a non-polynomial and complex-valued function, defined by

g (x 1 , x 2 ) = e i2πη 1 x 1 -π 2 η 2 1 r 2 0 e i2πη 2 x 2 -π 2 η 2 2 r 2 0 = g 1 (η α ) e ig 2 (η α ,x α ) with : g 1 (η α ) = e -π 2 η 2 α r 2 0 g 2 (η α , x α ) = 2πη α x α , (38)
where g(x 1 , x 2 ) is absolutely continuous on a domain D with smooth boundary described counterclockwise due to the corespreading procedures from eq. ( 31) at the dislocation contours, for each Fourier mode η α > 0 and r 0 0. In eq. ( 38), the pre-exponential function g 1 (η α ) decays rapidly, while g 2 (η α , x α ) is a bivariate function with respect to x α . For the case of compact dislocation cores, i.e., with r 0 = 0 in eq. ( 38), the corresponding integral is ordinarily highly oscillatory because of the presence of the exponential terms of imaginary arguments. In such problems where the multivariate integrand is oscillatory close to the dislocation cores, the numerical evaluation can be difficult for large amplitudes of frequencies or half-wave numbers, where a prohibitive number of points is needed in standard quadratures based on high order Gauss-Legendre integration rules (Hamming, 1973, Laghrouche and[START_REF] Laghrouche | Short wave modelling using special finite elements[END_REF]. Nevertheless, for spreading dislocation cores, the prescribed displacement field becomes infinitely differentiable over D, without oscillatory regimes, as illustrated for interfacial dislocations [START_REF] Vattré | Semicoherent heterophase interfaces with core-spreading dislocation structures in magneto-electro-elastic multilayers under external surface loads[END_REF] and heterogeneously distributed loads at the external surfaces (Vattré and Pan, 2021). With an excellent cubature convergence, the singularity-free property of the present displacement boundary conditions lead tractable the Gauss-Legendre quadrature rules enriched by the Green integral formula for the present theory of arbitrarily-shaped and curved dislocation loops.

The nodes {τ j } M D j=1 and the corresponding weights {µ j } M D j=1 are also associated with the Gauss-Legendre quadrature rules, for which the elements in eq. ( 38) are defined by

x i 1 (τ j , τ m ) = 1 2 x i 1 (τ j ) -τ m + 1 2 x i 1 (τ j ) + with : x i 1 (τ j ) = 1 2 x i+1 1 -x i 1 τ j + 1 2 x i+1 1 + x i 1 x i 2 (τ j ) = 1 2 x i+1 2 -x i 2 τ j + 1 2 x i+1 2 + x i 2 , ijm = 1 4 x i+1 2 -x i 2 x i 1 (τ j ) -µ j µ m , (39)
such that the overall number of quadrature nodes is bounded by (M D + 1) 2 N D . The Gauss-Legendre quadrature rules enriched by the Green formula avoid explicit partition of the polygonal domain D into triangles or quadrangles, as used in the finite element approaches in section 3, reducing intensively the computational time for the numerical integrations in ideal multilayered structures.

General expression for the non-singular Peach-Koehler force

The driving Peach-Koehler force per unit length F k (x ∂D ) at any point x ∂D = x 1 ∂D , x 2 ∂D , z s t on the boundary ∂D of the a dislocation loop embedded in the layer λ is related to the three-dimensional and non-singular stress field for r 0 0, as follows

F k (x ∂D ) = ε kil b j ξ l (x ∂D ) σ ij (x ∂D ) , (40)
where ε kil is the alternating tensor, and ξ l (x ∂D ) is the unit tangent vector of the dislocation at x ∂D . Solving eq. ( 29) with respect to the prescribed displacement discontinuity given by eq. ( 36) for the general case of simply-connected dislocation loops, the primary traction field tk (η α , z f ) at any z f in the reciprocal space is completed by the in-plane stress components σs ij (η α , z f ) to compute by Fourier inversion the Peach-Koehler force in eq. ( 40) at z f = z s . In the Fourier-transformed domain, the in-plane stress solutions are deduced from eq. ( 1) written in the reciprocal space, i.e.,

σs ij (η α , z f ) = -i2πη c ijkl m l ũk (η α , z f ) + c ijkl n l ũk,3 (η α , z f ) = [ σs 11 (η α , z f ) , σs 12 (η α , z f ) , σs 22 (η α , z f )] t , (41)
with i = {1, 2}, j = {1, 2}, and i ≤ j. Thus, eq. ( 41) depends on the primary displacement field ũk (η α , z f ) that is determined by eq. ( 29), only, for which the derivative with respect to x 3 is given by ũk,3

(η α , z f ) = [c ijkl n l n i ] -1 ( tj (η α , z f ) + i2πη c ijkl m l n i ũk (η α , z f )) , (42)
so that eq. ( 41) reads in vector-tensor form as

σs (η α , z f ) = -i2πη M 1 ũ (η α , z f ) + M 2 T -1 t (η α , z f ) + i2πη R t ũ (η α , z f ) , (43)
where both matrices M 1 and M 2 are explicitly defined by

M 1 =   c 11 m 1 + c 16 m 2 c 16 m 1 + c 12 m 2 c 15 m 1 + c 14 m 2 c 61 m 1 + c 66 m 2 c 66 m 1 + c 62 m 2 c 65 m 1 + c 64 m 2 c 21 m 1 + c 26 m 2 c 26 m 1 + c 22 m 2 c 25 m 1 + c 24 m 2   M 2 =   c 15 c 14 c 13 c 65 c 64 c 63 c 25 c 24 c 23   , (44)
which completes the total stress field σij (η α , z f ) generated by a dislocation loop in the Fourier-transformed space. The latter non-singular stress components are therefore determined into the physical domain by inverse Fourier transforms, similarly to eq. ( 3), as follows

σ ij (x 1 , x 2 , x 3 = z f ) = Re ∑ η α e -i2πη α x α σij (η α , z f ) = Re ∑ η α e -i2πη α x α   σs 11 (η α , z f ) σs 12 (η α , z f ) t1 (η α , z f ) σs 12 (η α , z f ) σs 22 (η α , z f ) t2 (η α , z f ) t1 (η α , z f ) t2 (η α , z f ) t3 (η α , z f )   , (45)
so that the non-singular Peach-Koehler force can be evaluated at x ∂D = x 1 ∂D , x 2 ∂D , z s t , by substituting the solutions from eqs. ( 29) and ( 43) as well as eq. ( 45) in eq. ( 40).

2.6. Comment on the limitations encountered with the preliminary theoretical analysis

Although the novel singularity-free field solutions for arbitrarily-shaped planar dislocation loops in multilayered anisotropic solids are relatively general in contrast with the elastic elements in standard dislocation dynamics simulations, the present work suffers from some basic limitations, namely the absence of motion and interaction of these dislocation loops with other dislocations [START_REF] Arsenlis | Enabling strain hardening simulations with dislocation dynamics[END_REF], Bulatov et al., 2006), precipitates [START_REF] Queyreau | Orowan strengthening and forest hardening superposition examined by dislocation dynamics simulations[END_REF][START_REF] Vattré | Orientation dependence of plastic deformation in nickel-based single crystal superalloys: discrete-continuous model simulations[END_REF], Gao et al., 2015b) and solutes (Aubry et al., 2015b), free surfaces [START_REF] Khraishi | Free-surface effects in 3D dislocation dynamics: formulation and modeling[END_REF]Zbib, 2002, Weinberger andCai, 2008) and grain boundaries [START_REF] Fan | The role of twimming deformation on the hardening response of polycrystalline magnesium from discrete dislocation dynamics simulations[END_REF], Cho et al., 2020). These features are important ingredients in gaining insights concerning the mechanical properties of materials, including work hardening in bulk solids as well as strain dislocation-mediated relaxation in polycrystalline metals under complex loading conditions.

In the next section, an advanced finite element framework using adaptive remeshing strategies is therefore developed to focus on more realistic boundary-value problems of dislocation loops in bounded crystals with complex internal microstructures, e.g., with randomly-oriented grains in presence of voids and heterophase inclusions. The numerical approach is implemented in the Z-set finite element platform [START_REF] Besson | Large scale object-oriented finite element code design[END_REF], and is strongly inspired from the fracture mechanics, which offers accurate solutions for crack propagation with explicit topological changes embedded in three-dimensional materials. Thus, the development for crack propagation, adaptive remeshing, as well as domain integral post-processing have been adapted for the present boundary-value problems with evolving dislocation loops. The excellent accuracy is evaluated by comparing the numerical stresses and driving forces with the previous theoretical solutions, reflecting i) the importance to formulate sophisticated field solutions to strengthen the numerical finite element techniques with various complex dislocation-related situations, as well as ii) the relevance to discontinuously track the discrete dislocation movements by efficient remeshing algorithms with a sequence of triangular meshes. Without resorting to the eigenstrain formalism or to the extended finite element method for dislocation loops, the numerical part is considered as the most comprehensive extension with respect to the preliminary and companion theoretical part as well as to the robust workhorse tool for future computational studies, thus providing new perspectives and enriching incomplete roadmaps recently addressed (van der Giessen et al., 2020).

Finite element analysis of discrete dislocation dynamics using adaptive remeshing

As illustrated in Fig. (1b), the linear elasticity problem consists of a bounded polycrystalline microsample with randomlyoriented grains in presence of a dislocation loop interacting with a heterophase inclusion in three dimensions. The fundamental Orowan mechanism leaving a dislocation loop around the bypassed particle in such a complex microstructural situation can also be modeled by the following finite element methodology using adaptive remeshing procedures. In contrast with Fig. ( 1a), the entire domain Ω is occupied by N grains Ω λ , with λ = 1, . . . , N , including internal grain boundaries with different shapes and characters. The grain boundaries as well as the material particles embedded in Ω λ are treated as impenetrable to the dislocation loops, so that the dislocation motion is confined inside the local grains. The solid Ω is also surrounded by a closed external boundary ∂Ω, where surface tractions and displacements are applied on disjoint portions, while the interior of the Ω λ th grain incorporates K λ discrete dislocation loops S λ κ , with κ = 1, . . . , K λ . The versatility of the conforming finite element analysis leads the first problem of Fig. (1a) to a particular case for the numerical methodology pictured in Fig. (1b). In the absence of body forces and thermal effects, the numerical approach is also capable of handling the same equilibrium equations of heterogeneous materials as in the previous section, also limited to quasi-static loads in the present work.

In the following, the variational formulation of the present boundary-value problem in presence of planar dislocation loops is presented in the context of the conventional finite element method, within which the meshes conform to the individual topological discontinuities. The calculation procedure of the driving Peach-Koehler forces exerted on the dislocation loops is then formulated by energetic arguments, while the discrete solution scheme is introduced using advanced and fully automatic remeshing techniques, specifically derived to investigate the three-dimensional interaction, propagation and coalescence of arbitrarily-shaped dislocation loops.

Variational formulation for the dislocation boundary-value problem

In three dimensions, the elastic boundary-value problem for a dislocated solid Ω = N λ=1 Ω λ , with closure ∂Ω, is formulated by the governing equations and boundary conditions for each subdomain Ω λ , as follows

               σ λ kj,k (x) = c λ kjml u λ m,lk (x) = 0 , ∀ x ∈ Ω λ σ λ kj (x) n ext j (x) = t N k (x) , ∀ x ∈ ∂Ω λ ∩ ∂Ω t u λ k (x) = u D k (x) , ∀ x ∈ ∂Ω λ ∩ ∂Ω u σ λ kj (x) n λ κ j + -= 0 , ∀ x ∈ S λ κ u λ k (x) + -= b λ κ k χ * (x) , ∀ x ∈ S λ κ , (46)
for which the first equation is defined by eq. ( 1), while the second and third equations are related to the prescribed Neumann and Dirichlet boundary conditions on the subsets ∂Ω λ ∩∂Ω t and ∂Ω λ ∩∂Ω u , respectively, with n ext (x) the outer normal vector at any point x that belongs to ∂Ω. In eq. ( 46), the two last boundary conditions are associated with the local traction and displacement requirements on each discrete dislocation loop S λ κ , which are described by the corresponding unit vector normal n λ κ to the local crystallographic planes and the Burgers vector b λ κ of S λ κ , as specified in a more general manner than in eq. ( 22). The variety of Burgers vectors in crystal types of body-centered cubic, face-centered cubic and hexagonal close-packed grains can also be incorporated adequately into the numerical calculations if required. In the finite element framework, the regularization linear function for dislocation cores in eq. ( 46) is parametrized by

χ * (x) = min (D(x, x ∂D )/r 0 , 1) , (47)
where D(x, x ∂D ) is the Euclidean distance between x and the position of the dislocation fronts ∂S λ κ at x ∂D , using the dislocation core parameter r 0 ≥ 0. Compared to the previous non-singular field theory, where the core-spreading procedure is obtained by convoluting the discontinuity displacement conditions with the Gaussian function in eq. ( 31) for straightforward mathematical operations into the reciprocal space, a much simpler linear function is conveniently used by eq. ( 47) in the present finite element formulation. Apart from the dislocation loops S λ κ , the traction stress and displacement fields are continuous across all internal surfaces, although physically motivated interface conditions [START_REF] Vattré | Determining the Burgers vectors and elastic strain energies of interface dislocation arrays using anisotropic elasticity theory[END_REF][START_REF] Vattré | Mechanical interactions between semicoherent heterophase interfaces and free surfaces in crystalline bilayers[END_REF][START_REF] Vattré | Elastic interactions between interface dislocations and internal stresses in finite-thickness nanolayered materials[END_REF]) could be included in the boundary-value problem.

The numerical computation of the strong form of the linear system in eq. ( 46) needs to agglomerate the equations upon a weak-form displacement integral equation for three-dimensional dislocation field analysis. Using the boundary conditions on dislocation loops as well as the Dirichlet boundary conditions, the trial function space of the kinematically admissible displacement fields U, with sufficient regularity, is introduced, i.e.

U = u λ k (x) u λ k (x) = u D k (x) , ∀ x ∈ ∂Ω λ ∩ ∂Ω u ∧ u λ k (x) + -= b λ κ k χ * (x) , ∀ x ∈ S λ κ , (48)
such that eq. ( 46) is formally equivalent to a single integral equation of the boundary-value problem of discrete dislocation loops, suitable for domain discretization techniques, as follows

N ∑ λ=1 ˆΩλ c λ kjml u λ m,l (x) v λ k,j (x) dΩ = N ∑ λ=1 ˆ∂Ω λ ∩∂Ω t t N k (x) v λ k (x) dΣ , ∀ v λ k (x) ∈ U 0 , (49)
where U 0 is the associated space for the test functions.

Domain integral evaluation of the configurational Peach-Koehler force

In the context of finite element analysis based on the previous variational formulation, the evaluation of the driving forces on the dislocation loop fronts is intended to be performed by path-independent integrals, as originally derived by [START_REF] Rice | Conserved integrals and energetic forces[END_REF] for crack problems in fracture mechanics, and then applied to dislocation growth [START_REF] Belytschko | On XFEM applications to dislocations and interfaces[END_REF], Gracie et al., 2007, 2008[START_REF] Oswald | An extended finite element method for dislocations in complex geometries: thin films and nanotubes[END_REF][START_REF] Huang | An en efficient 2D discrete dislocation dynamics-XFEM coupling framework and its application to crystal plasticity[END_REF]. The configurational Peach-Koehler force for dislocations in eq. ( 40) is also related to the dislocation energy release rate with respect to the newly formed loops due to a local and virtual perturbation of the dislocation fronts. According to the configurational mechanics from [START_REF] Eshelby | The elastic energy-momentum tensor[END_REF][START_REF] Eshelby | The force on an elastic singularity[END_REF][START_REF] Eshelby | The continuum theory of lattice defects[END_REF], the total potential energy Π of the externally loaded domain Ω containing a collection of dislocation loops S λ κ is defined as the difference of the strain energy W stored in the solid and the external load potential W ext , i.e.

Π = W -W ext = W bulk + W dis -W ext = 1 2 N ∑ λ=1 ˆΩλ c λ kjml u λ m,l (x) u λ k,j (x) dΩ + N ∑ λ=1 K λ ∑ κ=1 ˆSλκ σ λ kj (x) n λ κ j b λ κ k χ * (x) dΣ - N ∑ λ=1 ˆ∂Ω λ ∩∂Ω t t N k (x) u λ k (x) dΣ , (50)
for which W is decomposed into the bulk energy W bulk and the interfacial dislocation-related energy W dis , by virtue of the superposition principle. The total potential energy Π in eq. ( 50) is defined similarly as the Griffith functional for cracks in the variational theory of fracture [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF][START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF][START_REF] Karma | Phase-field model of mode III dynamic fracture[END_REF][START_REF] Sargado | High-accuracy phase-field models for brittle fracture based on a new family of degradation functions[END_REF], which depends on the material fracture properties and the surface Hausdorff measure. In contrast with the standard phase-field approaches to fracture, however, the discrete functional term W dis on S λ κ is not fully regularized over the subdomains Ω λ , so that the second integral in eq. ( 50) over S λ κ persists in the present formulation.

The energy release rate with respect to the shape perturbation of the dislocation loops is also formulated with respect to the variation of the potential energy Π, as follows

Π = Ŵbulk + Ŵdis -Ŵext , (51)
which also requires the rates of integrals defined on varying volume and surface domains, according to eq. ( 50). Using the formalism developed by [START_REF] Destuynder | Sur une interprétation mathématique de l'intégrale de Rice en Théorie de la rupture fragile[END_REF] for crack shape sensitivity and used here to compute the driving forces on dislocation loops, the shape of the body Ω is assumed to depend on a fictitious time τ through a continuum kinematics-type Lagrangian description. The unperturbed reference configuration Ω is also associated with τ = 0, while the motion of the body from Ω to Ω τ is represented by the following transformation mapping Φ(x, τ ), i.e.

∀ x ∈ Ω → x τ = Φ(x, τ ) ∈ Ω τ , (52)
where x and x τ are the position vectors of a material point in the reference and perturbed configurations, respectively, and Φ(x, 0) = x. In the neighborhood of the initial time, all time derivatives are determined at τ = 0, such that the first-order effect of infinitesimal perturbations of Ω ≡ Ω τ is considered only, thus ignoring higher-order terms. The corresponding initial transformation velocity θ(x) is also defined by

θ(x) = ∂Φ(x, τ ) ∂τ τ =0 , ( 53
)
while the dislocation shape sensitivity leads to specific requirements for both the morphing transformation Φ(x, τ ) and the initial velocity θ(x), as follows

(i) Ω τ ⊂ Ω ,
(ii) For τ ≥ 0, the Jacobian of Φ(x, τ ) must be strictly positive , (iii) The support of θ(x) is compact and localized closely to each dislocation front ∂S λ κ , so that θ(x) ≡ θ λ κ (x) ,

(iv) The prescribed tractions t N (x) are fixed and defined outside of the support of θ λ κ (x) ,

such that θ λ κ (x) represents the virtual growth of the dislocation loops S λ κ in the present infinitesimal extension approach. The Lagrangian derivative of a given field quantity f (x, τ ) is also defined by

f (x, τ ) = df (x, τ ) dτ τ = 0 = lim τ →0 f (x τ , τ ) -f (x, 0) τ = ∂f (x, τ ) ∂τ τ = 0 + f ,i (x, τ ) θ i (x) , (55)
and the corresponding material gradient of f (x, τ ) is given by

f ,i (x, τ ) = f,i (x, τ ) -f ,j (x, τ ) θ j,i (x) , (56)
with i = {1, 2, 3}. Furthermore, the derivatives of a generic volume and surface integrals read

FΩ = d dτ ˆΩτ f (x, τ ) dΩ = ˆΩ f (x, τ ) + f (x, τ ) θ i,i (x) dΩ = ˆΩ f (x, τ ) dΩ + ˆΣ f (x, τ ) θ n (x) dΩ FΣ = d dτ ˆΣτ f (x, τ ) dΣ = ˆΣ f (x, τ ) + f (x, τ ) θ i,j (x) δ ij -n Σ i n Σ j dΣ , (57)
respectively, as derived by [START_REF] Petryk | Time derivatives of integrals and functionals defined on varying volume and surface domains[END_REF], where δ ij is the Kronecker symbol, and θ n (x) = θ k (x) n Σ k , with n Σ k the unit normal to the surface differential element dΣ. By postulating the stationarity of Π, the variational relation in eq. ( 51) yields

Π = 0 ⇒ Ŵdis = -Ŵbulk + Ŵext , (58)
as a result of the present quasi-static morphing transformation Φ(x, τ ). On the other hand, the local Peach-Koehler force is defined by the line-integral requirement of the work on the displacement along the curved dislocation loop fronts ∂S λ κ , so that

N ∑ λ=1 K λ ∑ κ=1 ˆ∂S λκ F λ κ k (x) θ λ κ k (x) dl = -Ŵdis , (59)
exhibiting the fundamental relation between the driving forces on the dislocation loops and the corresponding energy release rate. According to eqs. ( 55), ( 56) and ( 57), the derivative of the elastic energy term W bulk with respect to τ in eq. ( 58) is given by

Ŵbulk = N ∑ λ=1 K λ ∑ κ=1 ˆΩλ 1 2 c λ kjml u λ m,l (x) u λ k,j (x) θ λ κ i,i (x) -c λ kjml u λ m,l (x) u λ k,i (x) θ λ κ i,j (x) + c λ kjml u λ m,l (x) ûλ k,j (x) dΩ , (60)
while the corresponding derivative of the external work term

W ext reads Ŵext = N ∑ λ=1 K λ ∑ κ=1 ˆ∂Ω λ ∩∂Ω t t N k (x) ûλ k (x) + t N k (x) u λ k (x) θ λ κ i,i (x) -θ λ κ i,j (x) n ext j n ext i dΣ , (61)
which completes the right-hand side of eq. ( 58). Since the derivative of the displacement field ûλ k (x) belongs to U 0 , the variational formulation of the present boundary-value problem in eq. ( 49) yields

N ∑ λ=1 ˆΩλ c λ kjml u λ m,l (x) ûλ k,j (x) dΩ = N ∑ λ=1 ˆ∂Ω λ ∩∂Ω t t N k (x) ûλ k (x) dΣ , (62)
so that, substituting eqs. ( 59), ( 60) and (61) in eq. ( 58) with eq. ( 62) and using (iv) in eq. ( 54), the total energy conservation law becomes

- N ∑ λ=1 K λ ∑ κ=1 ˆ∂S λκ F λ κ k (x) θ λ κ k (x) dl = N ∑ λ=1 K λ ∑ κ=1 ˆΩλ -1 2 c λ kjml u λ m,l (x) u λ k,j (x) θ λ κ i,i (x) + c λ kjml u λ m,l (x) u λ k,i (x) θ λ κ i,j (x) dΩ , (63)
which is widely used to compute the energy release rate for moving discontinuities, particularly successful in linear elastic fracture mechanics. Because eq. ( 63) is defined for any initial velocity field θ λ κ (x) that satisfies the four properties listed in eqs. ( 54), the standard finite element method enriched by specific adaptive remeshing techniques can be applied to compute the Peach-Koehler force per line length.

Finite element approximation of the non-singular elastic field solutions

This section deals with the discretization strategy adopted for the resolution of the variational relation in eqs. ( 49) and ( 63). The singularity of the elastic fields introduced by the local displacement jumps across the dislocation loops requires the use of a suitable approximation space [START_REF] Rajaram | Application of domain integral methods using tetrahedral elements to the determination of stress intensity factors[END_REF][START_REF] Wang | 3-D local mesh refinement XFEM with variable-node hexahedron elements for extraction of stress intensity factors of straight and curved planar cracks[END_REF]. Instead of employing the extended finite element method for dislocation loops [START_REF] Belytschko | On XFEM applications to dislocations and interfaces[END_REF], Gracie et al., 2007, 2008), the proposed technique based on adaptive meshing procedures becomes efficient for three-dimensional calculations of dislocation growth and allows i) to quasicontinuously describe general arbitrarily-shaped dislocation fronts by mixed-element meshes with a high accuracy with respect to the required boundary conditions as well as ii) to provide precise numerical elastic solutions of the mechanical equilibrium problem. In the following sections, the finite element nodes are denoted by the capital letters and are attached either to the superscripts or to the subscripts of the corresponding nodal quantities to avoid cumbersome notation.

Numerical discretization of the dislocation boundary-value problem solutions

As depicted in Fig. (2a), each dislocation loop S λ κ is discretized by a surface mesh composed of triangular elements from a standard Delaunay triangulation in the interiors of the dislocations, while quadrangular elements near the dislocation cores ∂S λ κ are created automatically to reduce the computational cost of the suitable approximation of the elastic field solutions. A conforming mixed-element mesh that contains tetrahedral, hexahedral, prismatic and pyramidal volume elements is subsequently built on the triangles and quadrilaterals of the discrete dislocations loops. Thus, a discretization scheme using three-dimensional isoparametric elements with in practice linear (but possibly quadratic) basis functions is illustrated in Fig. (2a) and is formulated as follows (i) The dislocation fronts ∂S λ κ are discretized into a set of ordered segments defined by a sequence of N λ κ vertices, such that ∂S λ κ ≡ [P λ κ I , P λ κ I+1 ], 1 ≤ I < N λ κ , for which P λ κ 1 = P λ κ N λκ for closed dislocation loops, or, ∂S λ κ ∩ ∂Ω = {P λ κ 1 , P λ κ N λκ } are two endpoints that belong to the external boundaries for dislocation loops interesting with the free surfaces;

(ii) The dislocated surfaces S λ κ , which are surrounded by the points {P λ κ I }, are characterized by double mesh entities formed by quadrangular and triangular elements with double interior nodes {M λ κ + J } and {M λ κ - J } on the upper and lower side of the dislocated surface, with 1 ≤ J ≤ N λ κ , respectively. The displacement jump is also prescribed to any couple of superimposed points, i.e., {M

λ κ + J } ≡ {M λ κ - J }, such that u λ κ k (M λ κ + J ) -u λ κ k (M λ κ - J ) = b λ κ k χ * (M λ κ + J ) , (64)
as stipulated by eq. ( 22). In contrast with the internal grain/heterophase boundaries, if a dislocation loop boundary reaches the external boundary, thus 64), so that no regularization is applied when a dislocation loop boundary reaches the external boundary in accordance with (i), leaving exactly a surface step of b λ κ on the free surfaces;

∀ M λ κ + J ∈ ∂Ω : χ * (M λ κ + J ) = 1 in eq. (
(iii) Around the discretized fronts ∂S λ κ , a few regular and concentric rings of elements of constant thickness r h = r 0 / m are constructed, where m is a user-specified integer and r 0 the regularization radius of the dislocation cores, as displayed by the spider-web mesh in red in Fig. (2a). The first layer of the spider-web mesh is also described by prismatic elements, while the last hybrid layer consists of five-node pyramidal and tetrahedral elements that ensure the conforming connection to the rest of the unstructured tetrahedral mesh domain Ω λ , as shown in blue in Fig. (2a).

In the fully conformal finite element mesh around the dislocation loops as well as the grain boundaries and multiple different heterophase domains, the displacement field at any point x is approximated by

u k (x) = N n ∑ I=1 N I (x) U I k , (65)
where N n is the total number of nodes in the computational mesh, N I (x) are the local basis functions associated with the node I, and U I k = u k (x I ) are the nodal displacements in each coordinate-direction with respect to the nodal coordinates x I . The external Dirichlet boundary conditions in eq. ( 46) as well as the displacement boundary condition in eq. ( 64) are therefore written with respect to the discretization form for each dislocation loop S λ κ as follows

U J k = u D k (x J ) , ∀ x J ∈ ∂Ω u , U J k -U I k = b λ κ k χ * (x I ) , ∀ (x I , x J ) ∈ S λ κ × S λ κ x I I<J = x J , (66)
respectively, with a global numbering of the mesh vertices ensuring that any node on the lower side of the dislocation loop surface is ranked lower than the corresponding node on the upper side. The integral equation of the dislocation boundary-value problem in eq. ( 49) is transformed into a global linear system of equilibrium equations, i.e.

K IJ km U J k = F I m with :        F I m = ˆ∂Ω N I (x) t N m (x) dΣ K IJ km = N ∑ λ=1 ˆΩλ N I ,j (x) c λ kjml N J ,l (x) dΩ , (67)
where F I m represents the vector components of the generalized forces, and K IJ km are the (3N n ) × (3N n ) components of the stiffness matrix of the three-dimensional discretized variational formulation of the dislocation boundary-value problem, also modified to take into account the prescribed displacement conditions in eq. ( 66), while N I ,j (x) is the gradient of the I th nodal basis function with respect to the j th -direction. The numerical resolution of the system in eq. ( 67) allows to determine the set of nodal displacement values U J k and to reveal the corresponding state of stress fields, as commonly carried out in standard finite element analysis.

Computational Peach-Koeler force using a finite element formulation

According to the properties of the transformation velocity θ λ κ in eq. ( 54), specially (iii), the integral equation of energy conservation in eq. ( 63) can be split with respect to each dislocation loop S λ κ , such that the discrete approximation of the individual Peach-Koehler force for each front ∂S λ κ can be derived, thus omitting the sum over κ. By virtue of eq. ( 54), because the virtual extension field θ λ κ (x) is defined on a compact and non-zero support close to the dislocation fronts, a one-to-one projection P, which maps any point x near enough ∂S λ κ to the corresponding point P λ κ (s) on the dislocation fronts ∂S λ κ , can be introduced as follows

∀ x ∈ Ω λ θ λ κ k (x) 0 , P : x → (s, r, ψ) ∂S λκ , (68)
whose Euclidean radial distance r between x and P λ κ (s) on ∂S λ κ is the smallest. The latter point is also parametrized in a local curvilinear coordinate system (s, r, ψ), where s is the curvilinear abscissa, and ψ is the polar angle between x -P λ κ (s) and m λ κ (s) = ξ λ κ (s) × n λ κ , which represents the unit and local outward-pointing direction vector of dislocation front onto the glide plane, with ξ λ κ (s) the local tangent vector onto the same plane, as illustrated in Fig. (2a). Thus, an explicit function q(r) that relates the virtual growth θ λ κ (x) to θ λ κ (s) is formally given by

θ λ κ k (x) = q(r) θ λ κ k (s) , with q(r) =    1 , if r ≤ r a , 0 , if r ≥ r b , 2 (r -r a )/(r b -r a ) 3 -3 (r -r a )/(r b -r a ) 2 + 1 , otherwise , (69)
where r a and r b are two parameters that fulfill the conditions in eq. ( 54) and are defined with respect to the representative local size of the mesh elements near the dislocation fronts, with r 0 ≤ r a < r b .

In the same fashion than the previous computation of the nodal displacement field in eq. ( 67), the discrete values of the Peach-Koehler forces can therefore be evaluated by the finite element equilibrium resolution. Since the mapping transformation P can be applied in the whole support of the virtual extension field θ λ κ (x) and because of the in-plane propagation of the present dislocation loops, the components of both quantities θ λ κ k (x) and F λ κ k (x) are therefore in plane and collinear to the vector m λ κ k (s). Thus, according to eq. ( 69), it follows that

P(x) = (s, r, ψ) ∂S λκ : ∀ x ∈ ∂S λ κ : F λ κ k (x) = F λ κ (s) m λ κ k (s) ∀ x ∈ Ω λ θ λ κ k (x) 0 : θ λ κ k (x) = q(r) θ λ κ k (s) = q(r) θ λ κ (s) m λ κ k (s) , (70)
where F λ κ (s) and θ λ κ (s) are the magnitude of the Peach-Koehler force and of the extension field of the dislocation front ∂S λ κ in the direction m λ κ (s), respectively. In the context of the conventional finite element technique, the scalar-valued quantities θ λ κ (s) and F λ κ (s) are interpolated along the dislocation fronts ∂S λ κ from the vertices {P λ κ I }, as follows

θ λ κ (s) = M λ κ I (s) Θ λ κ I F λ κ (s) = M λ κ I (s) F λ κ I , (71)
where M λ κ I (s) are the linear interpolation functions along ∂S λ κ , chosen for reasons of simplicity and stability, while Θ λ κ I and F λ κ I are the nodal values of the signed magnitude of the virtually extended dislocation loops and of the Peach-Koehler force along the discretized dislocation fronts, respectively. Since the variational formulation in eq. ( 63) is valid for any admissible field θ λ κ k (x) and a fortiori for any nodal restriction Θ λ κ I , the integral equation can also be discretized to obtain the following linear system, i.e.

S λ κ IJ F λ κ J = R λ κ I , (72)
where the N λ κ × N λ κ components of the mass matrix S λ κ IJ and the discretized second member components R λ κ I are defined by

               S λ κ IJ = ˆ∂S λκ M λ κ I (s) M λ κ J (s) dl R λ κ I = N ∑ λ=1 ˆΩλ -1 2 c λ khml u λ m,l (x) u λ k,h (x) δ ij + c λ kjml u λ m,l (x) u λ k,i (x) q ,j (r) M λ κ I (s) + q(r) M λ κ I ,j (s) m λ κ i (s) + q(r) M λ κ I (s) m λ κ i,j (s) dΩ , (73)
respectively, with 1 ≤ I ≤ N λ κ , 1 ≤ J ≤ N λ κ , and ,j the gradient of the corresponding explicit function with respect to the j th -direction. The nodal values of the configurational Peach-Koehler magnitudes F λ κ J along the dislocation fronts ∂S λ κ are therefore determined on all boundary nodes by inverting eq. ( 72) with eqs. ( 73) for any arbitrarily-shaped dislocation loops with arbitrary characters.

Dislocation propagation and coalescence with adaptive remeshing procedures

From the previous discrete values of the driving forces amplitudes F λ κ J applied to the dislocation loops with arbitrary shapes, the following empirical mobility law is used to describe the relation between force and velocity in the slip planes [START_REF] Hirth | Theory of dislocations[END_REF], i.e.

V λ κ J = B -1 F λ κ J m λ κ (s J ) , ( 74
)
where B is a viscous friction coefficient, while s J is the curvilinear coordinate of the discrete dislocation front node P λ κ J , so that V λ κ J is the nodal velocity vector used to describe the corresponding change in the dislocation topologies. By means of an explicit variable time increment, the mobility law is also applied to each node of all discrete dislocation loops S λ κ as long as the amplitude of the velocity falls below a given velocity threshold value, so that the dislocation loops can evolve towards equilibrium states until the convergence is reached.

The dislocation loop topologies are updated through a prescribed time step ∆t, such that each node P λ κ J of the dislocation loop fronts is explicitly moved by a local displacement ∆t V λ κ J . In the present finite element framework, the current time step is chosen in order to provide a maximum propagation of the fastest front node comparable to the local mesh size. Depending on the sign of the nodal Peach-Koehler magnitude F λ κ J , the dislocation loop expands in the direction m λ κ (s J ) for positive values and contracts for negative values, for which mesh-based surface deformation and morphing techniques are used to generate the newly morphed meshes using appropriate intermediate triangulated sections, as illustrated in Fig. (2b). Because the construction of the newly morphed triangulated meshes to describe the expansion and shrinkage of dislocation loops can lead to significant topological modifications, combining with inevitable short-range reactions between multiple dislocation loops as well as between dislocation loops and external and internal boundaries in realistic microstructures, special precautions are therefore pursued to capture the propagation and coalescence of multiple dislocation loops. In particular, two subsequent adaptive procedures can be performed under specific dislocation-based scenarios, as illustrated by the application examples in the next section, i.e.

(i) Short-range core reaction between intersecting dislocations is required to describe the collinear interactions and annihilation of two dislocation loops with opposite Burgers vector gliding in the intersecting slip planes. The numerical finite element treatment for the collinear interaction involves local topological criteria for coalescing surface generation that are based on the adaptive mesh sizes as well as on the element aspect ratio constraints relative to the active crystal slip systems, while meeting the displacement discontinuities at the double nodes on the dislocation surfaces; (ii) After the explicit coalescence of planar and non-planar dislocation surfaces, if the newly-triangulated surface intersects with other dislocation loop surfaces or domain boundaries, an efficient surface-surface intersection algorithm is applied to generate a conforming surface triangulation with respect to the positions of the discrete dislocations and of sub-domain boundaries of Ω.

Both procedures (i) and (ii) are fully automated at each propagation step and are topologically consistent with the threedimensional discretization scheme defined in section 3.3.1, especially with the three-dimensional conforming mixed-element mesh procedure that generates tetrahedral, hexahedral, prismatic and pyramidal volume elements on the triangles and quadrilaterals of each discrete dislocations loop present in the computational volume.

Application examples to various dislocation problems

Various application examples are studied, from a single circular shear dislocation loop embedded in a homogeneous and isotropic elastic material for validation purposes, to the propagation of the same dislocation loop in a polycrystalline micropillar interacting with grain boundaries and an infinitely hard precipitate, to the collinear annihilation of intersecting perfect dislocations in a voided face-centered cubic crystal. The present calculations are performed on copper, unless otherwise stated for heterophase material components. Copper is moderately anisotropic, where the elastic constants are defined by c 11 = 168.4 GPa, c 12 = 121.4 GPa, and c 44 = 75.4 GPa, while the magnitude of the Burgers vector is b = 0.25 nm. Special attention is devoted to determining numerically accurate elastic field solutions and singularity-free Peach-Koehler forces with respect to the theoretical analysis, without requiring time intensive computations in three dimensions, as quantified in Tab. (1) for the numerical application examples with evolving dislocation characteristics. The present finite element simulations are performed using a standard 12 core workstation with 32 GB of memory, for which roughly half of the computational time is dedicated to the meshing events, the same amount to the finite element solver, while the CPU time of the Peach-Koehler force post-processing is marginal. 4.1. Validation benchmarks by comparing the isotropic stress fields with existing analytical solutions

Figure (

3) shows the non-zero stress components produced by a circular shear dislocation loop in a homogeneous and isotropic linear elastic material, as illustrated in the schematics, for validation and comparison. The results are obtained by the present theoretical (solid lines) and numerical finite element (red crosses) solutions, as derived in sections 2 and 3, respectively, while the blue crosses illustrate the existing closed-form analytical solutions [START_REF] Hirth | Theory of dislocations[END_REF]. The isotropic material parameters of copper are determined using the Voigt averaging procedure, and the core-spreading parameter is chosen to be sufficiently small r 0 = b/20, for a consistent comparison with the analytics. In order to reveal the most significant discrepancies between the theoretical and numerical solutions, the results are plotted in the close vicinity of the dislocation cores at z = z s + b and z = z s + 10b, in Fig. (3a) and (b), respectively, for nanoscale shear dislocation loops located at z = z s = H/2, with a radius R = 40b. The stress profiles, which decrease rapidly with distance to the dislocation loop as -2 , are depicted along the x 1 -and x 2 -directions which are parallel and orthogonal to the Burgers vector b = b [1 0 0] t , respectively. The dimensions of the computational domain are L 2 × H = (320b) 3 , and the mesh is composed of ∼ 428k degrees of freedom. The present theoretical solutions are strictly identical to the existing analytical formulas, and are in excellent agreement with the finite element formulation at very short and long ranges. The oscillating characteristics of the stress fields at the aplomb of the dislocation cores are also very well reproduced by the numerical approach.

Shear dislocation loops in 3D homogeneous and heterogeneous anisotropic elastic materials

Using the previous calculation setup with the fully anisotropic elastic constants of copper, Fig. ( 4) shows an equally valuable comparison between the present theoretical framework and the finite element computations, for all components of the in-plane and normal stress field. In general, the numerical results are thus able to reproduce the discrepancies induced by the isotropic approximation of the elastic properties. On the other hand, Fig. ( 5) illustrates the anisotropic elastic field solutions produced by an elliptical shear dislocation loop with semi-axis ratio a 1 /a 2 = 2 with 2a 1 = 80b, embedded in the lower solid of a freestanding anisotropic bimaterial of total thickness H = 2a 1 . The horizontal dislocation loop with b = b [1 0 0] t is located in the 16 center of the heterogeneous structure, as indicated by the red double-headed arrows, while the heterophase interface is placed in the mid-plane between the upper free surface and the dislocated surface. The elastic constant c 11 of the upper solid has been fictitiously increased by a factor of four, involving a significant asymmetric distribution of the stress state with respect to the median surface at z = z s = H/2, as shown by the three-dimensional isosurface-type contours of the von Mises equivalent stress value of σ vM = 0.5 GPa for both theoretical and numerical analysis. In particular, the latter is conducted using two coarse and fine meshes, which consist of ∼ 152k and ∼ 428k degrees of freedom, respectively. Moreover, the corresponding component of the normal stress field σ 33 is plotted on each contour, for further comparison at long distances. The results are in close agreement and illustrate the presence of a non-zero von Mises stress state at the upper free surface, which is also emphasized by the fine mixed-element mesh, although the traction-free boundary condition at the upper surface is satisfied by both analyses. The von Mises stress field discontinuity is fairly consistent in magnitude for both calculations, while the continuity of the normal stress field across the heterophase interface is also fulfilled, as required by the boundary conditions for both theoretical and numerical calculations. In addition, the stress concentrations are reasonably well captured by the finite element simulations.

x1 / L -1.5 -1 0 1 1.5 3 -1.5 0 1.5 -3 x2 / L -1.5 -1 0 1 1.5 3 -1.5 0 1.5 -1.5 x1 / L -1.5 -1 0 1 1.5 1.5 Stress at z = zS + 10b (GPa) -1 0 x1 / L -1.5 -1 0 1 1.5 x1 / L -1.5 -1 0 1 1.5 x2 / L -1.5 -1 0 1 1.5 Stress at z = zS + 10b (GPa) x1 / L -1.5 -1 0 1 1.5 x2 / L -1.5 -1 0 1 1.5 -0.5 0.5 1 -1.
4.1.2. Non-singular Peach-Koehler forces on complex-shaped shear and prismatic dislocation loops According to the theoretical analysis in section 2, the core-spreading procedure is performed by convoluting the discontinuity displacement condition at the surface of the dislocation loop with the Gaussian function in eq. ( 31). As illustrated in Fig. (6a) for the specific circular shear dislocation case of Figs. ( 3) and ( 4) with three values of r 0 = {1b, 2b, 5b}, the displacement jump is spread out and becomes continuously defined over the dislocation cores with different amplitudes, whereas the jump is characterized by a piecewise rectangular function in the classical theory of dislocation loops with compact cores , i.e., for r 0 = 0. The corresponding non-singular in-plane σ 11 (x 1 , 0, z s ) and shear σ 13 (x 1 , 0, z s ) stress field components are depicted in Fig. ( 6b) and (c), respectively, using the fully anisotropic elasticity theory (solid lines) and the isotropic elastic approximation (dashed lines). Interestingly, the isotropic elastic solutions illustrate a higher magnitude of in-plane stress than the anisotropic elastic solutions, while σ 13 (±R, 0, z s ) is zero the dislocation fronts. The latter is not defined at the dislocation cores in classical elasticity.

x1 / L -1.5 -1 0 1 1.5 (a) x1 / L -1.5 -1 0 1 1.5 x1 / L -1.5 -1 0 1 1.5 x2 / L -1.5 -1 0 1 1.5 x1 / L -1.5 -1 0 1 1.5 x2 / L -1.5 -1 0 1 1.5 -1.5 x1 / L -1.5 -1 0 1 1.5 1.5 -1 0 x1 / L -1.5 -1 0 1 1.5 x1 / L -1.5 -1 0 1 1.5 x2 / L -1.5 -1 0 1 1.5 x1 / L -1.5 -1 0 1 1.5 x2 / L -1.5 -1 0
The associated singularity-free magnitude of the Peach-Koehler forces, which are obtained by the theoretical solutions from eq. ( 40) and the numerical finite element analysis by solving eq. ( 72) with eqs. (73), is plotted against the polar angle θ in Fig. (7a) for the circular shear dislocation loop. In particular, θ ≡ 0 mod 180 • (θ ≡ 90 • mod 180 • ) corresponds to the point on the dislocation front where the local character has a pure screw (edge) component. In comparison with the fully anisotropic elasticity calculations (solid lines), the significant deviations produced by the isotropic elastic approximation (dashed lines) are shown. Theoretical and numerical calculations indicate that the magnitude of the driving force is lowest for the local edge characters, while the isotropic elasticity predicts maximum values. The excellent agreement between the two approaches extends to the prediction of the camel hump-shaped force profiles (e.g., black solid lines) on either side of the pure screw component ±30 • when using anisotropic elasticity with smaller values for the regularization parameter r 0 . Furthermore, the Peach-Koehler forces that are associated with the elliptical shear dislocation configuration embedded in the free-standing anisotropic elastic bimaterial in show that the local screw components along the major axis experience the largest driving force in magnitude. However, the finite element results reveal some oscillations around the pure screw characters, which do not appear in the theoretical solutions.

The singularity-free Peach-Koehler magnitudes for prismatic dislocation loops with complex butterfly-and skull-shaped fronts are presented in Figs. (8a) and (b), respectively, for r 0 = 5b. The theoretical (numerical) solutions are shown as solid lines (with symbols), while the corresponding driving forces are drawn in pink along the contours on the right-hand sides, with and without the two-dimensional shear stress σ 12 (x 1 , x 2 , z s ) maps in the background for further comparison. The signed magnitudes of the Peach-Koehler forces are plotted against the polar angle θ, for which θ = 0 • corresponds to the points M in the schematics. In general, the very good agreement in terms of stresses and forces in sign and magnitude is also demonstrated, although slight deviations in direction are noticeable when the local radius of curvature changes drastically in sign. These discrepancies are mainly due to the different core-spreading schemes that have been appropriately adopted for mathematical convenience in each of the theoretical and numerical formulations. 74) with B = 10 -5 Pa.s. The magnitude of the self-force is plotted against the curvilinear coordinate s, for which s = 0 corresponds to the points M in the six associated snapshots of Fig. (9c). In the latter figure, the von Mises equivalent stress field is shown during the dislocation shrinkage, which is computed numerically with a core-spreading parameter r 0 = 1b. According to Tab. (1), the mesh of the computational volume for the initial configuration #1 is composed of ∼ 126k degrees of freedom, while the number of degrees of freedom is reduced to ∼ 72k for the configuration #6, and the average discrete time step between each remeshing event is 6.8 × 10 -15 s.

In accordance with the results in Fig. ( 7a), the negative driving force experiences the smallest magnitude for local edge characters, such that the computational contour of the loop does not remain circular during collapse, and the evolution of the dislocated surface becomes nonlinear from configuration #5, as indicated by the deviation from the dotted line and already described by the phase field approach [START_REF] Rodney | Phase field models and dislocations[END_REF][START_REF] Ruffini | Three-dimensional phase-field model of dislocations for heterogeneous face-centered cubic crystal[END_REF]. However, the present plots show more than the phase field results, as the importance of anisotropic elasticity is emphasized by promoting the appearance of faceted dislocation contours, as similarly observed with sharp-cornered dislocation configurations in anisotropic iron [START_REF] Fitzgerald | Shape of prismatic dislocation loops in anisotropic α-Fe[END_REF][START_REF] Aubry | Equilibrium shape of dislocation shear loops in anisotropic α-Fe[END_REF]. Interestingly, the self-force increases in magnitude and tends roughly to constant values from configuration #3 as the dislocation loop narrows.

Because the long-range elastic interaction and the short-range core reaction between intersecting dislocations is crucial for gaining insight into dislocation-based strain hardening mechanisms in materials, the collinear annihilation of two dislocation loops expanding on intersecting slip planes, is reported in Fig. (10a). Thus, two perfect dislocation shear loops with Burgers vectors b

1 1 = b/ √ 2 [1 1 0] t and b 1 2 = b/ √ 2 [1 1
0] t are introduced into the (111) and ( 11 1) slip planes, respectively, under a macroscopic applied stress σ 11 = 0.08c 11 . In particular, a much more complex configuration than the standard situation is described herein when the collinear reaction interacts with a spherical cavity in anisotropic copper, as illustrated by #0 with ∼ 56k degrees of freedom under periodic boundary conditions. The finite element formulation is able to account for the completely infinite contrast between phases, which remains a numerical challenge for traditional fast Fourier transform solvers with an accurate convergence rate. During the propagation and the collinear annihilation of the two dislocation loops with an average discrete time step of 8.5 × 10 -15 s at the nanoscale, the snapshots show the magnitude of the von Mises equivalent stress field and of the Peach-Koehler forces along the dislocation fronts. Under the external stress state, the two shear dislocation loops also expand, interact with each other, and annihilate in place where the opposite line directions meet at the intersecting glide planes. The remarkable versatility of the present conforming mixed-element mesh therefore allows for collinear annihilation of the mesh edges during the short-range core reaction from the configuration #3, without resulting dislocation lines as well as residual stress fields when the core reaction between both dislocation loops occurs at the product junction. The cavity is therefore symmetrically sheared by the individual cut surfaces, leaving two slip steps at the free surfaces of the spherical void, as stipulated in section 3.3.1, and also illustrated in the deformation configuration space in Fig. (10b). The configuration #8 features ∼ 229k degrees of freedom, in which the dislocation contour in the (111) glide plane is no longer symmetric about the [ 1 1 2] t -axis due to the influence of the stress field caused by the spherical void on the propagation and coalescence of the dislocation loops.

Dislocation dynamics of the Orowan precipitate bypass mechanism in a polycrystalline micropillar

Figure (11a) illustrates a large-scale three-dimensional finite element computation that cannot, to the knowledge of the authors, be achieved by existing numerical approaches in the broader literature. An anisotropic copper polycrystalline micropillar with 80 grains is automatically generated from the intersection of a cubical Voronoi tessellation with a representative pillar specimen, in which a shear dislocation loop with a Burgers vector b = b/ √ 2 1 1 0 t glides in the (111) slip plane of a specific host grain. The latter lies outside the microstructure, so that the outer grain boundary corresponds to the free surface of the computational sample. A high compressive strain of 7.1% is applied and maintained constant on one external face of the copper micropillar, for which the final configuration mesh is composed of ∼ 1007k degrees of freedom. The corresponding animation of the Orowan precipitate bypass mechanism is referred to as "Orowan bypass mechanism in a micropillar", computed in less than 20 hours with 291 adaptive remeshing events with an average discrete time step of 4.6 × 10 -13 s, according to Tab. (1).

Concluding remarks

The present work achieves a twofold objective by addressing the long-standing problem of arbitrarily-shaped dislocation loops in three-dimensional heterogeneous materials from a theoretical and numerical perspective. First, elastic field solutions for shear and prismatic dislocation loops in heterogeneous anisotropic multilayered structures are formulated in terms of the mathematically elegant Stroh formalism combining with the biperiodic Fourier-transform and dual variable and position techniques. A suitable Gauss-Legendre quadrature enriched by the Green integral formula is employed to cope with the dislocation loops with complicated simply-connected contours. In particular, the singularity-free Peach-Koehler force is derived using a core-spreading operation to regularize the prescribed displacement discontinuity, which is considered as a standard boundary condition of the dislocation boundary-value problem. On the other hand, a monolithic finite element formulation with advanced and flexible adaptive remeshing procedures is subsequently introduced to capture the configurational driving forces through dissipative energy considerations and to describe the dynamics of dislocation loops in realistic microstructures and structures. Excellent agreement between theoretical and numerical analyses is illustrated from simple circular shear dislocation loops to complicated butterfly-shaped prismatic dislocation loops in linear homogeneous isotropic solids and anisotropic elastic bimaterials. The numerous comparisons on complex configurations of dislocation loops also serve as a validation basis and benchmarks for dealing with sophisticated problems with evolving dislocation loops in three dimensions. The remarkable versatility and robustness of the finite element method allows for the consideration of the short-range core reaction between intersecting dislocation loops in interaction with precipitates (cavity and rigid inclusion), as well as the Orowan dislocation-precipitate bypass mechanism in a compressed micropillar of polycrystalline copper. The computationally efficient and robust adaptive mesh generation procedures for dislocation propagation, interaction, and coalescence are also capable of handling multiscale applications. The dislocation propagation by bowing around the inclusion in the polycrystalline sample exhibits three orders of magnitude in representative size scale between the micropillar and the nano-precipitate.

At first glance and in the current form, the finite-element framework should be considered as a computational tool to carry out calculations with several types of discontinuities, such as grain boundaries, free surfaces, dislocation loops and cracks, in multiphase finite material structures. The main interesting feature of the approach is to unify these discontinuities into a single finite-element entity to revisit the fundamental problems concerning the interactions between dislocation loops and cracks, in particular the emission of dislocations from crack fronts in three dimensions, as well as the interactions between dislocations and stress concentrations at grain boundaries and heterophase interfaces, especially the nucleation and emission of dislocation loops from the internal material boundaries. Although the computational approach undoubtedly opens many perspectives, also with close links to experiments, some extensions can be introduced. A current limitation is related to the use of a single regularization rule at the dislocation fronts, whether the dislocation loops are located in the core of the grains or near the internal interfaces. A more physics-based rule could be provided to offer a better description of the short-range elastic fields close to the grain boundaries to analyze the transmission of dislocation loops into neighboring grains, thus overcoming the current impermeability conditions. Furthermore, although the current simulations are performed on a workstation, the numerical framework could benefit from the robust iterative and domain decomposition solvers to handle the discretization of several tens of millions of unknowns [START_REF] Bovet | Two-level adaptation for Adaptive Multipreconditioned FETI[END_REF]. By the use of a parallel mesh generation algorithm for robust domain decomposition techniques, high-performance calculations with a hundred dislocation loops are anticipated to characterize standard dislocation microstructures with typical densities of 10 12 /10 14 m -2 in the 1-to-100 micrometer mesoscale range. Finite element calculations with hundreds of millions of degrees of freedom are therefore expected to achieve such numerical experiments for multiple dislocation loops in three-dimensional material structures. These subsequent boundary-value problems should be accompanied by consideration of additional dislocation junctions, such as the Lomer-Cottrell lock, the Hirth lock and the glissile junction as well as the implementation of the dislocation cross-slip mechanism and energetics, which are left for future investigations. In an extrapolation scenario, computations of several thousand dislocation loops on supercomputers could be carried out with the aim of better understanding dislocation-based strain hardening mechanisms in realistic structures at the macroscale. In accordance with the boundary conditions for impenetrable precipitates, the finite element mesh of the shear dislocation loop conforms to the precipitate and all intersection nodes that belong to the two defects are defined by a zero displacement jump. An Orowan dislocation loop is thus left around the infinitely strong precipitate.

  Figure1: (a) Ideal multilayered materials embedded by a single horizontal dislocation loop of arbitrary shape and character. Within linear elasticity theory in heterogeneous systems, the layers are treated as dissimilar and anisotropic elastic plates, defined by arbitrary thicknesses and parallel perfectly bounded internal interfaces. (b) Polycrystalline samples with arbitrary-oriented grains subjected to complex loading conditions, within which a discrete dislocation loop interacts with the internal grain boundaries as well as a heterogeneous precipitate. In such a realistic material configuration, the dislocation dynamics are treated by use of the finite element method combined with adaptive remeshing procedures.

Figure 2 :

 2 Figure 2: (a) The finite-element representation of arbitrarily-shaped dislocation loops with simply-connected contours ∂S λκ in three-dimensional adaptive mixed-element meshes using flexible conforming refinement. The displacement discontinuity is described by quadrangular and triangular elements with double nodes to prescribe the Burgers vectors attached to the dislocation loops. A spider-web mesh in red is used to define the regularization operation of the dislocation core around the fronts that are represented by a set of external vertices connected by structured boundary layers. (b) The finite element approach is capable of computing the singularity-free Peach-Koehler force at each external node and accordingly modeling the complex dislocation dynamics simulations based on local remeshing techniques.

Figure 3 :

 3 Figure3: Non-zero stress profiles obtained by the present theoretical (solid lines) and numerical finite element (red crosses) solutions and compared with analytical solutions (blue crosses) for a circular shear dislocation loop in a linear homogeneous and isotropic elastic material. The dislocation loop is located at z s , and the stress distributions are displayed at (a) z = z s + b and (b) z = z s + 10b, along the x 1 -and x 2 -directions that are parallel and orthogonal to the Burgers vector b = b [1 0 0] t , respectively.

Figure 4 :

 4 Figure 4: Non-zero stress profiles obtained by the present theoretical (solid lines) and numerical finite element (symbols) solutions for a circular shear dislocation loop in a linear homogeneous and anisotropic elastic material. The dislocation loop is located at z s , and the stress distributions are displayed at (a) z = z s + b and (b) z = z s + 10b, along the x 1 -and x 2 -directions that are parallel and orthogonal to the Burgers vector b = b [1 0 0] t , respectively.

Figure 5 :

 5 Figure5: Three-dimensional isosurface-type contours of the von Mises equivalent stress value of σ vM = 0.5 GPa produced by an elliptical shear dislocation loop in a free-standing anisotropic bimaterial with a ratio of semi-axes a 1 /a 2 = 2, as illustrated in the schematic in (a). The horizontal dislocation loop with b = b [1 0 0] t is located at the center of the bimaterial, as indicated by the red double-headed arrows, while the heterophase interface is placed in the midway plane between the upper free surface and the dislocation loop. The stress field component σ 33 is plotted on each contour in the (a) (x 1 , x 3 ), (b) (x 2 , x 3 ), and (x 1 , x 2 ) coordinate planes with respect to the theoretical and numerical finite element solutions. The latter solution is obtained by use of a coarse and fine mesh for comparison purposes.

Figure 6 :Figure 7 :

 67 Figure6: (a) The core-spreading displacement jump across a circular shear dislocation loop along the x 1 -direction, with three values of the dislocation core parameter r 0 . The corresponding singularity-free in-plane (b) normal σ 11 (x 1 , 0, z s ) and (c) shear σ 13 (x 1 , 0, z s ) stress field components obtained by the fully anisotropic elasticity theory (solid lines) and the isotropic elastic approximation (dashed lines). The plots are depicted exactly at z = z s .

Figure 9 :

 9 Figure 9: (a) Evolution of the magnitude of the Peach-Koehler forces on a shear dislocation loop with b = b [1 0 0] t at different time steps under the action of the self-stress contribution, only. The plots are displayed with respect to the curvilinear coordinate s. (b) Normalized surface of the shrinking shear loop as a function of time. (c) Six snapshots of the dislocation shrinkage with the associated change in the von Mises equivalent stress field, for which M corresponds to the magnitude where s = 0 in (a).

Figure 10 :

 10 Figure 10: (a) Collinear reaction between two dislocation shear loops with Burgers vectors b 11 = b/ √ 2 [1 1 0] t and b 1 2 = b/ √ 2 [110] t in the (111) and ( 11 1) slip planes, respectively, in interaction with of a spherical void in copper. The conforming mixed-element mesh allows to account for the collapsing mesh edges at the collinear annihilation, while the standard finite element is able to model the fully infinite contrast between the phases under the macroscopic applied stress σ 11 = 0.08c 11 . The snapshots show the corresponding amplitude of the von Mises equivalent stress field as well as the magnitude of the Peach-Koehler forces. (b) In the strain configuration space, the propagation of the shear dislocation loops produces two slip steps, which are left symmetrically at the free surfaces of the spherical cavity. Outside the free surfaces, neither resulting dislocation lines nor residual stress fields are created when the core reaction occurs due to the collinear annihilation of intersecting dislocation loops.

  Figure 11: (a) A shear dislocation loop around the bypassed heterogeneous precipitate in a monotonically deformed copper polycrystalline micropillar. The internal interfaces, such as the grain boundaries and the heterophase precipitate interfaces, are regarded as impenetrable barriers to the dislocation motion. The initial mesh consists of ∼ 193k degrees of freedom, within which the dislocation loop with a Burgers vector b 1 = b/ √ 2 1 1 0 t glides into the dense (111) slip plane, interacts with the long-range internal stresses, and intersects on a long distance the free surface of the sample. (b) The snapshots illustrate the advanced planar propagation of a dislocation loop in the polycrystalline micropillar by cutting completely the host grain in a realistic three-dimensional microstructure.The three-dimensional distribution of σ 33 (x 1 , x 2 , x 3 ) is mapped during the planar propagation, dislocation/precipitate interaction and self-coalescence.

Figure 12 :

 12 Figure 12: The enlarged inset #8 illustrates the regularization procedure on specific finite element nodes along the expanded dislocation front sliding in the host grain and locally interacting with the internal boundaries of the heterogeneous precipitate. The propagation step corresponds to the configuration #8 in Fig. (11b).In accordance with the boundary conditions for impenetrable precipitates, the finite element mesh of the shear dislocation loop conforms to the precipitate and all intersection nodes that belong to the two defects are defined by a zero displacement jump. An Orowan dislocation loop is thus left around the infinitely strong precipitate.

Figure 13 :

 13 Figure 13: From the initial dislocation loop embedded in a given grain of the polycrystalline copper micropillar with ∼ 193k degrees of freedom in (a) to the various propagation steps followed by the shear dislocation loop in (b), thus leaving residual dislocation edges around the bypassed heterophase precipitate. The final computational mesh involves ∼ 1007k degrees of freedom.

Table 1 :

 1 Computation time for the finite element application examples with evolving dislocation features, where the underlined quantities indicate average values.

	Application examples	Remeshing steps	Simulation discrete time step ∆t (10 -12 s)	Degrees of freedom number per step	Wall-clock time per step	Total computation wall-clock time
	Circular loop shrinkage (section 4.1.3)	39	0.013 > 0.0068 > 0.0030	126k > 82k > 49k	105s > 38s > 27s	24m30s
	Collinear interaction (section 4.1.3)	99	0.00049 < 0.0085 < 0.15	56k < 159k < 229k	65s < 134s < 171s	3h41m
	Polycrystalline micropillar (section 4.1.4)	291	0.018 < 0.46 < 51.1	193k < 410k < 1007k 83s < 238s < 668s	19h14m

specimen, while the opposite face is blocked, as displayed in Fig. (11a). At the grain scale, the Orowan bypass mechanism is described by the presence of the infinitely stiff, also elastically mismatched precipitate of arbitrary shape, for which the elastic constants are fictitiously multiplied by a factor of ten, with impenetrable boundaries and without consideration of cross-slip events. The internal grain boundaries are also considered as impenetrable barriers to dislocation motion, so that the dislocation loop is strictly confined to the host grain. The initial number of degrees of freedom associated with the full mesh of The quadrilateral plane elements of the regular spider-web mesh generated at the shear dislocation core are clearly distinguishable inside the dislocation loop near the mobile front (in red) that propagates in the (111) glide plane. For the corresponding finite element nodes, the regularization procedure defined by eq. ( 64) with eq. ( 47) is also applied, such that the displacement discontinuity is prescribed from 0 to b by means of double nodes, as evidenced by the two adaptive meshes in black and white. For the companion nodes that intersect the impenetrable boundaries of the heterogeneous precipitate, the generation of the regular spider-web mesh is not performed at the internal interfaces for simplicity, but the set of conforming dislocation nodes along the precipitate is defined by a zero displacement jump. Thus, an Orowan-like dislocation loop is left around the infinitely strong inclusion, providing a new route in understanding of the Bauschinger effect in realistic precipitation-strengthened material structures. Figure ( 13) summarizes the various stages of the dislocation loop propagation bypassing the inclusion in the polycrystalline