Green's function molecular dynamics Including finite heights, shear, and body fields

The Green ’ s function molecular dynamics ( GFMD ) method for the simulation of incompressible solids under normal loading is extended in several ways: shear is added to the GFMD continuum formulation and Poisson numbers as well as the heights of the deformed body can now be chosen at will. In addition, we give the full stress tensor inside the deformed body. We validate our generalizations by comparing our analytical and GFMD results to calculations based on the ﬁ nite-element method ( FEM ) and full molecular dynamics simulations. For the investigated systems we observe a signi ﬁ cant speed-up of GFMD compared to FEM. While calculation and proof of concept were conducted in two-dimensions only, the methodology can be extended to the three-dimensional case in a straightforward fashion.


PAPER • OPEN ACCESS
Green's function molecular dynamics: including finite heights, shear, and body fields

Introduction
Green's function molecular dynamics (GFMD) [1][2][3] is a boundary-value method allowing one to simulate the linear-elastic response of a solid to an external stress or, more generally, to a boundary condition. So far, GFMD has been used predominantly to describe either nonreflecting [4,5] and thermalizing [6][7][8] boundaries to which an atomistic region is coupled, or, as a tool to simulate the contact mechanics of solids with rough surfaces [9][10][11]. One advantage of GFMD is that it only necessitates knowledge of the displacements in the top layer of a solid and that effective interactions are block diagonal in Fourier space. Relatively large systems can therefore be simulated and be quickly relaxed. Typical system sizes in the context of contact mechanics range from 4096 × 4096 surface atoms on single CPUs [9] to O 10 10 5 5 ( )on supercomputers [11]. An additional, conceivable application consists in coupling GFMD to discrete dislocation dynamics (DDD) [12]. The idea is to use GFMD, instead of the finite-element method, to compute the image fields of dislocations in DDD. Towards this end, we generalize the GFMD method in the following ways: first, we consider the elastic response of a cubic or an isotropic body with arbitrary Poisson number and allow for lateral displacements as well as for shear tractions in addition to normal tractions. Second, we deduce the internal stresses for a given surface boundary condition and do this for solids of arbitrary height. The approaches pursued so far were limited to either normal displacements and normal tractions in the continuum formulation [11,13] or to the full atomistic Green's functions [2,9], which do not relate directly to the continuum limit. While the finite-width elastic continuum problem with shear was solved by Carbone and Mangialardi [14], their work did not put us into the position to deduce directly the Green's function coefficients needed for a numerical implementation. In a later work, in particular appendixA of [15], useful formulae for the GFMD simulations were stated, but unfortunately only for the frictionless case.
In this work, we present a solution for the Green's function of finite-height elastic slabs having the following advantages: the only required mathematical tools are partial derivatives, Fourier transforms, and linear algebra, i.e., there is no need to solve Fredholm integral equations. All equations needed to implement the approach into computer code are given explicitly in compact form. Moreover, our approach can be readily modified in various ways. For example, it should be straightforward to extend our solution strategy to layered materials, to materials with gradient or square-gradient corrections to the elastic energy, or to nonisotropic crystals-as long as they remain homogeneous within each plane. In fact, the most important equations for non-isotropic crystals are given and tested in this work. Lastly, we validate our solution against numerical data and moreover consider various limiting cases including that of very small slab heights or that of a vanishing shear modulus characteristic for fluids.
The coupling between GFMD and DDD will be presented in a separate manuscript.

General background
In this paper, we are concerned with the quasi-static loading, in which case the precise dynamical response of the simulated layer does not play a role, see [5,15,16] for generalizations from the static to the dynamic case. Moreover, we consider to load the surface of a body that is translationally invariant within the xy-plane, which, in principle, is allowed to be a gradient material as long as the gradients are normal to the xy surface plane. One can then write the linear stress-displacement relation u r s where G r¢ ab ( ) is the Green's function tensor, u r a ( ) is the α component of the displacement as a function of the (two-dimensional) in-plane coordinate r, and 3 s b is the traction in zdirection. In Fourier space, equation (1) reads In contact problems, one often knows the displacement of the bodies and wants to deduce the contact pressure, and thus, one usually does not need to evaluate the Green's function (tensor) itself, but its inverse. Thus, force calculations in GFMD simulations require one to evaluate The precise functional form of the (inverse) Green's function tensor depends on the elastic properties of the deformed material including its height. In its simplest form, usually used in the context of (continuum) contact mechanics [17], one is only interested in normal displacements induced by normal tractions applied to a semiinfinite, isotropic body. In this case, if the body is incompressible, 0.5 n = , all quantities in equation (3) can be considered as scalars and the equation reduces to )is the contact modulus, E being the Young's modulus and ν the Poisson number. However, as mentioned in the introduction, we wish to generalize GFMD to give the elastic response of a body with generic Poisson's ratio, and therefore, we can no longer use a scalar to describe surface displacement. Moreover, we intend to consider problems where contact loading is not restricted to be in normal direction, i.e. tractions and/or displacement can be applied in normal or tangential directions, everywhere on the body surface. Therefore, we can no longer rely on equation (4) and need to find a form for the Green's function tensor in equation (3), which has to depend on the Poisson's ratio and on the height of the slab. Towards this end, we next calculate the analytical solutions for the displacement in linearly elastic slabs of finite height, from which the stresses can be deduced in a straightforward fashion. The stress distribution underneath the contact are of particular interest in problems as fretting fatigue, to determine whether a possible tensile loading underneath the contact would give rise to crack nucleation and propagation.

Analytical solutions for the displacement in finite-height, linearly elastic slabs
We consider a linearly elastic body of cubic or higher symmetry in a slab geometry with a fixed bottom, i.e., the displacement reads x z u , ) as for semi-infinite solids. Moreover, we assume that no body forces are exerted, which implies the usual equilibrium condition , where r s ab ( ) is the stress at the point r inside the body and r ¶ º ¶ ¶ a a . For isotropic or cubic bodies the equilibrium condition is given by where we have restricted our attention to (1+1)-dimensional solids so that in-plane wavevectors are now scalars, and where the C ij ʼs denote coefficients of the elastic tensor in Voigt notation.
Assuming an in-plane undulation of the top layer with the real-valued wavenumber q, equations (5) and (6) can be solved with the factorization Thus, we obtain solutions for the displacements either oscillating exponential functions for b 1 < or purely exponential functions for b 1 > . The nature of the solution changes at b=1, which automatically holds for isotropic media as these satisfy the isotropy condition C C C 2 ) . The solutions for b=1 are proportional to qz exp  ( ), and, in addition, proportional to z qz exp  ( ), i.e., similar to those of critically damped harmonic oscillators. The decaying solutions can usually be ignored when the z-position of the top layer z m satisfies z q 1 m  but not for a finite-slab geometry. In the remainder of this section, we focus on the isotropic case, because this is a common approximation. In the result section, we also consider the case of a cubic solid violating the isotropy condition to demonstrate the correctness of our approach. At this point, it may suffice to state that metallic cubic crystals tend to have a relatively small shear modulus, in which case b 1 > , while non-metals rather correspond to b 1 < . Due to the nature of the differential equations (5) and (6), the solutions of the in-plane cosine transform of the lateral u 1 displacement field couples to the in-plane sine transform of the normal u 3 displacement, and vice versa. Thus, we can write Solutions satisfying the boundary condition x u , 0 0 = ( ) and the differential equation for isotropic media are then obtained after some algebra to satisfy where s C C 44 11 º . The latter ratio has allowed values of s 0 1 < < for stable, twodimensional isotropic solids [18]. The coefficients A 1,2 follow from equation (12) In summary, for an arbitrary surface displacement field x z u , m ( ), the in-plane Fourier transform is taken to yield q z u , m ( ). The real and imaginary parts can be associated with lefthand sides of equations (12) and (16), which allow one to determine the pertinent coefficients A 1,2 and B 1,2 for each wavenumber q by evaluating them at z z m = . The knowledge of these coefficients then allows one to deduce the displacement inside the body.
2.3. Finite-height-slab strain, stress, and energy density Not only the displacement but also the strain and thus the stress field on the surface or inside the body can be deduced as soon as the coefficients A 1,2 and B 1,2 have been determined for a given surface topography, obtained by perturbation of an initially flat surface. For reasons of simplicity, we first restrict our attention to the case of a perturbation by a single wave number q and B 0 1,2 = . The elements of the infinitesimal Cauchy's strain tensor (in Voigt notation) are then given by      ) , which indeed can be shown to be half E * . This implies that equation (36) is consistent with (4) in the limit of a frictionless semi-infinite contact.

Displacements in non-isotropic solids
To model a non-isotropic solid, we consider a cubic crystal with its (100) surface facing up. To make the comparison of continuum theory to full molecular dynamics (MDs) simulations as transparent as possible, we furthermore restrict ourselves to a simple cubic lattice in which each atom (which one may also see as a grid point) is connected to its nearest neighbors with a spring of stiffness k 1 and to next-nearest neighors with a spring of stiffness k k 2 1 = . Mechanical equilibrium of the springs is assumed at a distance a 0 between nearest neighbors and a 2 0 for next-nearest neighbors. With these definitions, it is readily seen that the elastic ) . Here, the stiffness is divided by a 0 so that the elastic constants have the usual units.
In our example, we consider a slab of height L L 2 is the lateral length of the domain, which is repeated periodically along the x-direction. The just-defined system is solved with a self-written MD code, in which individual atoms are also coupled to damping linear in velocity. The mass of atoms is set to unity, the time step to 0.1 and damping to 0.25. Two cases of displacements in the top layer are treated: normal loading for which u x z , 0 shear loading for which  Agreement between full MD and the analytical expressions for the displacement fields is clearly revealed in figure 1.

Displacement and stress fields in isotropic solids
As a benchmark problem to compare GFMD to FEM we here consider the indentation of an isotropic elastic two-dimensional slab by an array of flat rigid punches. Contact between punches and slab is taken to be fully sticking. The analysis is performed on a periodic unit cell with fixed bottom as shown schematically in figure 2(a).
Normal displacement is prescribed at the contact of length L x p : The slab is indented to u L 2.5 10 . For the finite-element analysis the slab is discretized using a uniform mesh of square elements. The number of degrees of freedom is n nnx nnz 2 dof =´, where nnx is the number of nodes in x-direction, and nnz the number of nodes in z-direction. For the GFMD simulation the surface is discretized using nx equispaced grid points, with nx=nnx. Contact between the rigid indenters and the slab is modeled through a hard-wall potential. The static solution is found in GFMD using damped dynamics as described in [11]. The damping factor η used in this simulations is  Following [19], the error in the norm obtained using FEM as a function of the degrees of freedom in the simulation can be written as: The exact solution, the asymptotic rate of convergence β, and the prefactor κ, are obtained by linearly fitting three data points corresponding to nnx=256, 512, 1024. As expected [19], given that the mesh refinement is uniform and the order of the interpolating polynomial for the shape functions is one, the asymptotic rate of convergence is found to be 0. x L 1 0.5 which allows for direct comparison with GFMD (see figure 5(a)). The order of convergence with respect to the discretization of the two methods is found to be the same, while the prefactors are favorable for GFMD, i.e. L 3.609 10 Figure 5(b) shows the simulation time as a function of the surface discretization. The results are all obtained on a single Intel Xeon(R) 3.10 GHz processor with 31.3 Gbytes of RAM. The GFMD simulations are found to be faster than FEM, and the computational advantage increases with increasing system size. In addition to this, a smaller number of grid points are needed in GFMD to obtain the same results as in FEM. For the simulation reported here, if we decide to tolerate an error e=0.005 in the L 2 norm of lateral surface displacement, the GFMD simulations require nx=128 grid points while the FEM simulations need nnx=1024 surface nodes. This results in a GFMD simulation being 1650 times faster than FEM.
It is to be noted here that: (1) we have not searched for the optimal meshing scheme to solve the finite-element problem, simply used a mesh with squared elements for ease of comparison with the equi-spaced surface grid points used in GFMD; (2) the speed of the finite-element simulation depends heavily on the solver used. Here we have used a direct sparse solver with skyline storage, where the time consuming step is the factorization of the skymatrix. The order of factorization in a skyline solver is generally O nnx B 2 (( ) ) [20] where B is the mean bandwidth, which cannot exceed nnx. We are aware that for large systems, an iterative solver would be much more cost-effective. In GFMD, the computational complexity scales with O nx nx nx log ( ), (3) the speed of the GFMD simulation depends on the choice of the damping factors, which in general are different in the x and z-direction. The optimal damping factor to obtain critical damping of the system depends on the loading, the height of the slab and the elastic constants. We here simply considered a single damping factor that would critically damp the modes in z-direction, and thus under-damp the modes in xdirection.

Conclusions
GFMD, a fastly converging boundary value method used to compute contact pressures and surface displacements of incompressible continuum semi-infinite solids, is here extended to apply to finite solids with generic Poisson's ratio and boundary conditions. Moreover, the body fields can now be computed analytically from the tractions and/or surface displacements. This extension allows the GFMD technique to provide the same information that can be obtained through the FEM, but with a significant gain in simulation time.
An additional advantage of GFMD is that for typical contact problems, where the contact area evolves during the simulation, the contact can be easily captured and simply described by means of an interaction potential. We have here used a hard-wall potential, but one can also model the bodies in contact explicitly and apply an interaction potential, as the Xu-Needleman [21], where the contact response in normal direction is coupled to that in tangential direction.
An interesting application that can be envisaged for the GFMD method, in virtue of the extensions presented in this paper, is the replacement of the FEM in discrete dislocation plasticity simulations of contact. This has the potential to significantly increase the time efficiency of the discrete dislocation plasticity calculations by that allowing to extend the applicability of such models to bodies of larger size, and with a realistic surface profile.