Green's function method for dynamic contact calculations

Resolving atomic scale details while capturing long-range elastic deformation is the principal difficulty when solving contact mechanics problems with computer simulations. Fully atomistic simulations must consider large blocks of atoms to support long-wavelength deformation modes, meaning that most atoms are far removed from the region of interest. Building on earlier methods that used elastic surface Green's functions to compute static substrate deformation, we present a numerically efficient dynamic Green's function technique to treat realistic, time-evolving, elastic solids. Our method solves substrate dynamics in reciprocal space and utilizes precomputed Green's functions that exactly reproduce elastic interactions without retaining the atomic degrees of freedom in the bulk. We invoke physical insights to determine the necessary number of explicit substrate layers required to capture the attenuation of subsurface waves as a function of surface wave vector. We observe that truncating substrate dynamics at depths that fall as a power of wave vector allows us to accurately model wave propagation without implementing arbitrary damping. The framework we have developed substantially accelerates molecular dynamics simulations of large elastic substrates. We apply the method to single asperity contact, impact, and sliding friction problems and present our preliminary findings.

Figure 1

Decomposition of DCGF method simulations. The topmost region is the MD region treated explicitly. The boundary region must be thicker than $r_c$ (gray, shaded circle) to account for all interactions with explicit atoms. Below the boundary, the total depth $N_{\text{tot}}$ is divided into $N_{\text{dyn}}$ layers treated using the DCGF method and $N_{\text{stat}}$ layers treated using CGF. The precise decomposition of $N_{\text{tot}}$ depends upon $q$.

Figure 2

Comparison of DCGF method and atomic computational time per time step for systems with $L\times L$ atoms in the boundary layer. The fully atomic system contains $L^3$ atoms and the computational time per time step scales similarly as $L^3$ (dashed line). DCGF method systems contain $L^2$ ($\mathrm{DCGF}+0$) and $9L^2$ ($\mathrm{DCGF}+8$) atoms in the MD domain. The corresponding computational time per time step scales as $L^{2}lnL$ (solid line).

Figure 3

Depth of the first layer where the kinetic energy of a pulsed surface wave drops below $e^{-3}\approx 0.05$ of the initial surface value as a function of wavelength $\lambda$ and damping coefficients indicated in the legend. The substrate is fully atomic with Kelvin damping between neighboring atoms, and the boundary layer contains at least $(\lambda /d_{nn})\times (\lambda /d_{nn})$ lattice sites. The solid line indicates attenuation depth equal to wavelength corresponding to Saint-Venant's principle (see text).

Figure 4

Center of mass displacement of a rigid disordered sphere under fixed load as a function of time for different DCGF and atomic systems. The reference system is fully atomic (black) with $128\times 128\times 129$ mobile lattice sites. $z_i$ is fixed across all simulations and $z_f$ is defined based on the fully atomic system. The denominator $z_f-z_i=-0.54\sigma$ under an applied load of magnitude $5\epsilon /\sigma$. A truncated atomic system (gray) with $128\times 128\times 33$ mobile lattice sites is included to emphasize the importance of depth on both dynamic and static properties. Agreement with the atomic result is incrementally improved by stacking additional lattice planes on top of the DCGF boundary layer. DCGF method systems contain $128\times 128$ lattice sites in the boundary layer with the number of additional lattice planes indicated in the legend.

Figure 5

Fractional error in the coefficient of restitution for the DCGF method compared to a fully atomic system with $128\times 128\times 129$ mobile lattice sites. DCGF method systems contain $128\times 128$ lattice sites in the boundary layer with the number of additional lattice planes indicated in the legend. The system is a rigid disordered sphere with variable initial velocity perpendicular to the boundary. The sphere rebounds away after contact. Bond length changes are measured in the atomic system, and always corresponded to the most compressed bond.

Figure 6

Kinetic friction per unit length imposed by the 1D sinusoidal potential moving along $\widehat{x}$ as a function of wavelength. The friction is normalized by the damping parameter and sliding speed. (a) Kinetic friction on 3D substrates with sliding speeds indicated in the legend. The solid line indicates $F_k\propto {\lambda }^{-1}$. (b) Kinetic friction on effectively 2D substrates with sliding speed $v=0.01d_{nn}/\tau$. The solid line indicates $F_k\propto {\lambda }^{-2}$.