Constitutive relations for plasticity of amorphous carbon

We deform representative volume elements of amorphous carbon obtained from melt-quenches in molecular dynamics calculations using bond-order and machine learning interatomic potentials. A Drucker-Prager law with a zero-pressure flow stress of $41.2$~GPa and an internal friction coefficient of $0.39$ describes the deviatoric stress during flow as a function of pressure. We identify the mean coordination number as the order parameter describing this flow surface. However, a description of the dynamical relaxation of the quenched samples towards steady-state flow requires an additional order parameter. We suggest an intrinsic strain of the samples as a possible order parameter and present equations for its evolution. Our results provide insights into rehybridization and pressure dependence of friction between coated surfaces as well as routes towards the description of amorphous carbon in macroscale models of deformation.

a-C is interesting not just for its wide range of applications but also because it forms an ideal network structure (Fig. 1a). Carbon atoms can be sp-(two neighbors), sp 2 -(three neighbors) or sp 3 -(four neighbors) hybridized. The pair-distribution function, shown in Fig. 1b, vanishes between the first and second neighbor peak, in contrast to metallic glasses 15 or even amorphous silicon. 16 This means a-C forms an ideal network; it is the only singlecomponent network-forming glass. Molecular dynamics has in the past been used to compute yield of polymer and network glasses. For example, Rottler & Robbins 17 showed that yield of polymer glasses described by bead-spring models follows a Drucker-Prager 18 or pressure-modified von-Mises law. Their model glasses yielded once the deviatoric (von-Mises) stress τ dev exceeded where p is the hydrostatic pressure. Similar behavior was found by Molnár et al. 19 for silicate glasses modeled with the BKS potential. 20,21 Since Eq. (1) looks like Amontons' friction law with an adhesive contribution, α is often called the internal friction coefficient.
Experimentally, Drucker-Prager-type behavior has been found for polymers, 22,23 foams 24 and metallic glasses. [25][26][27] Along the line of these prior simulation works, we deform representations of the network glass a-C in direct non-equilibrium molecular dynamics calculations. Specifically, we use simple shear (up to 100% strain, see Fig. 1c) and uniaxial compression (up to 50% see Fig. 1d) at an applied strain rate ofε = 10 9 s −1 to map out a representative portion of the yield surface. Our molecular dynamics calculations start from models of amorphous carbon consisting of ∼ 4000 atoms. These models are obtained by randomly placing the atoms inside a box of fixed volume, yielding model systems in a range of controlled densities ρ from 2.0 g cm −3 to 3.5 g cm −3 . All subsequent calculations are carried out at this fixed volume. We equilibrate these systems for 25 ps at 5000 K after which we quenched the system to 300 K with time constant 0.5 ps. The details of the quench protocol do not appear to matter as the system loses memory of its initial state during plastic deformation.
The quench protocol also does not affect structure and elastic properties of the samples, except for very slow quenches where the system may crystallize. [28][29][30] We use two interatomic force models that follow competing philosophies: The screened variant of the Tersoff III potential 31,32 (in the following denoted by Tersoff+S) and the Gaussian approximation potential 33 (denoted by GAP) as recently parameterized for a-C. 34 The former potential was designed to correctly describe bond-breaking processes 35 as those continuously occurring during plastic deformation; the latter machine-learning potential gives an accuracy comparable to density-functional theory within the local-density approximation 36 that was used to train it. Both potentials predict structure and mechanical properties within similar uncertainties of experimental measures 28-30 (e.g. see Fig. 1b). Temperature is controlled to 300 K using a Langevin thermostat with a relaxation time constant of 0.5 ps.
In the case of the simple shear deformation the thermostat was only applied in the direction perpendicular to the shear plane.
Both interatomic potentials yield a glassy disordered carbon network at the quench rates employed here, and their elastic properties are isotropic. 30 The isotropic nature of a-C implies that the yield surface can only depend on the principal stresses. From the principal stresses σ 1 , σ 2 and σ 3 in our simulations, we can calculate the first two invariants of the stress tensor, the hydrostatic pressure and deviatoric (von-Mises) stress Figure 1e shows τ dev as a function of the applied strain ε for three select cases: An initially linear (pseudo-) elastic response is followed by yield and then flow of the material at almost constant stress. The denser samples show shear-softening and we do not find an appreciable difference in the stress-strain response between simple shear and uniaxial compression.
Our simulations are carried out at constant volume and the hydrostatic pressure changes with applied strain. Figure 1f shows the deviatoric shear stress τ dev as a function of hydrostatic pressure p throughout our simulations. The pressure is constant at small applied strain where the material responds elastically. The nonzero pressure is a residue of the quenching process; we quench at constant volume and do not relax the simulation cell after the quench. The volume elements are under tensile (low density) or compressive (high density) stress. The hydrostatic pressure increases in all cases but then saturates as the material flows. This pressure increases because a-C expands in volume when plastically deformed.
Volume expansion has been previously reported in studies of wear of a-C 5 and diamond. 37,38 The reason for this expansion in volume is that shearing equilibrates the a-C's structure towards the structure of the liquid phase. 38  Figure 2b shows τ dev as a function of n, computed by counting neighbors within a cutoff of r c = 1.85Å where the pair distribution has dropped to zero (Fig. 1b). We find a linear dependence for both potentials, but with different slopes and different intercepts. Extrapolating τ dev (n) to τ dev = 0 we find that the GAP-glass loses rigidity at n = 2.4, the mean field prediction, while the Tersoff+S-glass loses rigidity at a higher mean coordination of n ≈ 2.8.
We believe that the difference between the two model glasses relates back to the idea of rigidity percolation. The limit n = 2.4 only holds for a continuous random network. For general networks, rings with more than 6 members are floppy and can form floppy regions within the material. 39 In Fig. 2c ring statistics 42 are shown for two systems at ρ ≈ 2.75 g/cm 3 .
The Tersoff+S structure contains notably more rings with sizes between 8 and 11 and those rings are floppy. Figure 2d shows the fraction of rings with sizes larger than six as a function of the mean coordination number n in the structures. At coordination numbers of 3.8 and above both potentials agree very well, but below Tersoff+S contains a much higher fraction of large, floppy rings towards the coordination where the whole system becomes floppy.
In the inelastic regime, our simulations show a drop of the mean coordination number n with applied strain (Fig. 3a): The material rehybridizes. 5,43,44 Atoms with lower coordination require more volume and hence the pressure during our constant-volume simulations rises.
This pressure is partially due to elastic deformation. The relaxed a-C systems follow a unique relationship between density and coordination number, ρ 0 (n). Figure 3b shows this relationship as obtained from the well-equilibrated simulations reported in Ref. 30. Similarly, the bulk modulus is shown in Fig. 3c to uniquely depend on density, B 0 (ρ 0 ). The pressure inside our simulation cell must therefore be given by with total volumetric strain We call ε int is the intrinsic (or residual) strain. With applied strain, the mean coordination number n decreases and ε int increases. The GAP trajectories also start at ε int ≈ 0, but n and ε int show less variation with applied strain than the Tersoff+S trajectories. The average values for ε int are lower than for Tersoff+S.
From the total volumetric strain ε V (open symbols in Fig. 3d), we see that the intrinsic strain is the dominant contribution to the overall volumetric strain in the system. The evolution of the hydrostatic pressure in our simulations can therefore be related to the evolution of the intrinsic strain during deformation. Our interpretation of the intrinsic strain is that deformation leads to a distortion of the atomic structure that changes its volume.
This distortion may be difficult to quantify in geometric terms, similar to the difficulty of The coordination number n alone is therefore not a sufficient order parameter for the description of the state of the material. A constitutive model for a-C requires the introduction of an additional state variable, for example the intrinsic strain ε int directly. As shown in Fig. 3d, the combined macroscopic state vector (n, ε int ) evolves towards a manifold of steady-state values that is shown by the dashed line in Fig. 3d and can be described by a functional relationship ε int,0 (n). The relaxation towards this steady-state behavior is for example shown in Fig. 1e and 3a. Assuming rate-independence with a characteristic relaxation strain ε c , an approximate evolution law for the state vector in the spirit of a relaxation time approximation is The target of the relaxation, the steady-state coordination number n 0 (ε int , n) depends on the current state, as can be directly seen in Fig. 3d. We can extract the steady-state behavior by following along the pathway of deformation in Fig. 3d. Given δn(ε int ) = dn/dε int as the slope of the evolution of n(ε int ) in this figure, we find n 0 as the solution of the nonlinear equation for each state (n, ε int ). Ingredients to this constitutive law are the tangent δn(ε int ), the steady-state intrinsic strain ε int,0 (n 0 ) and the relaxation constant ε c . Note that the relaxation constants for ε int and n in Eq. (6) could differ and would need to be determined from additional calculations not presented here.
The solid lines in Fig. 3d show a solution of this model for ε c = 0.1 within the order parameter space. As shown by the dashed lines in Fig. 3a, this solution describes the evolution of the coordination number n with applied strain ε well. It also serves as a partial explanation for the shear-softening behavior seen at high densities. The deviatoric stress τ dev drops (see Fig. 1a) because the coordination number decreases and this weakens the material. Using the linear τ dev (n) dependency shown in Fig. 2b, we obtain the dashed lines in Fig. 1e that qualitatively capture the response of the material.
Note that the set of equation presented here cannot describe the response to a change in the density that occurs along the dashed line in Fig. 3d. To describe this behavior, Eq. 7 must couple to the density ρ or the total pressure p, and additional calculations are required to extract an approximate mathematical description of this coupling required for a fully-formulated constitutive law.
Finally, we note that there are large differences between the behavior of the Tersoff+S and the GAP glass: Tersoff+S has a stronger tendency towards rehybridization. This means that for GAP, we cannot extract δn(ε int ) as n shows variation only by 0.05 and the resolution with which we can resolve changes in n depends on the total number of atoms in our unit cell. We expect that a similar picture emerges for GAP but are at present limited to ∼ 4000 atoms because of the computational cost of the GAP potential. Despite these differences in the structural changes of the material, the yield surface (Fig. 2a) appears independent of the choice of interatomic potential.
In summary, we find that steady-state flow of a-C is described by a Drucker-Prager law. Model glasses obtained from two different interatomic potentials collapse onto the same Drucker-Prager law, giving confidence to the extracted parameters. Our model glasses behave differently with regards to the evolution of the mean coordination number of the system (or alternatively, the numbers of sp 3 -, sp 2 -and sp-hybridized atoms). We can extract a constitutive relationship for these models that involves an intrinsic strain of these structures as an additional order parameter. These results are the first parameterization of the yield surface of a-C. They have relevance for understanding the rehybridization and friction of a-C surfaces in sliding contact that has been observed experimentally 2,43 and in simulations. 5,44 Our results also open a route for the development of constitutive models for macroscale calculations of plastic deformation or fracture in a-C.