Green's function molecular dynamics meets discrete dislocation plasticity

Modeling the contact mechanics of solid bodies assuming realistic surface roughness is highly relevant to tribology, the science of friction. This is a demanding task, because the height spectra of most surfaces, even polished ones, have non-negligible magnitude over several decades of wavelengths, typically from the atomic scale to several dozen or even hundreds of microns [3, 4], i.e., the root mean square heights are determined by the longest-wavelength height fluctuations (relevant for sealing and lubricant flow) while the root mean square gradients (relevant for typical contact stresses, at least in the elastic limit) are determined by the shortest wavelengths. In the past years, various theories of rough surface contact have been presented [5–10]. Persson's contact-mechanics theory [8–10] appears particularly promising to us, because it accounts for long-range elastic deformations, unlike traditional approaches such as those inspired by Greenwood and Williamson [6], who pursued local bearing models. In fact, Persson theory reproduces quite well experiments and numerical results on relative contact area, interfacial stress distribution functions, stress spectra, and gap distribution functions, although it requires a fitting parameter of order unity. However, the validity of Persson's analysis of plastic contacts [9] in the range where plasticity is size-dependent [11–14] has not yet been shown.

Numerical simulations of contact between elastic rough surfaces are made possible through several techniques. Of particular interest are the boundary-element based approaches such as BEM [15–17] and Green's function molecular dynamics (GFMD) [18–22], which calculate the response of a solid to contact loading by modeling only the surface. These methods are suitable to study large systems where the surface roughness is described by many orders of length scale [2, 20, 23]. In this work, we choose GFMD over BEM, since it is a simpler method, which does not involve solving the Fredholm integral equations. In GFMD, the equilibrium positions of the interfacial grid points are found by means of damped dynamics in Fourier space, where the individual modes decouple. This allows for large systems to be quickly brought to equilibrium. While the computational complexity of BEM scales with $O({n}^{2})$, where n is the number of discretization nodes, GFMD scales as $O(n\sqrt{n}\mathrm{log}(n)$. Furthermore, it is straightforwad in GFMD to employ non-holonomic boundary conditions and/or interaction potentials between the contacting bodies and hence it is cost-effective in studying problems where the contact area is not known a priori. Conventional FEM or BEM methods typically require several iterations as well as incremental updating of the boundary conditions in order to converge to a final contact area. Additionally, GFMD has the advantage that it can be extended to a multi-scale method [24] where the surface layer can be described atomistically and the substrate underneath be treated within the harmonic approximation. Pastewka et al [24, 25] have shown that the elastic Green's function can be quickly computed from interatomic potentials and a seamless coupling to the atomic region can be derived even for interactions going beyond nearest-neighbor interactions. GFMD has been successfully used to model the contact response of semi-infinite elastic solids [2, 26, 27], while plasticity was neglected.

There has also been much progress in the numerical study of elasto-plastic contacts. These, however, are either based on continuum plasticity [28–30], or, limited to single-wavelength roughness [31], or, based on brute-force all-atom approaches [32, 33]. To date, no studies of rough surface contact appear to have disseminated, in which (size-dependent) plasticity and long-range elasticity were both accurately modeled. Our work aims at building such a model.

Size-dependent plasticity in metal crystals has been successfully predicted by discrete dislocation plasticity (DDP) [1] simulations, which for simple problems, as the tensile response of free-standing metal films, give results in quantitative agreement with experiments [34]. In the context of contact mechanics, DDP has been used to study indentation of flat surfaced single crystals with single indenters [35, 36], periodic indenters [37, 38], as well as an indenter with self-affine surface modeled as a collection of Hertzian contacts [39]. The only DDP studies, in which the plastically deforming bodies were not approximated to be flat, involved the flattening of simple sinusoidal surfaces [40–42]. Results showed that the contact area is not continuous, but serrated, as a consequence of dislocations exiting the surface through discrete slip planes and leaving crystallographic steps. Due to the serrated nature of the contact area, the local contact pressure presented high spikes, at odds with the pressure profiles predicted by continuum plasticity. The study of size-dependent plasticity for realistic surface geometries is computationally very expensive. In order to being able to study indentation using realistic surface geometries, we here combine the accurate description of plasticity offered by DDP with the fastly converging elastic solution delivered by GFMD, in a modeling technique which we name Green's function dislocation dynamics (GFDD). This is not the first attempt to combine dislocation dynamics with a boundary element method [43–46]. However, for realistic surface geometries coming into contact, which require a fine discretization, GFDD should be more cost-effective. This is because the method relies on fast Fourier transform (FFT) and employs damped dynamics to quickly equilibrate large systems.

We begin the paper by briefly introducing DDP in section 2 and GFMD in section 3. We then present the methodology of the new model in section 4. The results obtained using the new GFDD model are compared with DDP in sections 5 and 6. Section 7 summarizes the advantages and potentials of the new method.

DDP is a numerical technique to solve boundary-value problems (b.v.p), which treats plasticity as the collective motion of discrete dislocations [1]. The dislocations are modeled as line defects in an elastic continuum. The solution at each time step of the simulation is obtained by the superposition of two linear elastic solutions: the elastic fields for dislocations in a homogeneous infinite solid, and the solution to the complementary elastic b.v.p., which corrects for the boundary conditions. The methodology is illustrated for the indentation of a single crystal by an array of flat rigid indenters in figure 1. The elastic dislocation fields are represented by a superscript (${}^{{\rm{d}}}$), the fields solving the complementary b.v.p. by a superscript ${(}^{\wedge })$. The elastic dislocation fields are given analytically. The dislocation in the repetitive cell and its periodic replicas are treated as an infinite array of dislocations. The solution to the b.v.p. is traditionally obtained using the finite-element method FEM. The rigid indenter is modeled implicitly by imposing boundary conditions on the deformable body.

Zoom In Zoom Out Reset image size

The total stress and displacement fields obtained at a given time increment are used, along with a set of constitutive rules, to describe the evolution of the dislocation structure. The constitutive rules are based on the Peach–Koehler force and control dislocation glide, nucleation, annihilation and pinning at obstacles. At the beginning of the simulation, the crystal is dislocation-free, but contains a density of point sources and obstacles that mimic Frank–Read sources and precipitates, respectively. The density of sources and obstacles is constant during the simulation. A dislocation pair is generated by a point source when the resolved shear stress acting on the source exceeds a critical nucleation strength, ${\tau }_{\mathrm{nuc}}$. The dislocation then glides on a plane with a velocity proportional to the Peach–Koehler force. Two dislocations with opposite Burgers vector annihilate when they come closer to each other than a threshold length, set to $6b$, where b is the magnitude of the Burgers vector. Whenever a dislocation meets an obstacle it gets pinned. It is released by the obstacle only when the Peach–Koehler force acting on the dislocation exceeds the critical strength of the obstacle, ${\tau }_{\mathrm{obs}}$. Dislocations can exit the domain through the free surface leaving behind a displacement step of $\pm b/2$ along the slip direction.

The aim of this work is to replace the FEM solution to the complementary b.v.p. with GFMD, while the constitutive rules that control the dislocation dynamics remain unchanged and similar to those proposed in [1].

GFMD is also a boundary-value method to study the elastic response of a body subjected to contact loading. In GFMD, only the surface of the deformable body is modeled explicitly and discretized using equi-spaced grid points as seen in figure 2. The grid points, under the influence of an external load, are displaced from their initial position, causing an increase in the areal elastic energy of the system. The new equilibrium positions are then calculated using damped dynamics in Fourier space. The advantage of damping the system in Fourier space is that the different modes describing the surface are uncoupled, and the stress–displacement relationship, when applying a load in z-direction, simply reads:

$$\begin{eqnarray}&&{\tilde{u}}_{\alpha }({\bf{q}})={\tilde{G}}_{\alpha \beta }({\bf{q}}){\tilde{\sigma }}_{z\beta }({\bf{q}}),\end{eqnarray} \tag{ 1 }$$

where ${G}_{\alpha \beta }({\bf{q}})$ is the Green's function tensor, ${u}_{\alpha }({\bf{q}})$ is the α component of the displacement and and ${\sigma }_{z\beta }({\bf{q}})$ is the traction in β-direction corresponding to wavenumber ${\bf{q}}$. The Green's function tensor depends on the elastic properties and size of the body.

Zoom In Zoom Out Reset image size

Here, the load is applied in an incremental manner by means of a rigid flat indenter, which is modeled by a hard-wall potential. To satisfy static equilibrium at each incremental change of the loading, the following condition must hold:

$$\begin{eqnarray}&&{\tilde{{\boldsymbol{F}}}}^{\mathrm{ext}}(q)+{\tilde{{\boldsymbol{F}}}}^{\mathrm{el}}(q)+{\tilde{{\boldsymbol{F}}}}^{\mathrm{if}}(q)=0,\end{eqnarray} \tag{ 2 }$$

where ${\tilde{{\boldsymbol{F}}}}^{\mathrm{ext}}(q)$ is the external force; ${\tilde{{\boldsymbol{F}}}}^{\mathrm{if}}(q)$ is the interfacial force ensuring the non-overlap constraint, which are imposed 'by hand' after each time step, in real space, ${\tilde{{\boldsymbol{F}}}}^{\mathrm{el}}(q)$ is the elastic restoring force that can be written as:

$$\begin{eqnarray}&&{\tilde{{\boldsymbol{F}}}}^{\mathrm{el}}(q)/{A}_{0}={\tilde{{\boldsymbol{G}}}}^{-1}(q)\tilde{{\boldsymbol{u}}}(q)={{\rm{\nabla }}}_{\tilde{u}}{v}_{\mathrm{el}},\end{eqnarray} \tag{ 3 }$$

where A₀ is the total surface area and ${\tilde{{\boldsymbol{G}}}}^{-1}(q)$ is the inverse Green's function, which can be evaluated from the areal elastic energy density ${v}_{\mathrm{el}}$. The areal elastic energy was derived for a slab with deformable top and fixed bottom in previous works [47, 48]. In section 4.1, we extend the derivation to the case where the bottom surface can be exposed to arbitrary stress, displacement, or mixed boundary condition. This is necessary for the coupling to dislocation dynamics, as should become clear in the following section.

The damping force has the form:

$$\begin{eqnarray}&&{\tilde{{\boldsymbol{F}}}}^{\mathrm{damp}}(q)=\eta ({\tilde{{\boldsymbol{u}}}}^{n-1}(q)-{\tilde{{\boldsymbol{u}}}}^{n}(q)),\end{eqnarray} \tag{ 4 }$$

where η is the damping factor, chosen such to critically damp the slowest mode, i.e., the center-of-mass mode, corresponding to q = 0, for quick convergence. Various further improvements can be applied to speed up convergence, such as, mode-dependent masses or damping.

The damping force is used in the position-Verlet algorithm to solve for the displacement fields at each increment $(n+1)$,

$$\begin{eqnarray}&&{\tilde{{\boldsymbol{u}}}}^{n+1}(q)=2{\tilde{{\boldsymbol{u}}}}^{n}(q)-{\tilde{{\boldsymbol{u}}}}^{n-1}(q)+({\tilde{{\boldsymbol{F}}}}^{\mathrm{el}}(q)+{\tilde{{\boldsymbol{F}}}}^{\mathrm{ext}}(q)+{\tilde{{\boldsymbol{F}}}}^{\mathrm{damp}}(q)){\tau }^{2},\end{eqnarray} \tag{ 5 }$$

where τ is the non-dimensional discrete time step used in the simulation.

The hard-wall potential is employed at the end of each iteration to ensure there is no inter-penetration, i.e., in real space,

$$\begin{eqnarray}&&{z}_{\mathrm{punch}}(x)\geqslant {z}_{\mathrm{substrate}}(x),\end{eqnarray} \tag{ 6 }$$

where ${z}_{\mathrm{punch}}$ and ${z}_{\mathrm{substrate}}$ are the z-coordinates of the punch and substrate surface, respectively as seen in the figure 2. Notice that the method is not bound to use a hard-wall potential. Finite interactions can be accounted for, however, we have here chosen for a hard-wall potential for the sake of comparison to conventional DDP.

When the surface equilibrates to the final deformed configuration, the body fields are calculated from the discrete surface fields using closed-form analytical solutions [48].

The GFDD method is based on the same decomposition concept as used in DDP except that, now, the image fields are found using GFMD instead of FEM. The methodology is schematically presented in figure 3 for the case of a layer indented using a flat rigid punch. Note that GFDD can solve generic boundary value problems where both top and bottom boundaries are arbitrarily partitioned into traction and displacement boundaries.

Zoom In Zoom Out Reset image size

When solving the complementary b.v.p., both tractions and displacements caused by the dislocations on the top and bottom boundary of the body need to be simultaneously prescribed at a given time increment. Tractions are imposed in Fourier space as described in the previous section by using $\hat{{\boldsymbol{t}}}(q)$ as external force ${\tilde{{\boldsymbol{F}}}}^{\mathrm{ext}}(q)$ before stepping forward in time. Here $\hat{{\boldsymbol{t}}}(q)$ is the Fourier transformation of the discontinuous function $\hat{{\boldsymbol{t}}}({\boldsymbol{r}})$ at point ${\boldsymbol{r}}=(x,z)$:

$$\begin{eqnarray}\,\hat{{t}}({\boldsymbol{r}}) & = & {t}({\boldsymbol{r}})-{{t}}^{{\rm{d}}}({\boldsymbol{r}})\qquad \mathrm{if}\,{\boldsymbol{r}}\in {S}_{{t}},\\ \hat{{t}}({\boldsymbol{r}}) & = & 0\qquad \qquad \qquad \ \mathrm{if}\,{\boldsymbol{r}}\in {S}_{{u}},\end{eqnarray} \tag{ 7 }$$

where S_t and S_u are traction- and displacement-prescribed boundaries, respectively. Unlike tractions, displacements are imposed in real space by setting the equilibrium position of the hard-wall to the required position:

$$\begin{eqnarray}&&\hat{{\boldsymbol{u}}}({\boldsymbol{r}})={\boldsymbol{u}}({\boldsymbol{r}})-{{\boldsymbol{u}}}^{{\rm{d}}}({\boldsymbol{r}})\qquad \mathrm{if}\ \,{\boldsymbol{r}}\in {S}_{u}.\end{eqnarray} \tag{ 8 }$$

The hard-wall is not applied to the boundaries where traction is prescribed ${S}_{{t}}$.

Notice that the elastic energy required for the evaluation of the Green's function, which is needed in the calculation of the restoring elastic force, was derived in [48] but only for the case of an isotropic slab with an undulated top layer and a fixed bottom. This allows a b.v.p. to be solved for a mixed boundary condition at the top, however, with the restriction of the bottom displacement to be zero. Here, however, we need to impose mixed boundary conditions also at the bottom. To this end, the areal elastic energy is required for a solid with both top and bottom undulation as seen in figure 4. This will be dealt with in the subsequent section.

Zoom In Zoom Out Reset image size

4.1. Elastic energy of an elastic layer loaded at both surfaces

A linear elastic isotropic body in a slab geometry is considered. The equilibrium condition for the case of no body forces can be written as ${\partial }_{\alpha }{\sigma }_{\alpha \beta }({\boldsymbol{r}})=0$, where ${\sigma }_{\alpha \beta }({\boldsymbol{r}})$ is the stress at point (x, z) represented by vector ${\boldsymbol{r}}$ and ${\partial }_{\alpha }\equiv \partial /\partial {r}_{\alpha }$. This can be written as:

$$\begin{eqnarray}[{C}_{11}{\partial }_{x}^{2}+{C}_{44}{\partial }_{z}^{2}]{u}_{x}({\boldsymbol{r}})+({C}_{44}+{C}_{12}){\partial }_{x}{\partial }_{z}{u}_{z}({\boldsymbol{r}}) & = & 0,\\ [{C}_{11}{\partial }_{z}^{2}+{C}_{44}{\partial }_{x}^{2}]{u}_{z}({\boldsymbol{r}})+({C}_{44}+{C}_{12}){\partial }_{z}{\partial }_{x}{u}_{x}({\boldsymbol{r}}) & = & 0,\end{eqnarray} \tag{ 9 }$$

where ${C}_{{ij}}$ denotes the coefficients of the elastic tensor. The in-plane wavenumber q is a scalar for the two-dimensional body considered here.

It is shown in appendix A that for the system of differential equations (9) the solutions of the in-plane cosine transform of the lateral u_x displacement field couples to the in-plane sine transform of the normal u_z displacement and vice versa. Thus, we can write:

$$\begin{eqnarray}{u}_{x}^{{\rm{c}}}(x,z) & = & \cos ({qx}){\tilde{u}}_{x}^{{\rm{c}}}(q,z),\\ {u}_{z}^{{\rm{s}}}(x,z) & = & \sin ({qx}){\tilde{u}}_{z}^{{\rm{s}}}(q,z).\end{eqnarray} \tag{ 10 }$$

Solutions satisfying the following boundary conditions:

$$\begin{eqnarray}{\tilde{u}}_{x}^{{\rm{c}}}(q,0) & = & {u}_{x}^{\mathrm{bot}},\\ {\tilde{u}}_{z}^{{\rm{s}}}(q,0) & = & {u}_{z}^{\mathrm{bot}},\\ {\tilde{u}}_{x}^{{\rm{c}}}(q,{z}_{{\rm{m}}}) & = & {u}_{x}^{\mathrm{top}},\\ {\tilde{u}}_{z}^{{\rm{s}}}(q,{z}_{{\rm{m}}}) & = & {u}_{z}^{\mathrm{top}}\end{eqnarray} \tag{ 11 }$$

and equation (9) are then obtained to satisfy

$$\begin{eqnarray}&&\left[\begin{array}{c}{\tilde{u}}_{x}^{{\rm{c}}}(q,z)\\ {\tilde{u}}_{z}^{{\rm{s}}}(q,z)\end{array}\right]=\left[\begin{array}{cccc}{h}_{1}(q,z) & {h}_{2}(q,z) & -{h}_{3}(q,z) & -{h}_{4}(q,z)\\ {h}_{3}(q,z) & {h}_{4}(q,z) & {h}_{5}(q,z) & {h}_{6}(q,z)\end{array}\right]\left[\begin{array}{c}{A}_{1}\\ {A}_{2}\\ {A}_{3}\\ {A}_{4}\end{array}\right]\end{eqnarray} \tag{ 12 }$$

with

$$\begin{eqnarray}{h}_{1}(q,z) & = & (1-r)\cosh ({qz})+{rqz}\sinh ({qz}),\\ {h}_{2}(q,z) & = & \sinh ({qz})+{rqz}\cosh ({qz}),\\ {h}_{3}(q,z) & = & r({qz}\cosh ({qz})-\sinh ({qz})),\\ {h}_{4}(q,z) & = & {rqz}\sinh ({qz}),\\ {h}_{5}(q,z) & = & (1+r)\cosh ({qz})-{rqz}\sinh ({qz}),\\ {h}_{6}(q,z) & = & \sinh ({qz})-{rqz}\cosh ({qz}),\end{eqnarray} \tag{ 13 }$$

where $r=\tfrac{1-s}{1+s}$ and $s={C}_{44}/{C}_{11}$. A_i can be found by applying the boundary conditions in equation (11):

$$\begin{eqnarray}&&\left[\begin{array}{c}{A}_{1}\\ {A}_{2}\\ {A}_{3}\\ {A}_{4}\end{array}\right]=\displaystyle \frac{1}{f(q,{z}_{{\rm{m}}})}\left[\begin{array}{cccc}\displaystyle \frac{f(q,{z}_{{\rm{m}}})}{(1-r)} & 0 & 0 & 0\\ {k}_{1}(q,{z}_{{\rm{m}}}) & {k}_{2}(q,{z}_{{\rm{m}}}) & {k}_{3}(q,{z}_{{\rm{m}}}) & {k}_{4}(q,{z}_{{\rm{m}}})\\ 0 & \displaystyle \frac{f(q,{z}_{{\rm{m}}})}{(1+r)} & 0 & 0\\ {k}_{5}(q,{z}_{{\rm{m}}}) & {k}_{6}(q,{z}_{{\rm{m}}}) & {k}_{7}(q,{z}_{{\rm{m}}}) & {k}_{8}(q,{z}_{{\rm{m}}})\end{array}\right]\left[\begin{array}{c}{\tilde{u}}_{x}^{{\rm{c}}}(q,0)\\ {\tilde{u}}_{z}^{{\rm{s}}}(q,0)\\ {\tilde{u}}_{x}^{{\rm{c}}}(q,{z}_{{\rm{m}}})\\ {\tilde{u}}_{z}^{{\rm{s}}}(q,{z}_{{\rm{m}}})\end{array}\right]\end{eqnarray} \tag{ 14 }$$

with

$$\begin{eqnarray}{k}_{1}(q,{z}_{{\rm{m}}}) & = & {{rqz}}_{{\rm{m}}}-\sinh ({{qz}}_{{\rm{m}}})\cosh ({{qz}}_{{\rm{m}}}),\\ {k}_{2}(q,{z}_{{\rm{m}}}) & = & -\displaystyle \frac{r}{1+r}({\cosh }^{2}({{qz}}_{{\rm{m}}})-1+r{({{qz}}_{{\rm{m}}})}^{2}),\\ {k}_{3}(q,{z}_{{\rm{m}}}) & = & {h}_{6}(q,{z}_{{\rm{m}}}),\\ {k}_{4}(q,{z}_{{\rm{m}}}) & = & {h}_{4}(q,{z}_{{\rm{m}}}),\\ {k}_{5}(q,{z}_{{\rm{m}}}) & = & \displaystyle \frac{r}{1-r}({\cosh }^{2}({{qz}}_{{\rm{m}}})-1-r{({{qz}}_{{\rm{m}}})}^{2}),\\ {k}_{6}(q,{z}_{{\rm{m}}}) & = & -{{rqz}}_{{\rm{m}}}-\sinh ({{qz}}_{{\rm{m}}})\cosh ({{qz}}_{{\rm{m}}}),\\ {k}_{7}(q,{z}_{{\rm{m}}}) & = & -{h}_{4}(q,{z}_{{\rm{m}}}),\\ {k}_{8}(q,{z}_{{\rm{m}}}) & = & {h}_{2}(q,{z}_{{\rm{m}}}),\end{eqnarray} \tag{ 15 }$$

and

$$\begin{eqnarray}f(q,{z}_{{\rm{m}}}) & = & {\cosh }^{2}({{qz}}_{{\rm{m}}})-{({{rqz}}_{{\rm{m}}})}^{2}-1\\ & = & \displaystyle \frac{1}{2}\{\cosh (2{{qz}}_{{\rm{m}}})-2{({{rqz}}_{{\rm{m}}})}^{2}-1\}.\end{eqnarray} \tag{ 16 }$$

Similarly, the in-plane sine transform of u_x and cosine transform of u_z can be obtained from:

$$\begin{eqnarray}&&\left[\begin{array}{c}{\tilde{u}}_{x}^{{\rm{s}}}(q,z)\\ {\tilde{u}}_{z}^{{\rm{c}}}(q,z)\end{array}\right]=\left[\begin{array}{cccc}{h}_{1}(q,z) & {h}_{2}(q,z) & {h}_{3}(q,z) & {h}_{4}(q,z)\\ -{h}_{3}(q,z) & -{h}_{4}(q,z) & {h}_{5}(q,z) & {h}_{6}(q,z)\end{array}\right]\left[\begin{array}{c}{B}_{1}\\ {B}_{2}\\ {B}_{3}\\ {B}_{4}\end{array}\right]\end{eqnarray} \tag{ 17 }$$

with

$$\begin{eqnarray}&&\left[\begin{array}{c}{B}_{1}\\ {B}_{2}\\ {B}_{3}\\ {B}_{4}\end{array}\right]=\displaystyle \frac{1}{f(q,{z}_{{\rm{m}}})}\left[\begin{array}{cccc}\displaystyle \frac{f(q,{z}_{{\rm{m}}})}{(1-r)} & 0 & 0 & 0\\ {k}_{1}({{qz}}_{{\rm{m}}}) & -{k}_{2}({{qz}}_{{\rm{m}}}) & {k}_{3}({{qz}}_{{\rm{m}}}) & -{k}_{4}({{qz}}_{{\rm{m}}})\\ 0 & \displaystyle \frac{f(q,{z}_{{\rm{m}}})}{(1+r)} & 0 & 0\\ -{k}_{5}({{qz}}_{{\rm{m}}}) & {k}_{6}({{qz}}_{{\rm{m}}}) & -{k}_{7}({{qz}}_{{\rm{m}}}) & {k}_{8}({{qz}}_{{\rm{m}}})\end{array}\right]\left[\begin{array}{c}{\tilde{u}}_{x}^{{\rm{s}}}(q,0)\\ {\tilde{u}}_{z}^{{\rm{c}}}(q,0)\\ {\tilde{u}}_{x}^{{\rm{s}}}(q,{z}_{{\rm{m}}})\\ {\tilde{u}}_{z}^{{\rm{c}}}(q,{z}_{{\rm{m}}})\end{array}\right].\end{eqnarray} \tag{ 18 }$$

From equation (12), the strains are calculated as:

$$\begin{eqnarray}{\tilde{\epsilon }}_{{xx}}^{{\rm{s}}}(q,z) & = & -q{\tilde{u}}_{x}^{{\rm{c}}}(q,z),\\ {\tilde{\epsilon }}_{{zz}}^{{\rm{s}}}(q,z) & = & {\partial }_{z}{\tilde{u}}_{z}^{{\rm{s}}}(q,z),\\ {\tilde{\epsilon }}_{{xz}}^{{\rm{c}}}(q,z) & = & {\partial }_{z}{\tilde{u}}_{x}^{{\rm{c}}}(q,z)+q{\tilde{u}}_{z}^{{\rm{s}}}(q,z).\end{eqnarray} \tag{ 19 }$$

Stresses are then obtained as usual through Hooke's law:

$$\begin{eqnarray}&&{\sigma }_{i}={C}_{{ij}}{\epsilon }_{j}.\end{eqnarray} \tag{ 20 }$$

Gathering all contributions to the elastic energy leads to

$$\begin{eqnarray}{v}_{\mathrm{el}} & = & \displaystyle \frac{{C}_{11}q}{2}\,[s{\tilde{\epsilon }}_{{xz}}^{{\rm{c}}}(q,{z}_{{\rm{m}}}){\tilde{u}}_{x}^{{\rm{c}}}(q,{z}_{{\rm{m}}})+\{{\tilde{\epsilon }}_{{zz}}^{{\rm{s}}}(q,{z}_{{\rm{m}}})+(1-2s){\tilde{\epsilon }}_{{xx}}^{{\rm{s}}}(q,{z}_{{\rm{m}}})\}{\tilde{u}}_{z}^{{\rm{s}}}(q,{z}_{{\rm{m}}})\\ & & +s{\tilde{\epsilon }}_{{xz}}^{{\rm{c}}}(q,0){\tilde{u}}_{x}^{{\rm{c}}}(q,0)+\{{\tilde{\epsilon }}_{{zz}}^{{\rm{s}}}(q,0)+(1-2s){\tilde{\epsilon }}_{{xx}}^{{\rm{s}}}(q,0)\}{\tilde{u}}_{z}^{{\rm{s}}}(q,0)]\\ & = & \displaystyle \frac{{C}_{11}q}{2}\left[\begin{array}{cccc}{\tilde{u}}_{x}^{{\rm{c}}}(q,0) & {\tilde{u}}_{z}^{{\rm{s}}}(q,0) & {\tilde{u}}_{x}^{{\rm{c}}}(q,{z}_{{\rm{m}}}) & {\tilde{u}}_{z}^{{\rm{s}}}(q,{z}_{{\rm{m}}})\end{array}\right]\left[\begin{array}{c}{H}_{{ij}}(q,{z}_{{\rm{m}}})\end{array}\right]\left[\begin{array}{c}{\tilde{u}}_{x}^{{\rm{c}}}(q,0)\\ {\tilde{u}}_{z}^{{\rm{s}}}(q,0)\\ {\tilde{u}}_{x}^{{\rm{c}}}(q,{z}_{{\rm{m}}})\\ {\tilde{u}}_{z}^{{\rm{s}}}(q,{z}_{{\rm{m}}})\end{array}\right]\end{eqnarray} \tag{ 21 }$$

with

$$\begin{eqnarray}{H}_{11}(q,{z}_{{\rm{m}}}) & = & {H}_{33}(q,{z}_{{\rm{m}}})=(1-r)\displaystyle \frac{\sinh ({{qz}}_{{\rm{m}}})\cosh ({{qz}}_{{\rm{m}}})-{{rqz}}_{{\rm{m}}}}{f(q,{z}_{{\rm{m}}})}\\ -{H}_{12}({{qz}}_{{\rm{m}}}) & = & {H}_{34}({{qz}}_{{\rm{m}}})=\displaystyle \frac{1-r}{1+r}\displaystyle \frac{(1-r){\sinh }^{2}({{qz}}_{{\rm{m}}})-2{({{rqz}}_{{\rm{m}}})}^{2}}{f(q,{z}_{{\rm{m}}})},\\ {H}_{13}(q,{z}_{{\rm{m}}}) & = & (1-r)\displaystyle \frac{{{rqz}}_{{\rm{m}}}\cosh (q,{z}_{{\rm{m}}})-\sinh ({{qz}}_{{\rm{m}}})}{f(q,{z}_{{\rm{m}}})},\\ -{H}_{14}(q,{z}_{{\rm{m}}}) & = & {H}_{23}(q,{z}_{{\rm{m}}})=(1-r)\displaystyle \frac{{{rqz}}_{{\rm{m}}}\sinh ({{qz}}_{{\rm{m}}})}{f(q,{z}_{{\rm{m}}})},\\ {H}_{22}(q,{z}_{{\rm{m}}}) & = & {H}_{44}(q,{z}_{{\rm{m}}})=(1-r)\displaystyle \frac{\sinh ({{qz}}_{{\rm{m}}})\cosh ({{qz}}_{{\rm{m}}})+{{rqz}}_{{\rm{m}}}}{f(q,{z}_{{\rm{m}}})},\\ {H}_{24}(q,{z}_{{\rm{m}}}) & = & -(1-r)\displaystyle \frac{{{rqz}}_{{\rm{m}}}\cosh (q,{z}_{{\rm{m}}})+\sinh ({{qz}}_{{\rm{m}}})}{f({{qz}}_{{\rm{m}}})}.\end{eqnarray} \tag{ 22 }$$

The complete elastic energy density containing the complex Fourier transform of the displacement with wavenumber q reads

$$\begin{eqnarray}{v}_{\mathrm{el}}=\displaystyle \sum _{q}\displaystyle \frac{{C}_{11}q}{2}\left[\begin{array}{cccc}{\tilde{u}}_{x}^{* }(q,0) & {\tilde{u}}_{z}^{* }(q,0) & {\tilde{u}}_{x}^{* }(q,{z}_{{\rm{m}}}) & {\tilde{u}}_{z}^{* }(q,{z}_{{\rm{m}}})\end{array}\right]\left[\begin{array}{c}{H}_{{ij}}(q,{z}_{{\rm{m}}})\end{array}\right]\left[\begin{array}{c}{\tilde{u}}_{x}(q,0)\\ {\tilde{u}}_{z}(q,0)\\ {\tilde{u}}_{x}(q,{z}_{{\rm{m}}})\\ {\tilde{u}}_{z}(q,{z}_{{\rm{m}}})\end{array}\right]\end{eqnarray} \tag{ 23 }$$

with

$$\begin{eqnarray}\left[\begin{array}{c}{H}_{{ij}}(q,{z}_{{\rm{m}}})\end{array}\right]=\left[\begin{array}{cccc}{H}_{11}(q,{z}_{{\rm{m}}}) & -{{\rm{i}}{H}}_{12}(q,{z}_{{\rm{m}}}) & -{{\rm{i}}{H}}_{13}(q,{z}_{{\rm{m}}}) & -{{\rm{i}}{H}}_{14}(q,{z}_{{\rm{m}}})\\ {{\rm{i}}{H}}_{21}(q,{z}_{{\rm{m}}}) & {H}_{22}(q,{z}_{{\rm{m}}}) & -{{\rm{i}}{H}}_{23}(q,{z}_{{\rm{m}}}) & -{{\rm{i}}{H}}_{24}(q,{z}_{{\rm{m}}})\\ {{\rm{i}}{H}}_{31}(q,{z}_{{\rm{m}}}) & {{\rm{i}}{H}}_{32}(q,{z}_{{\rm{m}}}) & {H}_{33}(q,{z}_{{\rm{m}}}) & -{{\rm{i}}{H}}_{34}(q,{z}_{{\rm{m}}})\\ {{\rm{i}}{H}}_{41}(q,{z}_{{\rm{m}}}) & {{\rm{i}}{H}}_{42}(q,{z}_{{\rm{m}}}) & {{\rm{i}}{H}}_{43}(q,{z}_{{\rm{m}}}) & {H}_{44}(q,{z}_{{\rm{m}}})\end{array}\right].\end{eqnarray} \tag{ 24 }$$

The body fields are obtained through the closed-form analytical expressions in equations (12) and (17). The displacement fields hence obtained are compared with those obtained by FEM in figures 5 and 6. The relative difference in the displacement field obtained using both methods are found to be below 0.25%. It has to be noted that periodic boundary conditions in GFMD are intrinsically enforced through the periodicity of the Fourier transforms. In FEM, periodicity can be imposed using various methods, including the penalty method, as done in this study, or Lagrangian multipliers, but is never exact.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The elastic energy density in equation (23) is extended to the case of a semi-infinite half-space in appendix B. This opens up the possibility of modeling the plastic contact response of a semi-infinite body using dislocation-dynamics simulations. This is beneficial to study the plastic response of a body under contact loading without the effect of its bottom, i.e., an increase in contact pressure caused by dislocations piling up at the bottom of the body.

In this section, the new GFDD model is compared to DDP when computing the image fields for the simplest case scenario: a single dislocation pinned in an isotropic slab with Young's modulus E = 70 GPa and Poissons ratio $\nu =0.33$. The magnitude of the Burger's vector is $b=0.25\,\mathrm{nm}$.

Simulations are carried out for a unit cell with the bottom fixed and a traction free top surface, i.e.,

$$\begin{eqnarray}{u}_{x}(x,0) & = & {u}_{z}(x,0)=0,\\ {\sigma }_{{xz}}(x,{z}_{{\rm{m}}}) & = & {\sigma }_{{zz}}(x,{z}_{{\rm{m}}})=0.\end{eqnarray} \tag{ 25 }$$

The stress distribution obtained using GFDD is compared with DDP in figure 7. The displacements at the top surface, where tractions are zero, and the tractions at the bottom surface, where displacement are zero, are shown in figures 8 and 9.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It is found that the tractions at the surface obtained using GFDD suffer from ringing, also known as Gibb's phenomenon (see figure 10). This is because the discontinuities in the displacement imposed at the surfaces cause the higher harmonics of the traction to have higher amplitudes than the lower harmonics. To remove the ringing artifacts the results displayed in figure 9 are obtained after multiplying the traction ${\boldsymbol{t}}(q)$ with a sinc function $\mathrm{sinc}({{qa}}_{0})$, where a₀ is the discretization length. This is equivalent to convolving in real space the traction suffering from ringing with a rectangular box of unit height and width equal to the discretization length. Notice, however, that ringing affects only surface tractions, not surface displacements, from which the body fields are calculated.

Zoom In Zoom Out Reset image size

This benchmark problem is used to compare GFDD to classical DDP. The simulations are carried out for a unit cell that is indented by a rigid punch as in figure 11.

Zoom In Zoom Out Reset image size

6.1. Boundary-value problem

The indentation is prescribed by specifying the normal displacement rate along the contact of length ${L}_{x}^{{\rm{p}}}$:

$$\begin{eqnarray*}&&\dot{{u}_{z}}(x,{z}_{{\rm{m}}})=-{\dot{u}}_{z}^{0},\quad x\in \left[\displaystyle \frac{{L}_{x}-{L}_{x}^{{\rm{p}}}}{2},\displaystyle \frac{{L}_{x}+{L}_{x}^{{\rm{p}}}}{2}\right].\end{eqnarray*}$$

A sticking contact is modeled in DDP by taking the lateral displacement u_x = 0 in the contact region. In GFDD, the lateral movement in the contact region is constrained horizontally through the hard-wall potential. The non-contact part of the top surface of the unit cell is taken to be traction-free,

$$\begin{eqnarray*}&&{\sigma }_{{xz}}(x,{z}_{{\rm{m}}})={\sigma }_{{zz}}(x,{z}_{{\rm{m}}})=0.\end{eqnarray*}$$

This is achieved in GFDD by letting the contact points relax to equilibrium without any constraints. Finally, the bottom of the unit cell, z = 0 is fixed:

$$\begin{eqnarray*}&&{u}_{x}(x,0)={u}_{z}(x,0)=0.\end{eqnarray*}$$

In GFDD, this is implemented by constraining the lateral and normal motion of grid points at the bottom surface using the hard-wall potential.

6.2. Choice of parameters

Calculations are carried out for crystals with aspect ratio $a={z}_{{\rm{m}}}/{L}_{x}=1/2$, contact fractions ${L}_{x}^{{\rm{p}}}/{L}_{x}=1/12$ and ${L}_{x}=12\,\mu {\rm{m}}$. The elastic constants are chosen to represent aluminum: the Young's modulus is E = 70 GPa and Poissons ratio $\nu =0.33$.

For the DDP simulations, the slab is discretized using a uniform mesh of square elements. The number of degrees of freedom is ${n}_{\mathrm{dof}}=2{nnx}\times {nnz}$, where nnx is the number of nodes in x-direction, and nnz the number of nodes in z-direction. For the GFDD simulations, the surface is discretized using nx equi-spaced grid points, with nx = nnx.

In GFDD, the center-of-mass mode is critically damped or slightly under-damped for quick convergence. The damping factor η is

$$\begin{eqnarray}&&\eta \propto \displaystyle \frac{1}{\tau \sqrt{{nx}}},\end{eqnarray} \tag{ 26 }$$

where τ is the non-dimensional time step shown earlier in equation (5) that is used in the position-Verlet algorithm in GFMD to solve for the unknown displacement fields and it is taken to be τ = 0.25. The number of iterations used to reach convergence scales as ${n}_{\mathrm{it}}\propto \sqrt{a\,{nx}}$.

The dislocations can glide on three sets of parallel slip planes, with slip plane orientations: $0^\circ ,\,60^\circ$ and 120^◦ to the top surface. The discrete slip planes are spaced at $200b$ where $b=0.25\,\mathrm{nm}$ is the length of the Burger's vector. The Frank–Read sources and obstacles are randomly distributed in the crystal with a density ${\rho }_{\mathrm{nuc}}=30\,\mu {{\rm{m}}}^{-2}$ and ${\rho }_{\mathrm{obs}}=30\,\mu {{\rm{m}}}^{-2}$ in an initially dislocation free crystal. The strength of the sources follows a Gaussian distribution with mean strength ${\tau }_{\mathrm{nuc}}=50\,\mathrm{MPa}$ and standard deviation of $10\,\mathrm{MPa}$. The critical time for nucleation is ${t}_{\mathrm{nuc}}=0.1\,\mathrm{ns}$. The strength of the obstacles is taken to be $150\,\mathrm{MPa}$. Dislocations of opposite sign in the same slip plane annihilate when the distance between them is below ${L}_{\mathrm{ann}}=6b$. The time step in both DDP and GFDD simulations is taken to be ${\rm{\Delta }}t=2.5$ ns.

6.3. A simple dislocation dynamic simulation: a single Frank–Read source

In this section, a simple problem is considered where a rigid flat punch indents a crystal containing a single Frank–Read source as shown in figure 12(a). In order to observe appreciable plastic deformation in the material, the magnitude of the Burger's vector is magnified four times, i.e., $b=1.0\,\mathrm{nm}$. The mean contact pressure obtained using both methods is displayed in figure 12(b). While the flat rigid punch indents the layer, the source keeps generating dislocation dipoles, causing periodic kinks in the pressure–displacement curve. The difference in mean contact pressure in figure 12(b) is not seen to the naked eye.

Zoom In Zoom Out Reset image size

The surface fields obtained using both methods at ${u}_{z}=0.01\mu {\rm{m}}$ are shown in figure 13. The displacement steps formed due to the exiting of dislocations can be clearly seen close to $x/{L}_{x}=0.6$.

Zoom In Zoom Out Reset image size

Note that in figure 13(a) the curves for the u_z displacement overlap since in z-direction displacement boundary conditions are imposed. The difference between the curves representing u_x stems from the numerical difference in calculating the resolved shear stress acting on the source and the location of the dislocations in the two numerical schemes. The calculated stress field depends not only on boundary conditions when the simulation is elastic, but also on the location of other dislocations in the crystal when there is plasticity. Therefore, the differences builds up with increasing dislocation density. However, relative differences between contact pressure claculated in figure 13(b) with two methods remain below $0.45 \%$.

6.4. Dislocation dynamics simulation with many sources and obstacles

In this section the indented crystals contain a density of Frank–Read sources ${\rho }_{\mathrm{nuc}}=30\,\mu {{\rm{m}}}^{-2}$ and a density of obstacles ${\rho }_{\mathrm{obs}}=30\,\mu {{\rm{m}}}^{-2}$. The simulations are carried out with DDP and GFDD on crystals containing the same realization of sources and obstacles, i.e., the location as well as the strength of sources and obstacles are identical. Figure 15 shows the stress state and dislocation distribution at final indentation depth. Note that there is no one-to-one correspondence between the dislocations and therefore also not in terms of stress distribution. This is not surprising given that a tiny difference in the evolution of the dislocation structure, like a small delay in the nucleation of a dislocation or in the formation of a junction, would trigger an avalanche of differences in the following dislocation dynamics [49]. The overall features such as the shear bands emitted by the contact are captured by both methods in the same way. This is also testified by the mean contact pressure in figure 14(a) for the simulation presented in figure 15 and for a different realization. While DDP and GFDD do not produce identical mean pressures as a function of displacement for a given realization of Frank–Read sources, differences tend to be larger within one method from one realization to the next. In figure 14(b) are presented the average between the three realizations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

6.5. Simulation time

The computational complexity involved in solving the elastic b.v.p. using GFDD is only $O({nx}\sqrt{{nx}}\,\mathrm{log}{nx})$ [48], while it is $O({{nx}}^{2}{B}^{2})$ [50] in DDP, where B is the mean bandwidth of the stiffness matrix, which cannot exceed nx.

The time consuming part in 2D dislocation dynamics is the calculation of the resolved shear stress ${\tau }_{\mathrm{res}}$ at the location of objects, i.e., sources, dislocations and obstacles. In DDP, this requires searching the element where the object is located and subsequently calculating the stress and interpolating it to the location of the object. This procedure scales as $O({{nx}}^{2})$. In GFDD, instead, the resolved shear stress can be evaluated directly at the points of interest, i.e., at dislocations, sources and obstacles, which requires a smaller computational effort, scaling with O(nx). This is because the body field is calculated based on surface displacements $\tilde{{\boldsymbol{u}}}(q)$ using ${nx}/2$ modes.

The simulation time required for elasticity and for the calculation of the resolved shear stress is displayed in figure 16, and shows in both cases how the computational advantage of using GFDD increases with increasing discretization. It has to be noted that we use a skyline solver for FEM. The computational time using this solver was found to be of the same order as that of iterative solvers. The simulations performed on a single Intel Xeon(R) 3.10 GHz processor with 31.3 Gbytes of RAM.

Zoom In Zoom Out Reset image size

The dislocation dynamics is computed using the same algorithm in both methods and takes therefore the same amount of time and resources. The time required for the dynamics is independent of the discretization and increases with dislocation density. For the DDP simulations performed here, the time consumed by the dislocation dynamics is a negligible fraction of the time required to compute the resolved shear stress. GFDD is thus computationally more efficient than DDP independently of dislocation density.

It has to be noted that the maximum number of surface nodes chosen for this study is only 2⁹. If one intends to study the contact response of a realistic self-affine surface where roughness scales over three orders of magnitude ranging in scale from 50 nm to 100 μm [3], the surface has to be discretized by at least 2¹³ points (or even more depending on the thermodynamic and continuum limit [2]). For such large systems, and small strain simulations, as typically used in dislocation dynamics, the computational complexity for elastic b.v.p. was found to always dominate in comparison to the dislocation dynamics. In this case the computational advantage of GFDD becomes even more appreciable.

Notice also that the benchmark problem chosen in this work involves a constant contact area. For contact problems where the area is not constant, DDP becomes even slower since finding the correct contact area by means of the FEM requires many iterations as well as updating the boundary conditions at each time increment. GFDD is inherently impervious to such issues since it employs an interaction potential between the contacting surfaces.

In this work, we propose a modeling technique, GFDD, which combines GFMD with DDP. We firstly extended the existing GFMD model such that it can simulate an elastic layer with arbitrary loading at both the top and bottom surfaces. To this end we derived the areal elastic energy for the case of an isotropic layer with sinusoidal loading at both ends. In addition, we derived the body fields required to capture the evolution of the dislocation structure. The results obtained using GFDD are compared with conventional DDP for a benchmark problem: periodic indentation of a single crystal by flat punches.

The mean contact pressure during indentation using the two methods is found to differ less than two different realizations using the same method. Here by realization is intended a given initial distribution of dislocation sources and obstacles. The differences between the two methods stems from the evaluation of the fields using different discretizations: GFDD discretizes only the surface, DDP also the body.

The new GFDD model has various advantages compared to classical DDP. First, it is faster and opens up the possibility of studying realistic rough surfaces by exploiting a larger number of degrees of freedom. Next, GFDD employs an interaction potential between the contacting bodies, and does not involve time-consuming algorithms to keep track of the evolution of the contact area. Also, the periodicity in GFDD is intrinsically enforced through Fourier transforms, making it a better candidate than DDP to study contact problems by exploiting the periodicity of the unit cell on which the analysis is performed. Obviously, this is also a limitation of the GFDD model, which is currently not suitable to study non-periodic problems. Extension of the model to overcome this limitation seems an interesting avenue for future research. Additionally, the GFDD model has the potential to serve as a platform for multi-scale modeling where the surface has an explicit atomistic description and the bulk can be treated as a dislocated continuum.

This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement no. 681813). LN also acknowledges support by the Dutch National Scientific Foundation NWO and Dutch Technology Foundation STW (VIDI grant 12669).

It was shown in [48] that assuming an in-plane undulation of the top surface of an isotropic layer with a real-valued wavenumber q, equation (9) can be solved with the factorization

$$\begin{eqnarray}&&{u}_{\alpha }(x,z)={u}_{\alpha }^{0}\exp ({\rm{i}}{q}{x})\exp (\pm {qz}).\end{eqnarray} \tag{ A.1 }$$

Substituting equation (A.1) in the equilibrium condition (9), it can be rewritten in Fourier space as

$$\begin{eqnarray}{C}_{11}{\tilde{u}}_{x}(q,z)-{C}_{44}{\tilde{u}}_{x}(q,z)-{\rm{i}}({C}_{12}+{C}_{44}){\tilde{u}}_{z}(q,z) & = & 0\\ {C}_{11}{\tilde{u}}_{z}(q,z)-{C}_{44}{\tilde{u}}_{z}(q,z)-{\rm{i}}({C}_{12}+{C}_{44}){\tilde{u}}_{x}(q,z) & = & 0.\end{eqnarray} \tag{ A.2 }$$

Equation (A.2) then reduces to

$$\begin{eqnarray}&&{\tilde{u}}_{x}(q,z)={\rm{i}}{\tilde{u}}_{z}(q,z).\end{eqnarray} \tag{ A.3 }$$

This implies that

$$\begin{eqnarray}&&{\tilde{u}}_{x}^{c}(q,z)-{\rm{i}}{\tilde{u}}_{x}^{s}(q,z)={\rm{i}}{\tilde{u}}_{z}^{c}(q,z)+{\tilde{u}}_{z}^{s}(q,z).\end{eqnarray} \tag{ A.4 }$$

This explains how the solutions of the in-plane cosine transform of the lateral u_x displacement field couples to the in-plane sine transform of the normal u_z displacement, and vice versa.

For the limiting case in which the height of the slab tends to infinity, $H(q,{z}_{{\rm{m}}})$ can be written as:

$$\begin{eqnarray}{H}_{11}({{qz}}_{{\rm{m}}}\gg 1) & = & {H}_{33}({{qz}}_{{\rm{m}}}\gg 1)=\displaystyle \frac{2s}{1+s},\\ -{H}_{12}({{qz}}_{{\rm{m}}}\gg 1) & = & {H}_{34}({{qz}}_{{\rm{m}}}\gg 1)=\displaystyle \frac{2{s}^{2}}{1+s},\\ {H}_{13}({{qz}}_{{\rm{m}}}\gg 1) & = & 0,\\ -{H}_{14}({{qz}}_{{\rm{m}}}\gg 1) & = & {H}_{23}({{qz}}_{{\rm{m}}}\gg 1)=0,\\ {H}_{22}({{qz}}_{{\rm{m}}}\gg 1) & = & {H}_{44}({{qz}}_{{\rm{m}}}\gg 1)=\displaystyle \frac{2s}{1+s},\\ {H}_{24}({{qz}}_{{\rm{m}}}\gg 1) & = & 0.\end{eqnarray} \tag{ B.1 }$$

In the limiting case of short wave-numbers, we find

$$\begin{eqnarray}{{qH}}_{11}({{qz}}_{{\rm{m}}}\ll 1) & = & {{qH}}_{33}({{qz}}_{{\rm{m}}}\ll 1)=\displaystyle \frac{{C}_{44}}{{z}_{{\rm{m}}}},\\ -{{qH}}_{12}({{qz}}_{{\rm{m}}}\ll 1) & = & {{qH}}_{34}({{qz}}_{{\rm{m}}}\ll 1)=0,\\ {{qH}}_{13}({{qz}}_{{\rm{m}}}\ll 1) & = & -\displaystyle \frac{{C}_{44}}{{z}_{{\rm{m}}}},\\ -{{qH}}_{14}({{qz}}_{{\rm{m}}}\ll 1) & = & {{qH}}_{23}({{qz}}_{{\rm{m}}}\ll 1)=0,\\ {{qH}}_{22}({{qz}}_{{\rm{m}}}\ll 1) & = & {{qH}}_{44}({{qz}}_{{\rm{m}}}\ll 1)=\displaystyle \frac{{C}_{11}}{{z}_{{\rm{m}}}},\\ {{qH}}_{24}({{qz}}_{{\rm{m}}}\ll 1) & = & -\displaystyle \frac{{C}_{11}}{{z}_{{\rm{m}}}}.\end{eqnarray} \tag{ B.2 }$$

In the case of bottom undulation going to zero ${\boldsymbol{u}}(q,0)\to 0$, we recover the elastic energy for a fixed bottom derived in [48].