Development of the Concurrent Multiscale Discrete-Continuum Model and Its Application in Plasticity Size Effect

Abstract

A concurrent multiscale model coupling discrete dislocation dynamics to the finite element method is developed to investigate the plastic mechanism of materials at micron/submicron length scales. In this model, the plastic strain is computed in discrete dislocation dynamics (DDD) and transferred to the finite element method (FEM) to participate in the constitutive law calculation, while the FEM solves the complex boundary problem for DDD simulation. The implementation of the whole coupling scheme takes advantage of user subroutines in the software ABAQUS. The data structures used for information transferring are introduced in detail. Moreover, a FE mesh-based regularization method is proposed to localize the discrete plastic strain to continuum material points. Uniaxial compression tests of single crystal micropillars are performed to validate the developed model. The results indicate the apparent dependence of yield stress on sample size, and its underlying mechanisms are also analyzed.

1. Introduction

The effect of sample size or shape on material deformation behavior [1,2,3] has received attention for decades. For instance, the sample size of concrete influences its kinetics of external sulfate attack [4]. At the microscale, it has been reported that the deformation mechanism of crystal materials is different from that of bulk materials [5,6,7]. It is recognized that the strength of micro/submicron crystals depends not only on the material itself but also the sample/grain size. Many efforts have been made to study the size effects of crystal materials at the micro/submicron scale in order to improve material performance for high-quality products [8,9,10,11]. The size effect is induced by many physical mechanisms, such as dislocation motion, nucleation, and starvation [12,13]. It remains difficult to capture the full picture of the deformation behavior of crystal material at the micro/nano-scale only by experimental tests. Therefore, it is necessary to employ a suitable simulation method to reveal the underlying deformation mechanisms of crystal materials at the micro/nano-scale.

Computational simulation methodologies have been continuously developed to predict material deformation behaviors. According to the difference in simulation scale, the main simulation methods can be broadly classified into the finite element method (FEM) [14,15], crystal plasticity finite element methods (CPFEMs) [16,17], discrete dislocation dynamics (DDD) [8,18,19], and molecular dynamics (MD) [20]. Based on empirical constitutive equations, the FEM shows a great advantage in dealing with various complex boundary problems at the macroscale. Compared to the classic FEM, CPFEMs consider the material microscopic plastic shear deformation by employing critical parameters that represent the microstructural state of materials beneath the materials’ macroscopic deformation. Typical examples are the slip systems-based phenomenological CPFEMs and the dislocation density-based CPFEM. Moreover, by considering the density of geometrically necessary dislocations (GNDs) in the constitutive equations, it can be used to predict the size effect of crystal materials because GNDs are related to the strain gradient and has been considered in the constitutive equation [21,22,23,24]. However, CPFEMs cannot provide an explicit microstructure evolution during the deformation process, because of the empirical constitutive equation used. To obtain insight into the physical deformation mechanism, a microscale simulation model that considers the intrinsic microstructure characteristics of crystal materials is needed. MD simulation directly depicts the interaction of atoms. In recent decades, it has been widely employed to investigate the mechanical properties of single crystals and nanopolycrystalline materials (e.g., nano-polycrystal silicon [25] and copper [26]) by means of nanoindentation or many other tests [11]. However, MD simulation is only suitable for nanoscale specimens due to the extremely heavy computational load. Indeed, the plastic deformation of crystals results from its dislocation activities. To overcome those limitations, DDD has been proposed by many researchers to capture the size effect and uncover its underlying mechanisms. It omits the direct atomic interactions and takes dislocation as the intrinsic carrier of plastic strain. DDD modeling is often employed together with the FEM to deal with complex boundary conditions. The mostly used coupling schemes are the discrete continuous model (DCM) [27] and superposition method (SPM) [28].

To solve the boundary value problem of a finite crystal, van der Giessen and Needleman [28] first proposed the concept of the SPM. In the SPM, the total stress field in a finite crystal is the sum of two stress fields. One is the analytical stress field due to the existence of dislocations in an infinite unbounded crystal and the other is the complementary elastic stress field calculated by finite element analysis in a finite crystal. The method is able to capture the short-range dislocation stress field accurately and has been improved by incorporating a 3D dislocation mechanism [13,29] and dislocation climb [30]. Nevertheless, the calculation of the analytical stress field between all existing dislocations is time-consuming. In addition, there is no explicit introduction of the plastic strain that arises from dislocation activities. The DCM is another hybrid approach coupling DDD to the FEM. The numerical calculation in the DCM is straightforward and has clear physical sense. The moving dislocation generates plastic strain along its gliding path, and then the plastic strain is localized to the Gaussian integration points of the FEM. By making use of the plastic strain, the equilibrium stress is computed in the FEM module and is used to drive the movement of discrete dislocations.

Since Lemarchand [27] first proposed the DCM in 2001, it has been widely used to study the plastic behavior of materials at micron/submicron length scales, for example, the effect of hard coating [31], low-amplitude cyclic loading [32], temperature [33], and shock loading [34] on the plastic deformation of micropillars. The dislocation mechanisms causing the channel size effects of superalloys [12], incipient plastic instability of nanoindentation [35], and strain bursts of micropillars [10] have been investigated using the DCM. Recently, Cui et al. [31] studied the confined plasticity in coated micropillars using the DCM and found that the operation of a single arm source plays a key role in the plastic deformation of coated pillars. Santos-Güemes et al. [36] studied the interaction mechanism of dislocation and precipitates using a developed multi-scale DCM. The intrinsic hardening mechanism of a nickel-based superalloy and its dependence on the interfacial dislocation network and lattice mismatch were reported [37]. By introducing a dislocation climb model into the DCM, Liu et al. [38] investigated the high-temperature anneal hardening behavior of an Au submicron-pillar.

Most studies are focused on the application of the DCM in plastic mechanisms at micron/submicron scales; however, less attention has been paid to its numerical algorithms for improved computational efficiency. Zbib et al. [39] proposed a simple method of regularizing the plastic strain induced by dislocations to the material points of their localized elements. Liu et al. [40] introduced the Burgers vector density function to apportion the discrete plastic strain to the corresponding material points. Vattré [41] and Cui [42] also discussed the detailed numerical algorithm of localizing the discrete strain to continuum material points in a 3D DCM simulation. In addition, the strain regularization method in a 2D DCM was developed by Huang [30] by considering the interaction of dislocation with its neighboring inter-phase boundaries and other dislocations. However, the computational load remains a challenge in computational simulation. There is always a balance of numerical accuracy and computational efficiency.

In this paper, a concurrent multiscale model coupling DDD to the FEM is developed to investigate the plastic mechanism of materials at micron/submicron length scales. A FE mesh-based regularization method is proposed to improve the computational efficiency. The implementation details of how to apply the DCM principle to couple DDD with FEM calculation modules during the deformation process are introduced. The optimized data structure and the achievement of the simultaneous coupling scheme by subroutines are presented for the first time. The description of the 2D-DCM is provided in Section 2. Section 3 discusses transfer methods of plastic strain and stress among DDD and FEM modules. The user subroutines of Abaqus used in the multiscale framework are described in detail in Section 4. In Section 5, a compression test of a single crystal is performed using the developed method, and it is followed with some main conclusions drawn in Section 6.

2. Description of DCM Coupling Framework

As schematically illustrated in Figure 1, the multiscale framework consists of two modules, i.e., DDD and FEM module. The key point of the 2D-DCM is that the plastic strain calculated by DDD directly participates in the constitutive law of the FEM. The stress tensor at simulation step n, ${\sigma }_n$, is expressed as,

$${\sigma }_n=[C^e]:\left({\epsilon }_n^t-{\epsilon }_{n-1}^p\right)$$

(1)

where ${\epsilon }_{n-1}^p$ is the plastic strain at the last step by DDD, ${\epsilon }_n^t$ is the total strain at simulation step $n$, and $C^e$ is the elastic moduli tensor. The explicit computation of movement, evolution, and interaction between dislocations and other defects (e.g., obstacles, dislocation, and sources) is performed in DDD. In addition, the calculation of the equilibrium stress field in the FEM module and the coupling algorithm of both modules are discussed.

2.1. Two-Dimensional Dislocation Dynamics

For the 2D plane strain problem discussed in this paper, only edge dislocations are considered. The effect of dislocation climb is also negligible compared with that of dislocation glide. The driving force applied on dislocations can be expressed as

$$f_i=n_i·{\sigma }_i·b_i,$$

(2)

where $n_i$ is a unit normal to the slip plane, $b_i$ is the Burgers vector, and ${\sigma }_i$ is the stress tensor at the position of dislocation i, which is evaluated in the FEM.

Assuming that the edge dislocations can reach a stable velocity instantly, the glide velocity of dislocation $i$ is

$$v^i=\frac{f^i}{B},$$

(3)

where B is the viscous drag coefficient.

Static point sources are distributed on the slip plane randomly to mimic the operation of Frank–Read sources. A dislocation dipole is generated once the shear stress at the point source exceeds its critical strength, ${\tau }_{nuc}$, during a time period. The initial distance $L_{nuc}$ between two newly generated dislocations is

$$L_{nuc}=\frac{\mu }{2\pi \left(1-\nu \right)} \frac{b}{{\tau }_{nuc}},$$

(4)

where $\mu$ is the shear modulus and $\nu$ is the Poisson ratio of crystal material.

Two dislocations with opposite Burgers vectors annihilate each other when their distance is less than the critical annihilation distance $L_e$, which is taken as 6b [43,44].

The motion of dislocations can be impeded by various obstacles, such as small precipitates and forest dislocation. The point obstacles are modeled in 2D-DDD to consider the effect of the obstacles existing in real crystals [43,45]. If one obstacle lies on the gliding path of a dislocation, the dislocation is pinned by the obstacle when its resolved shear stress is less than the strength of the obstacle. Dislocations can be de-pinned or bypass the obstacle if its shear stress exceeds the obstacle strength.

2.2. Calculation in Finite Element Module

The standard FEM is employed to solve the boundary value problem of a finite crystal as follows:

$$\left[K\right]\left\{U\right\}=\left\{f^a\right\}+\left\{f^p\right\},$$

(5)

$$\left\{f^p\right\}={{\displaystyle \int }}_V^{ }\left[C^e\right]{\epsilon }^p\left[B\right]dV,$$

(6)

where [K] is stiffness matrix, $\left\{U\right\}$ is the nodal displacement vector, $\left\{f^a\right\}$ is the applied force vector, $\left\{f^p\right\}$ is the force vector from plastic strain, [B] = grad[N], and [N] is the shape function vector. The DDD module contributes the plastic strain ${\epsilon }^p$ to the unified continuum mechanics framework.

3. Plastic Strain–Stress Transfer in Coupling Scheme

3.1. Plastic Strain and Stress Transfer Scheme

In this computational simulation model, the plastic strain is explicitly calculated in the DDD module and is then transferred to Gaussian integration points of the FEM module for the constitutive calculation. The plastic strain is induced by the glide of dislocation, which can be calculated as:

$${\epsilon }_i^p=\frac{A_i}{2V}\left(n_i\otimes b_i+b_i\otimes n_i\right),$$

(7)

where $A_i$ is the swept area of the dislocation i and V is the representative volume.

The plastic strain ${\epsilon }_i^p$ obtained from the DDD module should be located at the corresponding integration points as the eigen-strains. Different regularization methods have been proposed to solve this problem. In the early work by Lemarchand [27], it was assumed that a slab with finite thickness h surrounding the glide plane represents the elementary slip events. The elementary volume associated with each Gauss integration point is taken as the representative volume. The plastic strain distributed at each integration point depends on the intersection volume between the slab and elementary volume. Based on the work of Lemarchand [27], Liu [40] proposed a weight function (i.e., the Burgers vector density function) to distribute the total plastic strain to each integration point of finite element. Vattré [41] and Cui [42] also made efforts to improve the regularization method of plastic strain. For the 2D computational model, Huang [30] presented a detailed regularization procedure, which is similar to that of Liu [40]. A local domain along the slip plane of dislocation was defined and half of its width was set as $a$, which is analogous to the dislocation core spreading radius. The integration points in the local domain would be allocated with plastic strain, and the allocation proportion of each integration point was decided by its distance to the slip plane of dislocation. This regularization method can relieve a high gradient of plastic strain among neighboring elements. However, it causes a heavy computational burden, especially when numerous dislocations exist in the simulation model, due to the distances between each integration point, and the slip plane should be calculated.

An efficient regularization method is proposed in this paper by considering the characteristics of meshes of the FEM model. Elements in the FEM model are divided into elementary domains (i.e., integration point domain), which are associated with each integration point in the elements. As shown in Figure 2, a dislocation glides along its slip plane from D to D’. These elementary domains intersecting with the gliding path are regarded as ‘core domains,’ which are shown as the red domains in Figure 2. If the distance from the integration point of an elementary domain to a core domain is less than the given value $a$, the elementary domain k is defined as a ‘secondary domain.’ The hatched and meshed areas represent the secondary domains in Figure 2. The parameter $a$ has the same meaning and evaluation as the one used in Huang’s method. The core and secondary domains together compose the domain ${\mathsf{\Omega }}_i$. The localized plastic strain of integration point j, ${\epsilon }_{ij}^p$, in the domain ${\mathsf{\Omega }}_i$ is given according to a weight function,

$${\epsilon }_{ij}^p=\frac{{\omega }_j}{{{\displaystyle \sum }}_j{\omega }_j}\frac{A_i}{2V_{ }}\left(n_i\otimes b_i+b_i\otimes n_i\right),$$

(8)

where ${\omega }_j$ is the distribution coefficient for plastic strain caused by dislocation i on integration point $j$. The value of ${\omega }_j$, which is selected as the Burgers vector spreading function developed by Cai [46], can be calculated as

$${\omega }_j=\frac{1}{{\left(r_j^2+a^2\right)}^{3.5}},$$

(9)

where $r_j$ is the distance between the integration points of the j-th elementary domain and its nearest core domain. For the integration point $O$ shown in Figure 2, the distance is equal to the length of $OB$, and it replaces the length of $OA$ used by Huang [47]. The simplification of $r_j$ makes it easier to obtain the distribution coefficient. Instead of solving the formula for the distance from point to line, it only needs the distance between the integration points of the core and secondary domain. Once the core domains are obtained, the secondary domains around them can be quickly captured according to the numbering order of integration points and elements. Figure 3 shows the numbering order of integration points in different elements. The new method omits the determination of the gliding path equations of each dislocation by replacing the glide path with the core domain. Thus, with the present method, it is computationally efficient to capture the influenced domain, ${\mathsf{\Omega }}_i$, of each moving dislocation and calculate the distribution coefficient through directly solving the distance formula of two points.

After FEM calculation, there are two steps to transfer the stress field from the FEM module to the DDD module, i.e., (i) stress extrapolation from integration points to nodes by the inverse shape function, and (ii) stress interpolation from nodes to dislocations by the shape function of the element [41]. In finite element analysis, the numerical integration is carried out at the Gaussian points and the stress value at the integration point is the most accurate. Therefore, in this work, the stress of dislocations is determined by the integration point domain where dislocations reside. The stress tensor of a dislocation is the same as its nearest Gaussian integration point. The division of various elements into integration point domains is illustrated in Figure 3. Each domain contains one integration point and any point in the domain has the shortest distance to the contained integration point compared to other integration points in the parent element. The number of integration points in the parent element is determined by the degree of the integrated polynomial terms in the element’s stiffness matrix and the desired accuracy of results. According to the number of integration points required, their locations in the parent element can be obtained by

Table 1. Gauss numerical integration constants.

Number of Integration Points (n)	$\mathbf{Integration} \mathbf{Point} \mathbf{Location} {\mathit{x}}_{\mathit{i}}$
1	$x_1=0$
2	$x_1=1/\sqrt{3}$$, x_2=-1/\sqrt{3}$
3	$x_1=\sqrt{0.6}$$, x_2=0$$, x_3=-\sqrt{0.6}$

.

3.2. Coordinate System Conversion Involved in Coupling Scheme

In order to complete the data exchange between DDD and the FEM module, the transformation of the coordinate system should be taken into account.

The calculation in the DDD module, such as glide force, dislocation velocity, and plastic strain, is finished on the local coordinate system whose axes are consistent with the slip direction and normal direction of the slip plane, as shown in Figure 4a. The FEM module is constructed in the sample coordinate system, unlike that of the DDD module. To achieve the coordinate system transformation, the transformation matrix, $\left[R\right]$, is calculated by representing the vector $OA$ in the two coordinate frames, as illustrated in Figure 4b. Then, the vectors, {g}, and tensors, [M], can be transformed with rules as follows:

$$\left\{g^{\prime }\right\}=\left[R\right]\left\{g\right\},$$

(10)

$$\left[M^{\prime }\right]=\left[R^T\right]\left[M\right]\left[R\right],$$

(11)

where $\left\{g^{\prime }\right\}$ and $\left[M^{\prime }\right]$ are the transformed vector and tensor, respectively. Superscript T means matrix transpose.

4. Subroutines for ABAQUS

The multiscale framework is built by the secondary development of finite element software ABAQUS. In this section, the implementation details of the DCM model in ABAQUS are introduced.

4.1. User Subroutines in the Multiscale Framework

Abundant user subroutines are provided in ABAQUS, which allow advanced users to customize a variety of ABAQUS capabilities [48]. The general workflow for running the DCM using ABAQUS is illustrated in Figure 5. Three user subroutines are used for the DCM model, including URDFIL for results processing, UEXTERNALDB for external file manipulation, and UMAT for the physics-based constitutive law. The functions of these subroutines in the developed model are discussed below.

A user subroutine called URDFIL is written particularly as a stress relay from the FEM to DDD module. The subroutine performs two functions: it can access the results file for stress tensors of Gaussian integration points, nodes coordinates, elements, and its relative integration point numbers at the end of the increment, and it transfers this new information to the DDD module by writing them to user-defined COMMON blocks.

The DDD-relative calculations (e.g., the velocity of dislocations, dislocation multiplication, and annihilation) are implemented by calling the user subroutine UEXTERNALDB. This subroutine provides flexible points-in-time where other software or codes are able to communicate with the finite element analysis program (i.e., Abaqus/Standard), and data information stored in external files can be read in.

After user subroutine UEXTERNALDB is called, a user-supplied subroutine DDCODE is implemented to specify the DDD-relative calculations. It is constantly called at the beginning of each increment (i.e., LOP = 1) to update the state of dislocations, obstacles, and Frank–Read sources and to calculate the plastic strain induced by dislocation gliding.

There are ten self-defined subroutines, STRAINSUB, SOURCESUB, RFRESUB, ABAIDSUB, DISPSUB, NEWCODSUB, FORCESUB, OBSSUB, PNPLOY, and ANNISUB consisting of the dominant subroutine DDCODE. These subordinative subroutines implement the dislocation mechanisms by operating on the data structures of dislocations, dislocation sources, and obstacles. The data structures used in the model are introduced in the next section. The calling sequences of these subroutines are shown in Figure 5. The subroutine FORCESUB is used to calculate the force applied on dislocations using the stresses transferred from the FEM module, by Equation (1). Once the force on dislocations is obtained, their velocity and displacement can be given by Equation (2), as performed in subroutine DISPSUB. The transformation matrix and equations have been provided in Section 3.2. According to the displacement of dislocations, new coordinates of dislocations are calculated in subroutine NEWCODSUB. The subroutine OBSSUB determines if dislocations are pinned or unpinned by obstacles and it provides the dislocation numbers that are pinned. The subroutine ANNISUB decides whether the dislocation annihilates with other dislocations according to the relative positions of dislocations. The subroutine SOURCESUB, which follows the dislocation-generation rules discussed in Section 2.1, generates a dislocation dipole once all the conditions (i.e., the critical strength and nucleation time) are met. The subroutine REFSUB eliminates the dislocations that annihilate or escape from the simulation model and re-number the remaining and newly generated dislocations. The subroutine ABAIDSUB, which calls subroutine PNPLOY, identifies which element and integration point domain the dislocations belong to. The subroutine PNPLOY is called to find whether a point is inside or outside a polygon. The last one subroutine STRAINSUB calculates plastic strain induced by each moveable dislocation and distributes them to corresponding integration points by the FE mesh-based method introduced in Section 3.2.

The plastic strain tensors at each integration point are saved in a three-dimensional array. It is important to notice that all the subroutines used in the paper is with FORTRAN language as most user subroutine interfaces in the manual of ABAQUS use it. The FORTRAN language, unlike other programming languages, stores arrays in column-major format, as shown in Figure 6a. That is to say, in a two- or multi-dimensional array, the left-most array index varies the most rapidly and the right-most array index varies the most slowly. Thus, the index of the first dimension (from left to right) represents the components of plastic strain and the component order is the same as that in UMAT, as illustrated in Figure 6b. The second and third indices are the integration point and element number where the plastic strains are distributed, respectively. This storage means of plastic strain allows the generation of simple and efficient machine code, to improve the computing and running speed. The storage format of the multi-dimensional array and the illustration of the plastic strain array are shown in Figure 7.

Moreover, the subroutine UEXTERNALDB is called once at the start of the analysis (LOP = 0). This is to obtain the reference state for the DDD module. The DDD reference state consists of the initial information of dislocation, obstacle, and Frank–Read source, such as the coordinates, critical strength, and nucleation time, before loading is applied. The data structure of the reference state is given in Section 4.2. The initial information is stored in external files. At the beginning of the analysis, subroutine UEXTERNALDB reads data from these files and uses them to initialize the DD module.

The subroutine UMAT, allowing user-defined mechanical material behavior, implements the above physics-based constitutive law (Equation (1)) by using plastic strain calculated by DDD. When the subroutine is called, the parameters such as stresses, strain increments, and state variables (SDVs) are directly passed into UMAT, while the material parameters defined in the input file are accessed by the variable PROPS. With the variables passed in, the UMAT subroutine updates the stresses and SDVs at the end of the increment and generates the Jacobian matrix of the constitutive model. An important note is that the definition of the material Jacobian affects only the convergence rate, not the simulation results obtained [48].

4.2. Data Structure for the Coupling Framework

The data structures of dislocations, obstacles, and dislocation sources are shown in Figure 8. It can be seen that the basic data structure of a dislocation consists of identity number, location information, its sign (i.e., plus or minus one, representing positive or negative dislocations, respectively), and the state information. The information of location includes the coordinate value, the located slip system, and the element and integration point number where the dislocation resides, as shown in Figure 8a. The preliminary calculation with user subroutines generates the stress, displacement, and new coordinates of the dislocation $i$, and then the state information (i.e., illustrating whether the dislocation is pinned, annihilated, or newly generated) is obtained. Moreover, as can be seen from Figure 8b, the data structure related to dislocation sources and obstacles is much simpler compared to the dislocations. In addition to the location information, the critical strength and time are added for sources and only strength is added for obstacles, as shown in Figure 8c. During the deformation process, these data are read and updated continuously by the self-defined subroutines mentioned in the above section. As shown in Figure 8a, for the dislocation, its relative element and integration point number are updated by the subroutine ABAIDSUB, and the driving force, displacement, and new coordinates are calculated and updated by subroutines FORCESUB, DISPSUB, and NEWCODSUB, respectively. The signs of annihilating, newly generating, and pinning are determined by subroutines ANNISUB, SOURCESUB, and OBSSUB. The data structure of sources is mainly visited and updated by subroutine SOURCESUB. Similarly, the data structure of obstacles is accessed by subroutine OBSSUB.

The initial data structures of dislocations, sources, and obstacles (i.e., the data structures before loading is applied) are utilized to describe the reference state of the DDD module. They are generated by MATLAB and written into text files that are visited at the start of ABAQUS/Standard analysis to initialize the DDD module.

5. Uniaxial Compression Simulation for Single Crystal Micropillar

5.1. Computational Model of Micropillars

For verification, the developed simulation model is used to study the deformation behaviors of crystal materials under uniaxial compression. As shown in Figure 9, micropillars with four different widths are built with a constant aspect ratio of height H to width L (H/L = 2). An isotropic aluminum single crystal is considered here with Poisson’s ratio $v=0.3$, shear modulus $\mu =26 \mathrm{Gpa}$, Burgers vector $b=0.25 \mathrm{nm}$, and drag coefficient $B={10}^{-4} \mathrm{Pa} \mathrm{s}$ [28,44,45].

For face-centered crystal Al, the representative material parameters used in DDD are the same as the ones used in other research [28,30,43,45]. Two slip systems with the interaction angle of $\pm {60}^{°}$ relative to the horizontal direction are considered with the spacing distance of two neighboring slip planes being 10b. This approximates to the active slip systems of the face-centered crystal in plastic strain problems. The x-axis and y-axis represent the crystal directions [101] and [010], respectively. Static dislocation sources are distributed on these slip planes randomly with strength satisfying the Gaussian distribution with mean value 50 MPa and standard deviation 10 MPa. The critical time for dislocation nucleation is set as 10 ns. The obstacles in the model are with mean strength 150 MPa and 20% standard deviation.

To mimic the compression experiment, a fixed compressive strain rate $\dot{\epsilon }=1000/\mathrm{s}$ is imposed for all samples. Time steps of $1 \mathrm{ns}$ for DD and $10 \mathrm{ns}$ for the FEM are used here [49]. In all the simulations, the bottom boundaries of specimens are fixed as:

$$u_x=u_y=0, at Y=0,$$

(12)

The upper surface is subjected to a velocity $v_x$ such that

$$v_y=L\dot{\epsilon }, at Y=L,$$

(13)

DCM simulations are performed with the boundary conditions and material parameters as described above. The simulation results are discussed below.

5.2. Effect of Sample Size on Crystal Strength

The computationally obtained $\sigma ~\epsilon$ curves are plotted in Figure 10 for four different simple sizes $L=\left\{0.3, 0.5,0.8,1.0\right\} \mathsf{\mu }\mathrm{m}$. To avoid the potential effect of defect distribution (i.e., the distribution of dislocations, sources, and obstacles) on the simulation results, for each size, three simulations are carried out with different defect distributions. It is found that the smaller the sample size, the higher the strength, as observed in experimental investigation [5]. This validates that the present simulation model is capable of capturing the microscale plastic response. Due to numerous dislocations being greatly active at L = 1.0 μm, the stress–strain curves oscillate poorly and even fail at small strains. A similar phenomenon was also found in the compressing simulation of copper micropillars [50]. From these stress–strain curves, we compute the nominal yield stress ${\sigma }_{0.2}$, which is the stress value at 0.2% plastic strain offsets, and plot it in Figure 11a. The inset of Figure 11a shows the acquiring of the normal yield stress for the 0.3 μm-wide micropillar. The results indicate that the yield stress increases with decreasing pillar size. The error bars denote the range of the values of ${\sigma }_{0.2}$.

In addition, the secant strain hardening modulus, $H_{sec}$, is defined as the ratio of stress to plastic strain from ${\epsilon }_p=0.1\%$ to ${\epsilon }_p=0.3\%$, as follows:

$$H_{sec}=\frac{{\sigma }_{\epsilon =0.003}-{\sigma }_{\epsilon =0.001}}{0.003-0.001},$$

(14)

Then, the average value of $H_{sec}$ over three realizations for a particular pillar size, $H_{ave}$, is calculated and plotted in Figure 11b. It shows that the average strain hardening modulus decreases with sample size and that larger pillars are with the smaller scatter of $H_{ave}$ values. In addition, the ratio of surface area to volume, which is reduced to the ratio of length to surface area for the plane strain problem, is calculated for each pillar. The inset of Figure 11b shows that the ratio decreases with sample size, although all the samples have the same length–width ratio. It indicates that dislocations in small samples have a high possibility to escape from the sample. The escape of dislocations from the free surface can result in the dislocation starvation effect [49] in small pillars and the hardening of materials.

6. Conclusions

In this paper, a concurrent multiscale simulation method is developed to shed light on the material plasticity at the microscale. In this method, the DDD module looks to the FEM module to solve the boundary conditions and the FEM module turns to DDD for plastic strain, which determines the constitutive law of the FEM.

The key issue in the DDD-FEM coupling scheme is how to reasonably regularize the discrete plastic strain to integration points of continuum mechanics. A novel FE mesh-based regularization method is proposed by using a Burgers vector density function and the characteristics of a finite element mesh in ABAQUS. It shows high accuracy and considerably improves the efficiency by avoiding the tedious calculations.

To achieve the concurrent computational model, the DDD module is programmed into Fortran codes, and then it is incorporated into the Abaqus analysis by calling user subroutines provided by Abaqus to determine the constitutive law. The stress field is stored in COMMON blocks and transferred to DDD codes. The realization of the multiscale model capitalizes on the secondary development of FEM software Abaqus, and the implementation details are presented for the first time.

The developed multiscale framework is able to capture the plastic behavior of microscale crystals. Using this multiscale method, uniaxial compression tests of micropillars are performed to study the size effect on plasticity and its underlying mechanisms. The observed secant strain hardening of the modulus and yield stress might be due to the high probability of dislocations in smaller size to escape from the free surface.

This multiscale framework can be potentially employed in studying the plastic behavior of materials with complex microstructures and deeply revealing its underlying mechanisms. In future work, the developed simulation model will consider a variety of complex dislocation motions, such as dislocation climb, interactions between dislocation and grain boundary, and dislocation junction formation.