Computer Modeling in Engineering & Sciences

DOI: 10.32604/cmes.2021.016756

ARTICLE

Geometrically-Compatible Dislocation Pattern and Modeling of Crystal Plasticity in Body-Centered Cubic (BCC) Crystal at Micron Scale

Dedicated to Professor Karl Stark Pister for his 95th birthday

Yuxi Xie and Shaofan Li*

Department of Civil and Environmental Engineering, The University of California, Berkeley, CA 94720, USA
*Corresponding Author: Shaofan Li. Email: shaofan@berkeley.edu
Received: 24 March 2021; Accepted: 03 June 2021

1  Introduction

Crystal plasticity is a very common physical phenomenon. The origin of crystal plasticity is crystal dislocation motion and evolution. However, the crystal plasticity is resulted from aggregated dislocation motions, i.e., motion of dislocation pattern, rather than the motion of a single dislocation. It is because not only dislocation itself but also mutual interaction and evolution of assemble dislocations influence plastic behavior in crystalline materials. Therefore, the massive number of dislocations assembling in some specific regions in crystalline materials so that studying statistical variables from the aggregation is the only feasible way.

Based on various topological grouping of dislocations, some dislocation patterns have specific names. For example, we usually name the “wall-shape” dislocations as dislocation cells, the “strip-shape” dislocations as veins and mutually crossed “grid-shape” as dislocation labyrinth. According to past studies on crystal plasticity, many plastic phenomena in crystalline materials are strongly influenced by those micro-defects.

Recently, Li et al. [1–4] have proposed and developed a Multiscale Crystal Dislocation Dynamics (MCDD) method in which a physics-informed mesh is proposed and an atomistic-informed constitutive model and strain gradient-based atomistic kinematic rule are applied. Differen from phenomenological models, this method is an atom-level first principle method, i.e., there are no non-physical parameters to be fitted and inelastic behavior of crystalline materials are modelled by crystal lattice-dependent mesh and nonlinear kinematic rule so that this method can directly reflect essences of plasticity in crystalline materials. Besides, compared to Molecular Dynamics (MD) or Density-functional Theory (DFT), the crystal lattice-dependent mesh is a coarse-grained meshing method which can be scaled up to form mesoscale or even microscale models. Therefore, computational efficiency of the MCDD is much higher than MD and DFT. More importantly, a novel physical concept of geometrically compatible dislocation pattern (GCDP) is proposed to link the crystal lattice-dependent mesh and real topological structures of dislocation patterns in crystalline materials at least in early phase of crystal plasticity. This is because the initial dislocation formation or crystal defect always occurs on a lattice plane as forms of lattice complex (simplicial complex), it is reasonable to postulate that their aggregation will form a dual lattice complex, or CW complex [5,2]. In algebraic topology, a CW complex is an abbreviation for “Closure-finite Weak topology” complex, which is a kind of the Hausdorff space that can be partitioned into many open cells. In a CW complex, a simplex σk is homeomorphic to Bk,

Bk={x∈Rk|||x||2≤1}(1)

In this work, we use the concept of CW complex to construct geometrically compatible dislocation patterns, and such as a class of specific geometric objects, geometrically compatible dislocation patterns, have only been occasionally mentioned or reported in the literature, e.g., [5,6] but without systematical studies.

Tantalum, as lustrous transition metal, is a highly corrosion-resistant material, and it is widely used as minor components in many alloys, as part of the refractory metals group. In current industry applications, tantalum is mainly used in following four fields: capacitors, chemical industry, alloy and hard metals [7,8,9]. Due to its excellent corrosion resistance to oxide layer, tantalum can be applied to components manufactured for installations under corrosion attack, such as heat exchangers, heating spark plugs, coolers, pipes and coatings for vessels in chemical engineering. In superalloys, i.e., alloys based on nickel, cobalt and iron-nickel, introduction of tantalum improves the corrosion resistance, increases grain finesse, makes carbon stable and improves hardness and strength. Tantalum oxide has high dielectric coefficient and good stability. Therefore, tantalum can be manufactured as powder, wire, foil and coatings in capacitors. Besides, the hard metals industry consumes about 30% tantalum production. Tantalum is used in mixed carbides addition which can improve the erosion resistance, hot strength, hardness, ultimate bending strength, toughness, shock resistance and grain refinement.

In present work, we employed the multiscale crystal defect dynamics, or the geometrically-compatible dislocation pattern dynamics, to investigate the crystal plasticity of a special BCC crystal–α-Tantalum at micron scale (up to 3 μm × 3 μm × 6 μm) to provide further validation of MCDD method. Furthermore, we have used MCDD to investigated the size effect of crystal plasticity in the study of crystal plasticity of single crystal α-Ta.

The rest of the paper is organized as follows: based on algebraic topology, construction of crystal lattice-dependent mesh and resemblance between the physics-informed mesh and real dislocation patterns will be discussed in Section 2. Then the atom-level first principled constitutive model for single Tantalum is presented in Section 3. In Section 4, the strain gradient based nonlinear kinematic rule are demonstrated. And the numerical examples are provided in Section 5. We conclude the study with summaries and comments in Section 6.

2  MCDD Model for BCC Crystal

In this section, we first introduce the dual lattice tessellation to partition a BCC crystal into a finite element mesh of lattice process zones, which was proposed in [1]. Then we shall discuss how to use discrete exterior calculus and algebraic topology to offer mathematical description of process zone tessellation [2]. Finally by introducing the notion of geometrically-compatible dislocation pattern, we give the physical meaning of the lattice process zone, which are the generic dislocation patterns [3,4].

2.1 Dual Lattice Process Zone Model

In discrete dislocation dynamics (DD), the smallest unit of dislocation dynamics is the dislocation segment [10]. Different from dislocation dynamics, the smallest unit of dislocation pattern dynamics is discrete dislocation pattern. To start with, we need first to construct dislocation pattern elements in a crystal lattice space.

To identify the generic dislocation pattern for a given lattice structure, we need first to construct a Voronoi cell for an representative atom: placing planes normal to line segments formed by two neighboring atoms of lattice complex at midpoints. These truncation planes will form a convex polyhedron around the representative atom, and the convex polyhedron is called the Wigner-Seitz cell in crystallography for Bravais lattices. Using the terminology of crystal chemistry, we may call it as the Voronoi-Dirichelt Polyhedron (VDP). However, VDP is a bigger set of unit cells that includes the Wigner-Seitz cell. This is because different choice of representative atoms will lead to different shape of VDP. For our current BCC model, we choose the nearest 8 atoms and the second nearest 6 atoms as representative atoms. Therefore, the corresponding VDP is a truncated octahedron or the standard Wigner-Seitz cell for the BCC lattice as shown in Fig. 1a.

Figure 1: (a) A BCC Wigner-Seitz cell and (b) One third-order process zone and 14 s-order process zone

The complete steps to have dual lattice process zone tessellation:

   (1) Based on crystal lattices of whole atoms in 2D or 3D space, we can have correspondingdual lattice;

   (2) Then we shall scale the VDPs down, put them at atom positions and term it as the highestorder process zone element;

   (3) It is found that the whole two-dimensional or three-dimensional space cannot be coveredby the scaled VDP. Therefore, corresponding different shape of elements are used to fill inthe gaps among scaled VDP. Therefore, for BCC crystal lattice, Fig. 1b shows two typesof prism elements (hexagonal cross section and square cross section). Figs. 2a–2f show theconnections of prism element-to-truncated octahedron element and prism element-to-wedge element. Figs. 2g–2l show the connections of prism element-to-tetrahedron element.

Figure 2: (a) There are six slip planes and 36 first-order process zone elements ((a)--(f)); and 24 bulk tetrahedron elements ((g)--(l)) in a BCC super dual-lattice unit cell

Therefore, in a BCC dual-lattice unit, we have (1) one 3rd order process zone, i.e., scaled down truncated octahedron element; (2) fourteen 2nd order process zone elements including 8 hexagonal prism elements and 6 square prism elements; (3) thirty-six 1st process zone elements; (4) twenty-four tetrahedron elements.

2.2 Exterior Calculus Representation on BCC Lattice Complex

Inspired from Ariza et al. [11], we also follow the concepts and terminology in algebraic topology [12,13] and give quantitative description of super dual lattice complex in BCC crystals.

Different from Ariza and Ortiz's work, the present work is about dislocation pattern dynamics, and we cannot use simplicial complex p-cell to describe discrete dislocation pattern units. However, we still found that the mathematical theory of algebraic topology and exterior differential calculus can also applied to describe dislocation pattern dynamics in terms of p-cell of CW complexes. In this work, we introduce following cell index system to label points, edges, faces and volumes in different discrete dislocation pattern units, i.e., truncated octahedron elements, prism elements, wedge elements and tetrahedron elements. We use a unified symbol, τji(l(ℓ),m), to label to various geometric objects, which are explained as follows: Here i represents the dimension of the polytopal cell, subscript j denotes the type of three-dimensional (3D) polytopal element, argument l(ℓ) denotes the number sequence of the 3D polytopal element at the lattice site l, and m is the number sequence of the polytopal complex indicated within the polytopal element τj3(l(ℓ)). For example, we can define the local index for every vertices on τ33(l), i.e., τ30(l,m),m=1,2,…,24, and then we can define local index for every segment and facet cells on τ33(l), i.e., τ31(l,m),m=1,2,…,36, and τ32(l,m),m=1,…,14. We can also use τ20(l(ℓ),m),m=1,…,8 for square prism element to label indices on prism elements τ23(l(ℓ)),ℓ=1,2,…,14, or τ20(l(ℓ),m),m=1,…,14 for hexagonal prism element to label indices. If we follow this convention on labeling wedge elements and tetrahedron elements, we can define every local cell index for every 3D polytopal elements for a given lattice site l.

In the next we shall take τ21(l(5),1) as an example to illustrate the application of those labeling and those operators. As we know that in the label τji(l(ℓ),m), i represents the dimension of polytopal cell and j represents the which type of polytopal element. Therefore, we have superscript i = 1 represents edge element, subscript j represents prism element, l(5) represents No. 5 prism element in Fig. 6. 1 represents the No. 1 edge in No. 5 prism element as shown in Fig. 3.

Figure 3: Polytopal labeling

This labeling process is important in finite element method (FEM) formulation and implementation, because we need to use the exterior calculus symbol to construct the connectivity array of FEM code that maps the local indices into global indices.

By using boundary operator of exterior calculus, we can represent lines and faces. For example, as shown in Fig. 6, for central scaled truncated octahedron τ33(l), we have

∂τ31(l,1)=τ30(l,5)−τ30(l,1)∂τ32(l,1)=τ31(l,1)+τ31(l,5)+τ31(l,6)+τ31(l,4)−τ31(l,2)−τ31(l,3)(2)

and

∂τ33(l,1)=τ32(l,1)+τ32(l,2)+τ32(l,3)+τ32(l,4)+τ32(l,5)+τ32(l,6)+τ32(l,7)−τ32(l,8)−τ32(l,9)−τ32(l,10)−τ32(l,11)−τ32(l,12)−τ32(l,13)−τ32(l,14)(3)

which are shown in Fig. 4

Figure 4: Boundary operator illustration:from vertex to edge, from edge to surface, from surface to volume

Similarly, we can use co-boundary operator δ to help us understand the relations among different types of elements. For example, for the central scaled truncated octahedron τ33(l), we have

δτ30(l,1)=τ31(l,7)−τ31(l,1)−τ31(l,2)−τ21(l(7),1)−τ21(l(3),2)−τ21(l(14),3)δτ31(l,1)=τ32(l,1)−τ32(l,7)+τ22(l(7),1)−τ22(l(3),2)δτ32(l,1)=τ33(l)−τ23(l,7)(4)

which established the links from vertex to edge, edge to surface and surface to volume as shown in Fig. 5.

Figure 5: Co-boundary operator illustration: from vertex to edge, from edge to surface, from surface to volume

Figure 6: Lattice complex representation of the BCC lattice and indexing scheme: (a) 2nd-order process zone elements and (b) the prism element No. 5

Furthermore, we can also apply the cell complex notation defined above to represent other types of elements. For example, the prism element τ23(l(i)) (6). There are seven facets for NO.5 prism element, i.e., ∂τ21(l(5),j),j=1,2,…,7. We write the exterior representation by using boundary operator,

∂τ21(l(5),1)=−τ30(l−ε3,9)+τ30(l,4)∂τ21(l(5),2)=−τ30(l−ε3,10)+τ30(l,3)(5)

in which, τ30(l−εi,j) denote it is the j-th vertex (j = 1,2,…,24) of on the l + ɛi,i = 1,2,…,7 polytopal cell. Also, we have

∂τ22(l(5),1)=τ21(l(5),1)−τ21(l(5),2)+τ31(l,13)−τ31(l−ε3,5)(6)

∂τ23(l,5)=τ22(l(5),1)−τ22(l(5),2)+τ22(l(5),3)+τ22(l(5),4)−τ22(l(5),5)−τ22(l(5),6)−τ32(l,3)−τ32(l−ε3,6)(7)

2.3 Geometrically Compatible Dislocation Pattern

In Lyu et al.’s work [3], an innovative concept, Geometrically-compatible dislocation pattern (GCDP) that offers a physical foundation of pre-dislocation pattern mesh. As Lyu et al.’s work, in this section we give geometrically-compatible dislocation pattern formation of BCC crystal structure.

Since we know that body-centered cubic (BCC) crystal structure is NOT the closed-packed structure, the most likely slip plane should be {110}, which is the secondary most packed plane. Based on Figs. 7a and 7b, we project the DLPZ tilling of BCC structure into {110} plane which corresponds to the transformation from Figs. 7a to 7b and therefore we have 2D illustrative Geometrically compatible dislocation pattern as shown in Fig. 7c.

Figure 7: BCC (110) plane projection of DLPZ tilling. (a) 3D dual-lattice process zone unit cell of BCC crystal (b) The projection of 3D dual-lattice process zone unit cell onto (110) plane (c) 2D geometrically-compatible dislocation pattern

From the projection from 3D to 2D, we understand that 2D triangle bulk element is actually the projection of the first-order process zone element in 3D DLPZ tilling, rather than the projection of 3D tetrahedron bulk crystal element; the projection of the second-order process zone element (prism element) in {110} plane is the 2D rectangular element, and 3D third-order process zone element (truncated octahedron element) will be projected to be the 2D honeycomb element as shown in Fig. 7.

Since we have the two-dimensional geometrically-compatible dislocation pattern as in Fig. 7c, we can have the following scenario of dislocation pattern formation process in body-centered cubic materials. At the very beginning stage of the dislocation multiplication, statistical number of parallel slips have pre-existed in different slip systems along different atomic planes, and those parallel slips will converge to a bulk junction as shown in Fig. 8. From energetic perspective, maintaining geometric compatibility and local crystal structure integrity may be a necessity for deformed crystal, therefore these parallel slip planes in different slip systems must form a stable and geometrically-compatible dislocation networks. We call the networks as geometrically-compatible dislocation pattern (GCDP) (see Fig. 8c). Therefore, we identify the geometrically-compatible dislocation pattern as substructures that have the shape of the Wigner-Seitz cell of the original crystal lattice.

Figure 8: Illustration of a possible formation of a two-dimensional geometrically-compatible dislocation pattern: (a) Dislocations in slip planes, (b) Aggregated parallel slips in multiple slip systems converge to a bulk junction, and (c) A two-dimensional geometrically-compatible dislocation pattern

3  Atomistic-Informed Constitutive Relation

3.1 Embedded-Atom Potentials for BCC

The embedded-atom model (EAM) is widely used atomistic potential model for metals in molecular dynamics, which takes into account interactions between atom nucleus as well as electron density contribution. The EAM method has been extensively used in molecular dynamics for metallic materials,

Based on general formulation of EAM [14,15], the EAM potential could be written as:

Etot=12∑ijϕ(rij)+∑iF(ρ¯i),i,j=1,2,…nb(8)

in which, ϕ (rij) is the nucleus pair interaction energy, F is the electron embedding energy, ρ¯i is the electron density of atom i, which is defined as:

where ρj(rij) is the electron density of atom i from atom j, and rij = |rj − ri| is the distance between atom i and atom j. Since in present work we are modeling and simulating crystal plasticity in α-phase Tantalum, which is a BCC crystalline material, we shall discuss the detailed BCC EAM formulation in Section 5.

Based on EAM formulation and stress-work conjugate relation, we can have a fourth-order elastic tensor formulation:

C=∂2W∂E⊗∂E|E=0=1Ω0{∑j[12(ϕ′′−1rjϕ′)+F′(ρ¯)(∂2ρ∂r2−1rj∂ρ∂r)]×Rj⊗Rj⊗Rj⊗Rjrj2+F′′(ρ¯)(∑j∂ρ∂rRj⊗Rjrj)⊗(∑j∂ρ∂rRj⊗Rjrj)}.(9)

The macroscale constitutive relation is expressed in terms of several important stress measures as follows:

P=∂W∂F=1Ω0u∑j=1nb(F′(ρ¯)ρ′(rj)+12ϕ′(rj))rj⊗RjrjS=1Ω0u∑j=1nb(F′(ρ¯)ρ′(rj)+12ϕ′(rj))Rj⊗Rjrjσ=1Ωu∑j=1nb(F′(ρ¯)ρ′(rj)+12ϕ′(rj))rj⊗rjrj.(10)

in which P is the first Piola-Kirchhoff stress, S is the second Piola-Kirchhoff stress and σ is the Cauchy stress, Ωu0 is the volume of crystal unit cell in the referential configuration, ωu is the volume of crystal unit cell in the spatial configuration, R is the undeformed atom bond, r is the deformed atom bond, F is deformation gradient in the crystal unit cell and W is the strain energy density of the crystal unit cell, nb represents the number of neighboring atoms in the crystal unit cell.

To find the initial equilibrium state for many-body potential problem, we use the general or generic analytical form of EAM potential reported in literature. To find the initial equilibrium state, we first consider the initial equilibrium condition_, which is defined as the stress free state in an undeformed configuration,

S=12∑j=1nbϕ′(rj)Rj⊗Rjrj+F′(ρ¯)∑j=1nbρ′(rj)Rj⊗Rjrj=0.(11)

Based on Eq. (11), we can compute the equilibrium lattice constant a, which satisfies the initial equilibrium condition. Thus it becomes an optimization problem with one variable and we firstly find the derivative of norm of S about lattice constant a:

∂∥S∥F∂a=∂∥S∥F∂S:∂S∂a=2S:∂S∂a=0(12)

To solve this nonlinear equation, we use the Newton-Raphson method to linearize the above equation,

∂∥S∥F∂a(ai+1)=∂∥S∥F∂a(ai)+∂2∥S∥F∂a2Δai=0(13)

Therefore, we can obtain the i-th iteration result ai as:

Δai=−∂∥S∥F∂a(ai)/∂2∥S∥F∂a2(14)

the iterative solution stops when until δ ai < tol.

The analytical expression of ∂2∥S∥F∂a2 is given as:

∂2∥S∥F∂a2=∂∂a(2S:∂S∂a)=2∂S∂a:∂S∂a+2S:∂2S∂a2(15)

in which

+F′(ρ¯)ρ′(rj)+12ϕ′(rj)rj]Rj⊗Rj}(16)

+12rjϕ′′′(rj)+F′′(ρ¯)(ρ′(rj))2+F′(ρ¯)ρ′′(rj)+12ϕ′′(rj)}Rj⊗Rj.(17)

We use three hierarchical order strain gradients to model them (Eqs. (16), (17)),

G=∂2x∂X⊗∂X,H=∂3x∂X⊗∂X⊗∂X,K=∂4x∂X⊗∂X⊗∂X⊗∂X.(18)

where x represents deformed material point coordinates in spatial configuration, X represents the initial material point coordinates in referential configuration.

Then the bond vectors in corresponding dislocation pattern segments may be defined by the higher order Cauchy-Born model as follows:

rj=F⋅Rj+12G:(Rj⊗Rj),(19)

rj=F⋅Rj+12!G:(Rj⊗Rj)+13!H⋮(Rj⊗Rj⊗Rj),(20)

rj=F⋅Rj+12!G:(Rj⊗Rj)+13!H⋮(Rj⊗Rj⊗Rj)+14!K::(Rj⊗Rj⊗Rj⊗Rj),(21)

respectively. Subsequently, as a function of deformation gradients, the strain energy density Ww, Wp, and Wt of the corresponding dislocation pattern elements under consideration, which can be expressed as follows:

Ww=1Ω0vc[F(ρ¯)+12∑i=1nbϕ(|F⋅Ri+12G:(Ri⊗Ri)|)],(22)

Wp=1Ω0vc[F(ρ¯)+12∑i=1nbϕ(|F⋅Rj+12!G:(Rj⊗Rj)+13!H⋮(Rj⊗Rj⊗Rj)|)],(23)

+13!H⋮(Rj⊗Rj⊗Rj)+14!K::(Rj⊗Rj⊗Rj⊗Rj)|)].(24)

The constitutive relations in wedge, prism and truncated octahedron (tetradecahedron) elements can be further explicated as follows:

P=∂Ww(F,G)∂F, and Q=∂Ww(F,G)∂G;(25)

P=∂Wp(F,G,H)∂F,Q=∂Wp(F,G,H)∂G, and U=∂Wp(F,G,H)∂H;(26)

P=∂Wt(F,G,H,K)∂F,Q=∂Wt(F,G,H,K)∂G,U=∂Wt(F,G,H,K)∂H, andV=∂Wt(F,G,H,K)∂K.(27)

in which P is the first Piola-Kirchhoff stress, which is coupled with the deformation gradient F; Q is the second-order stress coupled with the 1st-order strain gradient G; U is the third-order stress coupled with the 2nd-order strain gradient H, and V is the fourth-order stress coupled with the 3rd-order strain gradient K.

The higher order stress tensor are expressed in a general form as:

Q=12Ω0∑j=1nb(12ϕ′(rj)+F′(ρ¯)ρ′(rj))rj⊗Rj⊗Rjrj,(28)

U=12Ω0∑j=1nb(12ϕ′(rj)+F′(ρ¯)ρ′(rj))rj⊗Rj⊗Rj⊗Rjrj.(29)

V=12Ω0∑j=1nb(12ϕ′(rj)+F′(ρ¯)ρ′(rj))rj⊗Rj⊗Rj⊗Rj⊗Rjrj.(30)

4  MCDD Finite Element Formulation

To formulate MCDD finite element formulation, we first consider the Hamilton principle in terms of displacement variation for a fixed time interval as:

δ∫t0t1(T−Wint)dt+∫t0t1δWextdt=0(31)

in which, Wint=∫VW(F,G,H,K)dV is the strain and higher order stain gradient energy density, T=∫V12ρu˙⋅u˙dV is the kinetic energy of the crystalline solid, Wext is the external potential energy.

The variation of the kinetic energy of the crystal is:

The internal virtual work of the crystal is:

and the external virtual work of the crystal is:

Substituting these formulas above into the Hamilton principle, we can have Galerkin weak form of Multiscale Crystal Defects Dynamics for a crystal as shown as follows:

∫Ω[ρu¨⋅δu+P:δF+Q⋮δG+U::δH+V:⋮δK]dV=∫Ωb⋅δudV+∫∂ΩtT¯⋅δudS(32)

Considering finite element discretization and interpolation,

where Bi(X) is the finite element shape function, and ui are the element nodal displacements, we can rewrite the form of considering different types of elements as:

∑e=1nelemp0{∫Ωp0eρ0u¨h⋅δuhdV+∫Ωp0eP:δFhdV}+∑e=1nelemp1{∫Ωp1eρ0u¨h⋅δuhdV+∫Ωp1e(P:δFhdV+Q⋮δGh)dV}+∑e=1nelemp2{∫Ωp2eρ0u¨h⋅δuhdV+∫Ωp2e(P:δFhdV+Q⋮δGh+U::δHh)dV}+∑e=1nelemp3{∫Ωp3eρ0u¨h⋅δuhdV+∫Ωp3e(P:δFhdV+Q⋮δGh+U::δHh+V:⋮δKh)dV}=∑i=03∑ei=1nelempi{∫Ωpieiρ0b⋅δuhdV}+∑e0=1nelemp0∫∂tΩp0e0T¯⋅δuhdS(33)

where b is the body force, Ωeipi is the ei element in the p-th process zone, ∂tΩeipi is the interception between the traction boundary and the boundary of the ei-th element in the i-th process zone ωpi, ρ0 is crystalline density in referential configuration. We use the notation nelempj,j=0,1,2,3 to represent the total number of elements for each type of process zone element.

Note that based on internal virtual work,

δWint=∫V{∂W∂F:δF+∂W∂G:δG+∂W∂H:δH+∂W∂K:δK}dV(34)

the successive integration by parts yields the following expression:

δW=−∫V∇X⋅{P−∇X⋅[Q−∇X⋅(U−∇X⋅V)]}⋅δudV+∫∂V{N⋅[P−∇X⋅(Q−∇X⋅(U−∇X⋅V))]}⋅δudS+∫∂V{N⋅(Q−∇X⋅(U−∇X⋅V))}:δFdS+∫∂V(N⋅(U−∇X⋅V))⋮δGdS+∫∂V(N⋅V)::δHdS(35)

In our current MCDD FLEM model, we assume that on the boundary of the crystalline solid the higher order stress effects are negligible, i.e.,

Q~≈0, where Q~=∇X⋅(U−∇X⋅V),∀X∈∂VtU~≈0, where U~=∇X⋅V,∀X∈∂VtV~≈0, where V~=V,∀X∈∂Vt(36)

This is a very coarse approximation, and further work is needed to consider higher order strain gradient traction boundary condition.

5  Numerical Simulations

5.1 Uniaxial Tension and Pure Shear Tests

In this section, we first present MCDD numerical results in uniaxial tension and pure shear tests. To model and simulate crystal plasticity in α-phase tantalum crystal, we adopt the higher order Cauchy-Born rule that utilizes the Embedded-Atom method potential discussed in [15].

The total energy E of the crystal solid can be written as:

E=12∑i≠jϕ(rij)+∑iFi(ρ¯i),(37)

where ϕ (rij) represents the pair energy between atoms i and j separated by a distance rij, and Fi is the embedding energy associated with embedding an atom i into a local site with an electron density ρ¯i, which is defined as:

ρ¯i=∑i≠jρj(rij)(38)

with ρj(rij) as the electron density at the site of atom i arising from atom j at a distance rij away.

We use the following pair potential ϕ(r) for α-Ta,

ϕ(r)=Aexp[−α(rre−1)]1+(rre−κ)20−Bexp[−β(rre−1)]1+(rre−λ)20,(39)

where re is the equilibrium spacing between nearest neighbors, A,B,α and β are four adjustable parameters, and κ and λ are two additional parameters for the cutoff. The electron density ρ(r) is defined as:

ρ(r)=feexp[−β(rre−1)]1+(rre−λ)20.(40)

The embedding energy function is expressed as follows:

F(ρ)=∑i=03Fni(ρρn−1)i,ρ<ρn,ρn=0.85ρeF(ρ)=∑i=03Fi(ρρe−1)i,ρn≤ρ<ρ0,ρ0=1.15ρeF(ρ)=Fe[1−ln(ρρsη)](ρρs)η,ρ0≤ρ(41)

The value of parameters used for α-phase Tantalum are taken from in [18]. Based on the EAM formulation, i.e., Eq. (9), we first calculate elastic modulus of the crystal and compare it with experimental results, and the results are listed in Appendix Table A1.

The MCDD model with the size of 100A∘×100A∘×100A∘ is shown as in Fig. 9. In this example, there are 576 tetrahedron elements, 1008 wedge elements, 154 prism elements (including 88 hexagonal prism elements and 66 square prism elements) and 15 truncated octahedron elements. We firstly apply the prescribed displacement/velocity loading condition on MCDD specimen under prescribed loading and unloading strain rate with 0.01 ps–1. The boundary conditions for MCDD simulations are as follows: The upper 10% portion of the specimen is applied with the prescribed displacement/velocity, and the unloading is initiated at the three different strain marks set at 10%, 15% and 20% with the same strain rate. The lower 10% portion of the specimen is fixed. The rest of the boundaries are traction and higher order traction free (see [4]).

Figure 9: MCDD modeling and simulation of a BCC crystal under uniaxial tension/compression loadings

Meanwhile, molecular dynamics (MD) simulation is conducted by using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS), which is a popular open-source molecular dynamics simulation software package. The MD model or specimen has 59582 atoms. The atomic pair potential for α-Ta is based the reference [15]. The same strain rate used in MCDD calculations is applied in MD simulations.

By choosing the unloading position at three different strain marks mentioned above, we have obtained the corresponding stress-strain relations in MCDD simulations. The obtained stress-strain relations are compared with those of MD simulations as shown in Fig. 10, respectively. From Fig. 10, one may find that in all three loading/unloading cases the stress-strain relations obtained from MCDD simulations and from MD simulations are consistent, which indicates the MCDD method is valid in simulating nanoscale crystal plasticity for single α-Ta crystal.

Figure 10: Stress-strain relations for three unloading strains: (a) 20%. (b) 15%. (c) 10% under loading and unloading

By considering higher order Cauchy-Born rule based higher order strain gradient formulations on wedge shaped dislocation patterns, prism shaped dislocation patterns, and truncated octahedron dislocation patterns, the MCDD method shows its ability to capture inelastic behavior and the path-dependent constitutive relation of α-phase tantalum in single crystal plasticity simulation. Fig. 11 shows the von-Mise stress contours on (100) and (010) planes. It is remarkable that the above crystal plasticity simulation is based on local “higher order hyperelastic strain gradient formulation”. Moreover, the above MCDD simulation has much computational efficiency in comparison with the MD model or specimen with the same size and boundary and loading conditions.

Figure 11: The von-Mise stress distributions in (100) and (001) planes at the maximum loading points: (a), (d) for 20% strain unloading; (b), (e); for 15% strain unloading, and (c), (f) for 10% strain unloading

Similar to the pure tension case, we also conducted pure shear test on an α-Ta specimen of 100A∘×100A∘×100A∘. The setup of the MCDD simulation and α-phase tantalum model are shown as Fig. 12. The displacement/velocity boundary conditions under prescribed pure shear deformation are similar to that of the uniaxial loading one, except that the tensile loading becomes shear loading. For the MD simulation, the boundary condition on x,y and z directions are respectively p,s and s in which p represents periodic and s represents non-periodic, shrink-wrapped.

Figure 12: MCDD model under pure shear loading

By prescribing three different unloading strains, we can observe overall stress-strain relation in shear in comparable with that of MD simulations with the same loading and unloading conditions as shown in Fig. 13.

Figure 13: Shear stress-strain relations for three different unloading strain conditions: (a) at 20%. (b) at 15%. (c) at 10% shear strain unloading

Based on the results presented above, we may conclude that under pure shear loading MCDD method also offers consistent results with that obtained from MD simulations. It may be noted that in Fig. 14 the stress corresponding to the higher shear strain (maximum strain point) leads to higher residual stress.

5.2 The Relationship between Size Effect and Cross-Slip

Conventionally, size effect is not considered in constitutive relation therefore the corresponding stress-strain curve is smooth and continuous. However, with sample size decreasing to micron or sub-micron scale, material strength shows increased compared to macro-scale sample size. Some scholars give some statistical relations to describe size effect on material strength reduce. But until now, people are not clear what mechanism leads to size-effect of micron-scale materials.

Cross-slip is atomistic level phenomenon observed commonly in crystal plasticity especially when dislocation density attains to a certain value. Capturing cross-slip in micro-scale is not trivial, not to mention that cross-slip is observed in sub-micron scale.

Figure 14: The von-Mise stress distribution respectively at the maximum strain point (a), (b), (c) and the residual strain point (d), (e), (f): (a), (d) for 20% strain unloading; (b), (e); for 15% strain unloading, and (c), (f) for 10% strain unloading

MCDD, as a multiscale method, has advantages in simulating crystal in micron-scale. Therefore, in this numerical example, we can find that cross-slip can be captured in numerical sample with different size. Furthermore, we give a preliminary conclusion that the appearance of cross-slip may lead to size effect in crystal plasticity.

Based on [16], we firstly constructed the same size MCDD model as comparison to validate our MCDD method. The MCDD model is shown as Fig. 15. The initial aspect ratio is 1:2:4, with dimension 420A∘×840A∘×1680A∘. The upper 10% portion of the model is velocity-prescribed corresponding to three different strain rate: 1.1 × 107 s–1, 5.55 s × 107 s–1 and 2.77 s× 108 s–1. The MCDD results are compared with MD simulation as Fig. 15.

From the figure above, we can find that the yielding stress of MCDD model is much lower than perfect crystal microstructure of MD simulation results and it is approximate to results of MD model with voids inserted. Based on the model above, we construct different size models to study relation between size effect and cross-slip.

From the Fig. 16, we can find that size effect in yielding stress: larger size model has lower yielding stress. Furthermore, we are able to capture cross-slip in sub-micron scale crystal. Based on Fig. 16, we can conclude that the appearance sequence of cross slip may lead to the appearance sequence of yielding stress, i.e., size effect.

Figure 15: MCDD model and stress-strain curve comparison with MD simulation

Figure 16: Stress-strain curve comparison and cross-slip observation

Based on such simulation result, we try to have the connection between size effect and cross-slip whose mechanism is following: initially when number of dislocation is few, crystal only has the primary slip plane. With strain increases, the number of dislocation in unit area (i.e., dislocation density) will increase. When dislocation density attains at a specific value, the second slip plane will be activated therefore the prerequisite of dislocation jumping (i.e., cross-slip) forms. Under this condition, dislocation will “jump” from one plane to another. Therefore, the cross-slip will be formed. Until now, MD simulation is hard to capture cross-slip but MCDD method can capure it. The reason may be that single or very few dislocations cannot activate second slip plane which is the prerequisite of cross-slip. But the MCDD method, as a multiscale method, simulates dislocation pattern which represents the ensemble or sub-network of dislocations. Therefore, MCDD method can capture that formation of secondary slip plane and “jumping” of dislocation.

Figure 17: (a) MCDD FEM model and (b) Comparison of the stress-strain relations obtained from specimens of four different sizes

Based on analysis above, we know that dislocation density will activate cross-slip, the appearance sequence of cross-slip may determines the appearance sequence of yielding stress of materials and the difference among appearance sequence of yielding stress is the size effect. By this conclusion, we establish the connection between dislocation density and size effect and offer a preliminary explanation to size effect in micron-scale crystal sample.

5.3 Size Effect in Micron Scale Crystal Plasticity

At nanometer or sub-micron scale, a very strong size effect of crystal plasticity has been reported in experiments, e.g., [17]. In the following, we report a preliminary study on size-effect of meso-scale plasticity observed in MCDD simulations of a compression test of α-Ta.

We constructed three different MCDD numerical specimens at micron scale, and they are square columns with fixed length-to-cross section width ratio of 2:1. The four specimens have the sizes: Lx× Ly × Lz = 1.0 μm × 1.0 μm × 2.0 μm, 2.0 μm × 2.0 μm × 4.0 μm, and 3.0 μm × 3.0 μm × 6.0 μm as shown in Fig. 17a.

In this micron-scale MCDD specimen, it has 5952 tetrahedron elements (bulk crystal elements), 4096 wedge elements (1st order process zone or wedge dislocation pattern elements), 2752 prism elements (2nd order process zone or prim dislocation pattern elements), and 367 truncated octahedron elements (3rd process zone or truncated octahedron dislocation pattern elements). We then conducted uniaxial compression tests with prescribed displacement boundary condition with the loading strain rate: 5× 10 − 6 ps−1.

The stress-strain curves of α-Ta micropillar are shown in Fig. 17b. The yielding stress exhibits an increase with decreasing the size of MCDD specimen. The similar tendency have been observed in Abad et al. [17] in experiments. The increase in yield strength of micropillar might have been attributed to the variance of surface energy with micropillar size [18] as well as the less chance for defects.

6  Conclusions and Discussions

In this work, we have further developed the recently proposed multiscale dislocation pattern dynamics by carrying out meso-scale discrete dislocation pattern dynamics computations and simulating crystal plasticity in α-Ta at the length scale up to 3.0 μm × 3.0 μm × 6.0 μm, which is beyond the capacity of conventional molecular dynamics simulations. Moreover, we have demonstrated that as a dislocation pattern dynamics, MCDD method can capture the size effect of crystal plasticity, and accurately capture cross-slip in BCC crystals. Based on the numerical simulation results discussed above, we have found that MCDD simulation has almost the same accuracy with that of molecular dynamics, but higher computation efficiency. This allows us to study mesoscale crystal plasticity at micron scale, i.e., 3.0 μm3 above, which is not reachable for MD simulations under current computer technology.

This work is highlighted by following two advances: (1) The simulated model is able to show accurate inelastic stress-strain curve comparable to MD results, and we can even capture dislocation nucleation and growth process with obvious crystal plasticity. There is a significant development in computational efficiency, which allows us to imagine further step: simulation of super micron-scale model and comparing with experiment results; (2) MCDD method shows the potential to be used to simulate the macroscale inelastic deformation process and fatigue damage, (3) MCDD successfully capture cross-slip which is hard to be observed in MD simulations. And having a preliminary relation between cross-slip and size effect lead us to have a preliminary explanation to mechanism of size effect and (4) MCDD model demonstrates that it can successfully capture size effect of up to 3.0 μm3 micropillar crystal model, which shows validity of MCDD method in simulating crystal plasticity.

However, the method still has some unresolved issues. For example, the thermal effect that is strongly influenced dislocation nucleation and growth has not been taken into account in the present dislocation pattern dynamics. Moreover for macro-scale dislocation pattern dynamics simulations, it will need massive MPI parallel computation to realize. For the above issues and more, we hope to deal with them in the subsequent work.

Funding Statement: The authors received no specific funding for this research.

Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.

Appendix A. Elastic constant of Tantalum

In Appendix Table A1, we listed elastic constants of α-Ta obtained from both experiments [19] as well as the predictions from the EAM calculation.

This work is licensed under a Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.