Molecular statics simulation of nanoindentation using adaptive quasicontinuum method

In this work, molecular statics is used to model a nanoindentation test on a two-dimensional hexagonal lattice. To this end, the QuasiContinuum (QC) method with adaptive propagation of the fully resolved domain is used to reduce the computational cost required by the full atomistic model. Three different adaptive mesh refinement criteria are introduced and tested, based on: (i) the Zienkiewicz– Zhu criterion (used for the deformation gradient), (ii) local atoms’ site energy, and (iii) local lattice disregistry. Accuracy and efficiency of individual refinement schemes are compared against the full atomistic model and obtained results are discussed.


Introduction
Nanoindentation is a commonly used testing procedure applied to small volumes of materials for measuring their micromechanical properties. Typically, a hard tip (i.e. indenter) with known mechanical properties is pressed into an examined sample of unknown mechanical properties. Loading force and penetration depth of the indenter are recorded during the loading and unloading stages, providing a basis for the estimation of the unknown mechanical properties.
Numerical models are typically used as a tool for better understanding the underlying phenomena, and to obtain detailed information about local mechanisms occurring below the indenter tip (such as dislocation nucleation, propagation, and interaction), which directly influence measured reaction force. To this end, both the indenter and specimen are typically modelled at the atomistic level using molecular statics or molecular dynamics, entailing high computational costs when realistic configurations and dimensions are used. The QuasiContinuum (QC) method (cf. e.g. [1]) is employed to simplify the full atomistic model, to reduce the associated computational costs, and to allow for modelling of realistic situations. This paper focuses on the predictive abilities of an adaptive QC methodology (recalled in Section 3) in combination with three types of error indicators/estimators for local mesh refinement compared against the full atomistic simulations. In particular, (i) the Zienkiewicz-Zhu error estimator (used for the deformation gradient), (ii) an indicator based on local atoms' site energy, and (iii) an estimator based on local disregistry profiles are tested for a simple twodimensional indentation test. The individual definitions are outlined in Section 3.3, whereas the accuracy and associated computational costs are discussed in Section 4.
The paper closes with conclusions and recommendations in Section 5.

Full atomistic model
Atomistic models based on molecular statics are characterized by an underlying lattice in combination with an interatomic potential. In this work, a twodimensional hexagonal lattice with lattice spacing d 0 is used, as shown in Fig. 1. Individual atom interactions are described by the Lennard-Jones (LJ) potential, defined as where r αβ = r β −r α 2 denotes the distance between two atoms α and β, r m denotes the distance at which the interaction energy reaches its minimum, and ε is the energy well depth. The total potential energy associated with the entire atomic structure is computed as a sum over all interactions, i.e.
where N Atm represents the number of atoms, and r is a column storing their positions. Because evaluation of the interatomic potential for all pairwise combinations is computationally expensive, and because long-distance interactions have negligible contributions to the total potential energy, a cut-off radius r cut is considered [1], beyond which interactions are neglected. Such a simplification introduces a discontinuity of φ αβ at r cut , which is removed by subtracting a linear function to assure zero value and zero slope of φ αβ at r cut . As sketched in Fig. 1 cut-off radius r cut = 2.5 d 0 is employed to provide next-to-nearest interactions.
In order to find a stress-free configuration, an initial relaxation is carried out on an ideal periodic lattice with spacing d 0 = r m , which results in a reduced lattice spacing d 0 = 0.9917496 r m used for constructing the initial system.
The geometry of employed indentation test is sketched in Fig. 2. The specimen domain is of the size 128d 0 × 128h 0 , contains 16, 862 atoms, and considers atoms near the bottom and both vertical edges as fixed, whereas the top edge is a free surface. The flat indenter is modelled at the atomistic level using the same hexagonal lattice as used for the specimen, but having infinite stiffness. Its geometry is specified through a width 11d 0 at the tip and two surfaces inclined by 60 • , as shown in Fig. 3 (left). The positions of all indenter atoms are prescribed in 80 uniform loading, and 80 uniform unloading steps, achieving the maximum indentation depth 8d 0 . The interaction strength between atoms of the indenter and the tested material is reduced by a factor of 0.55 (compared to the atoms of the tested material) to prevent tearing of the indented specimen during the unloading stage. The total potential energy of the atomistic system E(r) is minimized at each time step using the trust-region algorithm; for further details see e.g. [2].

Quasicontinuum Method
The Quasicontinuum (QC) method is a concurrent multiscale technique introduced in [3]. The key idea consists in combining the accurate but expensive atomistic description only in regions of high interest with a cheap continuum approximation elsewhere. The specimen domain therefore is divided into two parts: (i) the fully-resolved region, in which the full nonlocal atomistic model is used, seamlessly coupled with (ii) a coarse-grained continuum region, in which interpolation through triangular elements along with an efficient summation scheme is introduced.

Interpolation
The first step in the QC reduction is interpolation, which introduces the so-called repatoms through which the kinematic behaviour of the entire system is reconstructed according to where r rep is a column storing positions of all repatoms, interpolated through an interpolation matrix Φ associated with the adopted triangulation. The number of repatoms is typically much smaller compared to the number of all atoms, reducing thus the computational effort required. Two example triangulations associated with different sizes of the fullyresolved regions are shown in Fig. 4.

Summation
In the second QC reduction step, a so-called summation rule is introduced to avoid the necessity of visiting all atoms when assembling the total potential energy in Eq. (2). To this end, the site energies of all atoms situated inside a triangular element are approximated by the energy of only a few, or even one sampling atom and its corresponding weight factor w α , i.e.
vol. / In Eq. (4), φ α is the site energy of a sampling atom α, defined as

Molecular Statics Simulation of Nanoindentation
Although multiple possibilities are available in the literature to select sampling atoms and their corresponding weight factors, the central summation rule introduced in [4] is used hereafter (cf. Fig. 4).

Adaptivity
The area of high interest (i.e. the fully-resolved region and hence also the associated triangulation) can adaptively evolve at each time increment or iteration to accommodate dislocation movements, while retaining QC efficiency.
To this end, the fully-resolved region is sequentially updated as follows. At each time increment, the system is equilibrated for a fixed fully-resolved region. Using a selected refinement criterion (detailed below), the triangulation is checked and refined if required. Any additional atoms are added as repatoms, the interpolation mesh is updated, and equilibrium is restored. Such a procedure is repeated until the mesh refinement criterion is satisfied for all elements, proceeding subsequently to a new load increment.
In total three refinement criteria are tested: (i) The Zienkiewicz-Zhu error estimator (ZZ), as introduced in [5], used for the deformation gradient. That is, projection of the deformation gradient is used to estimate the local error inside each element. A threshold value ZZ tr is specified to determine elements to be refined, and atoms located in those elements to become repatoms.
(ii) The error indicator based on a local energy criterion uses the local atoms' site energy to determine elements to be refined. First, the site energy of each sampling atom is tested for the following condition where E tr is a selected critical energy threshold. Sampling atoms satisfying the condition of Eq. (6), as well as their neighbours within radius r ref , are labelled as critical atoms. Interpolation elements that contain at least one critical atom are refined, i.e. all atoms in such elements are added as repatoms.
(iii) The error indicator based on disregistry profile works in a similar way as the local energy criterion. The only difference consists in the condition for the selection of critical atoms, now specified as where D tr is a selected disregistry threshold, and ||r β − r α t || 2 denotes the distance of atom β from the closest position of its theoretical neighbouring atom α considered in the reference configuration. The positions of either all eighteen nearestand next-to-nearest-neighbour atoms or just the six nearest-neighbour atoms can be used to evaluate the condition of Eq. (7). All atoms within radius r ref are again labelled as critical atoms.

Results
Numerical simulations of the nanoindentation test, as described in Section 2, are carried out to analyse the adaptive propagation of the fully-refined region associated with movements of individual dislocations.
In the full atomistic model, initially all atoms are displaced elastically until a critical penetration depth is reached, which triggers nucleation of four dislocations positioned symmetrically (due to the symmetry of the problem) under both indenter edges. Upon further loading, these dislocations propagate along preferred lattice directions until reaching the boundary of the specimen. During unloading, however, when the last symmetric dislocation pair annihilates, one direction is preferred due to numerical round-off errors, resulting in a nonsymmetric final shape of the indent (see Fig. 3 (right)). Results of the full simulation are used as the reference solution against which the QC simulations are compared, and are summarized along with the corresponding performance in Table 1. The force-displacement diagram is shown in Fig. 6.
In QC simulations, movements of individual dislocations can be captured properly only inside the fully-resolved region. The adaptive algorithms introduced in Section 3.3 are used to this end, providing results summarized in Fig. 5 for the Zienkiewicz-Zhu error estimator and loading steps number 28 (penetration of the first symmetric dislocation pair outside of the original area of high interest), 39 (one step before nucleation of the second symmetric dislocation pair), 55 (one step before the nucleation of the third symmetric dislocation pair), and 80 (maximum penetration depth). The following cases for error estimators/indicators in QC methodology have been tested:    (1) two simulations with fixed fully-resolved regions of different sizes, see Fig. 4, referred to as QC small and QC large;

Initial/final Iterations Mesh
(2) one simulation with an adaptive area of high interest using the ZZ criterion with ZZ tr = 0.015, referred to as QC zz; (3) one simulation with an adaptive area of high interest based on local energy with refinement radius r ref = 5d 0 and threshold value corresponding to 6% change of the initial potential energy (for a typical internal atom this value corresponds to E tr = −2.95); and (4) one simulation with adaptive area of high interest based on lattice disregistry with threshold value D tr = 0.5d 0 , the nearest-neighbour option, and refinement radius r ref = 5d 0 , referred to as QC dis.
The simplest QC simulation with a small fixed fullyresolved area (QC small) provides a significant speedup of a factor of 16 compared to the full atomistic model. The initial elastic behaviour and nucleation of the first symmetric dislocation pair is captured accurately, cf. Fig. 6, whereas in later stages (i.e. once the first symmetric dislocation pair reaches the boundary of the fully-resolved region), the resulting behaviour becomes overly stiff. This is clearly visible in the corresponding force-displacement diagram.
A partial improvement can be achieved by enlarging the fully-resolved QC region (QC large), which, nevertheless, suffers from the same shortcomings. The overestimation of the force-displacement diagram is less significant, capturing more accurately the unloading branch. The shapes of local force peaks do not correspond, however, to the exact solution due to the obstructed dislocation movements. The achieved speed-up is approximately of the order of 10.
Clearly, only the adaptive schemes are able to capture properly the dislocation movement under the tip of the indenter. Note that all adaptive QC examples are initialized with the fully-resolved area of the same size as used for the QC small simulation. The ZZ criterion detects and refines elements near individual dislocations, allowing for their propagation throughout the entire specimen. This is reflected in the corresponding force-displacement curve in Fig. 6, where a significant improvement in accuracy compared to the approaches with fixed meshes can be observed.
Finally, the local lattice deflection criterion is used, which provides (in the case considered) practically the same results as the energy criterion (not shown). Generally, both criteria are not equivalent and differences can be observed in examples with different type of atoms, where the energy criterion is more robust at the interface. However, in case of a single crystal with one type of atom, for each value of disregistry threshold a corresponding value of energy threshold that provides virtually the same results can be found. Compared to the ZZ condition, local lattice deflection provides similar accuracy in the resulting reaction force, while requiring a slightly higher number of repatoms (cf. Fig. 7). In spite of this fact, its overall performance is slightly faster compared to ZZ (cf. Tab. 1), explained by the fewer number of mesh iterations needed. Unlike ZZ, the indicator based on the local disregistry provides two parameters, D tr and r ref , for controlling the mesh evolution. Whereas D tr locates the position inside the dislocation core region (see colours in Fig. 5), r ref reflects the size of the fully-resolved region considered around that position. Large values of r ref therefore provide fast propagation in the direction of the dislocation movement, but at the same time may require too many repatoms in the perpendicular direction. The value r ref = 5d 0 proved to provide a good balance.
The performance of all QC approaches is summarized in Tab. 1, where the number of repatoms, required number of solver and mesh iterations, and computational times are reported. All adaptive QC approaches provide accurate descriptions during the loading as well as unloading stages, resulting in acceptable estimates of the corresponding reaction force. Although the obtained gain in computation time is relatively small (about 25% compared to the full simulation), a more substantial speed-up is expected for