<figure><img src="image_1.png" /></figure>

HAL
open science

# Densification dependent yield criteria for sodium silicate glasses - An atomistic simulation approach

Gergely Molnar, Patrick Ganster, Anne Tanguy, Etienne Barthel, Guillaume Kermouche

<figure><img src="image_2.png" /></figure>

## To cite this version:

Gergely Molnar, Patrick Ganster, Anne Tanguy, Etienne Barthel, Guillaume Kermouche. Densification dependent yield criteria for sodium silicate glasses - An atomistic simulation approach. Acta Materialia, 2016, 111, pp.129-137. 10.1016/j.actamat.2016.03.053. emse-01353610

HAL Id: emse-01353610

https://hal-emse.ccsd.cnrs.fr/emse-01353610

Submitted on 8 Sep 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Densification dependent yield criteria for sodium silicate
glasses – an atomistic simulation approach

Gergely Molnár<sup>a,*</sup>, Patrick Ganster<sup>a</sup>, Anne Tanguy<sup>b</sup>, Etienne Barthel<sup>c,d</sup>,
Guillaume Kermouche<sup>a</sup>

<sup>a</sup>École de Mines de Saint-Etienne, Centre SMS, Laboratoire Georges Friedel
CNRS-UMR5307, 158 Cours Fauriel, 42023, Saint-Etienne, France

<sup>b</sup>Laboratoire de Mécanique des Contacts et des Structures, Institut National des Sciences
Appliquées de Lyon 18-20, rue des Sciences, 69621, Villeurbanne Cedex, France

<sup>c</sup>École Supérieure de Physique et de Chimie Industrielles de la Ville de Paris (ESPCI)
ParisTech, PSL Research University, Sciences et Ingénierie de la matière Molle, CNRS
UMR 7615, 10, Rue Vauquelin, F-75231 Paris Cedex 05, France

<sup>d</sup>Sorbonne-Universités, UPMC Univ. Paris 06, SIMM, 10, Rue Vauquelin, F-75231
Paris Cedex 05, France

Abstract

Silicate glasses are macroscopically brittle but ductile at the micron scale. This plastic response is complex: in open structure materials, such as amorphous silica, plastic yield results in significant densification. While, more compact structures (e.g. soda-silicate glasses) are known to suppress densification and promote shear flow. We have carried out atomic scale simulations to analyze the plastic response of a series of silicates with increasing sodium content. Quasi-static, multi-axial deformation tests were performed on large samples (≈ 10³ nm³). Their yield behavior was quantified at different stress states, by measuring permanent volume changes. Qualitative agreement was found between the response of modeled systems and experimental results. Strong coupling between plastic yield and densification was observed. Our results also suggest that sodium silicates may densify not only under hydrostatic compression but also upon shear at large strains. Based on these numerical results, we propose a general yield criterion for soda-silicate glasses in which density is an internal variable. As density increases, the elliptic yield surface (characterizing amorphous silicates with open structures) gradually

*Corresponding author
Email address: gmolnar.work@gmail.com (Gergely Molnár)

Preprint submitted to Acta Materialia

March 18, 2016
evolves into a Drucker-Prager-like model for fully densified samples.

*Keywords:* Silicate, Constitutive modeling, Yield strength, Multiscale modeling, Molecular statics

# 1. Introduction

Silicate glasses are used for many technical purposes, especially where stiffness and transparency are required. Silicate glasses are brittle on the macroscopic scale but ductile at the micron scale [1]. This plastic response is expected to be key to understand brittleness. However, it was soon found that there is an unusual feature in the plastic response of amorphous silicates. Open structure glasses exhibit irreversible volumetric strain upon compression: for amorphous silica this densification saturates at ca. 20 % at a hydrostatic pressure of about 20 GPa [2, 3, 4, 5]. Moreover, for technical silicate glasses, sodium oxide ($Na_2O$) is usually added to silica, along with other compounds. Indeed, sodium modifies the silica network and lowers the glass transition temperature for easier glass processing. But the addition of sodium also impacts mechanical ductility. Upon hydrostatic compression soda-lime-silicates show less densification [6], and at reduced pressures values [7]. This behavior results from the phenomena that sodium gradually fills up the open structures.

In fact, the relation between densification and plastic yield in silicate glasses is a complex question, and a continuing matter of debate. Ever since densification was observed, the respective contribution of irreversible volumetric strain and shear flow has been an issue [8, 9]. Experimentally, because larger samples break, nano-indentation was a useful tool. Unfortunately the resulting strain field of an indentation is very inhomogeneous, which helps little in the identification of a complex constitutive behavior. According to a rule of thumb, if pile-up is present – especially with very sharp indenters [10] – there is significant shear plasticity [11]. While if densification is dominant, the refractive index of the glass changes [12]. A large number of experiments was carried out, and the overall picture emerged: if densification is prevented, shear flow is enhanced. It was shown, that the addition of sodium [13, 14] and pre-densification [15, 16] decreases volumetric plasticity and increases shear flow.

However, for a more quantitative assessment of the competing processes, it is necessary to develop a continuum scale description of the governing
plastic yield. Given the complexity of tests such as indentation, the use of advanced numerical tools and especially finite element modeling is inevitable. Compared to the large number of papers reporting indentation results, constitutive models for plastic yield in silicate glasses are few and far [17, 18, 19, 20]. Moreover, these models must be calibrated using experimental results.

During indentation test force displacement curves are the usual indicators, however they are not very sensitive. To strengthen the experimental basis, a more elaborate identification of plasticity have been considered such as the measurement of residual stress (from crack patterns [21] or birefringence [22]) or density fields after indentation (from Raman scattering [23], luminescence [24] or chemical etching [25]). However, for a precise description other tests must be devised. The “simplest” one is probably the uniaxial compression test [26, 27].

Once due consideration is given to these difficulties, it appears that numerical experiments may offer insight into the mechanical response of such complex materials. In this paper, our aim is to develop yield criteria from atomistic simulations for amorphous silicates, specifically considering the impact of composition. Given the present state of development of molecular dynamics, we do not expect these calculations to provide a quantitative characterization. Rather, we hope that our simulation results provide an explicit functional description of yield surface. Which model can then be made quantitative for given actual glass compositions, based on the various micron-scale mechanical experiments mentioned above, such as diamond anvil [28], pillar compression [26] or indentation tests [15].

In a pioneering work, Schuh and Lund [29] derived a constitutive relation from atomistic simulations for metallic glasses. Their numerically calculated yield surface compared favorably with experimental results. Since then, many works have been dedicated to measure plasticity in metallic glasses [29, 30, 31, 32, 33], nano-crystalline metals [34, 35] and glassy polymers [36, 37]. Amorphous solids in general were studied [38, 39, 40], though an elaborate quantitative description of silicate glasses is still missing. Simulations are very useful in this context, not only because they can provide detailed information about the mechanical response, but also because the underlying rearrangement mechanisms can be examined at the atom scale.

The paper is structured as follows. Section 2 introduces the numerical methods with details on molecular dynamics and statics [41]. In section 3 results are given, first for simple hydrostatic compression and shear, then for the combined of the two. This is necessary to identify the constitutive
behavior with strong coupling between volumetric and shear deformations. We show that the yield process can be rationalized when the evolution of density is considered. Finally in section 4, a new form of yield function is proposed with the density as an internal variable.

## 2. Methods

Our aim is to simulate the mechanical response of amorphous silicate materials with open structures with increasing depolymerization.

### 2.1. *Atomic sample preparation*

The glass samples were generated by the random sequential placement of the atoms in a cubic simulation box with periodic boundary conditions. Molecular simulations were performed with LAMMPS [42] to equilibrate the liquid at 3000 K, to cool it using 10¹³ K/s quenching rate until ~ 0 K, and finally deform the samples. Atomic interactions were modeled with the empirical BKS potential [43] using the parameters of Yuan and Cormack [44]. The potential function was defined for the different pair interactions (e.g. Si-O, O-O or Na-O) with respect to the strong ionic/covalent Si-O bonds and the weaker but longer Na-O ones<fnref n="1" />.

To avoid the collapse of the atoms the short term potential was substituted with a repulsive function [45]. In this manner we have simulated three glasses xNa₂O-(100-x)SiO₂ with x = 5, 15, and 30 %mol which will be referred to as NSx5, NSx15 and NSx30 following Yuan and Cormack [44]. The system sizes were 67 041, 69 849 and 73 368 atoms respectively with a final simulation box length of 10 nm. Thanks to this relatively large box, finite size effects are minimized. All samples were then compared with neutron [46], Brillouin [47] scattering and NMR [48, 49, 50] experiments to validate the initial structure. Agreement was found within the precision of the experimental results as shown also by Yuan and Cormack [44]. Further information about the sample generation can be found in Ref. [51].

<footnote n="1">¹The cutoff values used for the potential function were fine-tuned in order to achieve the experimentally measured densities (2.238, 2.340 and 2.470 g/cm³) at ambient pressure. Therefore, we used a parameter *r*<sub>cut</sub> equal to 5.1, 5.9, 6.9 Å for NSx5, NSx15 and NSx30 samples respectively.</footnote>
## 2.2. Deformation scheme

To display a pressure dependent yield surface the main variables are pressure (p) and equivalent shear stress (k). Where p = -σm, and the average normal stress σm is calculated from the components of the diagonal of the Cauchy stress tensor (σ): σm = (σ1 + σ2 + σ3)/3, σi are the principal normal stresses.

Similarly, equivalent shear stress k can be expressed as: $k = \sqrt{J_2}$, where $J_2$ is the second invariant of deviatoric stress tensor $s$ ($s = \sigma - \sigma_m I$, $I$ is the identity matrix). $J_2$ can be given by the components of the Cauchy stress tensor (σ) as well:

$$J_2 = \frac{1}{6} ((\sigma_1 - \sigma_2)^2 + (\sigma_1 - \sigma_3)^2 + (\sigma_2 - \sigma_3)^2). \quad (1)$$

In addition, we have also studied the impact of the third stress invariant $J_3 = \det(s)$. Therefore, we will also use the Haigh-Westergaard (HW) stress space [52, 53, 54, 55, 56, 57, 58], where the three independent variables are:

* hydrostatic stress $\rho = -p\sqrt{3}$,

* deviatoric stress $s = k\sqrt{2} = \sqrt{2J_2}$,

* meridian angle (or Lode angle) - $\vartheta$.

The meridian angle ($\vartheta$) is defined as [59]:

$$\cos 3\vartheta = \frac{3\sqrt{3} J_3}{2 J_2^{3/2}}. \quad (2)$$

Note the opposite sign convention for pressure and hydrostatic stress: in compression pressure is positive while hydrostatic stress is negative. The third variable ($\vartheta$) in the HW stress space indicates whether the deviatoric stress tensor is predominantly tensile ($\vartheta = 0^\circ$), shear ($\vartheta = 30^\circ$) or compressive ($\vartheta = 60^\circ$).

The deformation was applied in a quasi-static way [41]. During both compression and tension the dimensions of the simulation box was reduced by a constant displacement step, while the positions of the particles were rescaled in a homogeneous way. After box deformation, a new equilibrium position was searched using the Polak-Ribiere conjugate gradient algorithm. The macroscopic stress tensor of the system was calculated according to Ref. [60]:
$$ \boldsymbol{\sigma} = -\frac{1}{V} \sum_{i} \left[ -m_i \boldsymbol{v}_i \otimes \boldsymbol{v}_i + \frac{1}{2} \sum_{i \neq j} \boldsymbol{r}_{ij} \otimes \boldsymbol{f}_{ij} \right], \quad (3) $$

where i and j are the atomic indices. The summation is over all the particles contained by the total volume (V). The first term in the bracket is the kinetic energy contribution which is considered zero because samples are in static equilibrium. The variables $v_i$ and $m_i$ are the velocity vector and the mass of atom i. The second term is the stress originating from the pairwise energy. This value can be calculated for each atom by summing the tensorial product of $\boldsymbol{r}_{ij}$ ($\boldsymbol{r}_{ij} = \boldsymbol{r}_j - \boldsymbol{r}_i$, where $\boldsymbol{r}_i$ and $\boldsymbol{r}_j$ are the position vectors of atom i and j) and the interatomic force vector ($\boldsymbol{f}_{ij}$) applied on atom i by atom j.

To evaluate the constitutive relations, permanent plastic strains were calculated as a function of maximum stress state in various loading/unloading cycles. In practice, the yield stress is often written as a function of the plastic strains. Here, by the inverse process, we register plastic strains as a function of applied stresses to determine the hardening functions.

The following protocol was executed to calculate the permanent volume variations ($\epsilon_V^{pl}$) at a given stress state (P) after loading and unloading. In Fig. A.1 the schema of the applied deformation is shown. In a first step, the simulation box ($\rho = 0$ GPa) was deformed hydrostatically ($\sigma_1 = \sigma_2 = \sigma_3$) to obtain the desired pressure value. Then the sides of the box were elongated or compressed separately to apply deviatoric stress at constant pressure and constant meridian angle (step 2) to reach the target load point (P). During step 3 the direction of the deformation was reversed to reduce the applied deviatoric stress to zero, maintaining the principal stress relations as in step 2. Finally in step 4 the pressure was relaxed to zero. The final box shape was compared to the original one to calculate the permanent volumetric strain ($\epsilon_V^{pl}(\rho, s)$) as a function of applied stresses. The experiment was carried out for various stress states ($\rho = [-25, 25]$ GPa, $s = [0, 12]$ GPa) to map the plastic response of the material accurately. Fig. A.1 also displays a schema of the expected mechanical response: the red curve shows the quasi elastic domain with only small permanent strain, then the material should yield and after a gradual hardening the maximum load bearing capacity should be reached.

These tests were carried out for all three compositions not only on the pristine (undensified) samples but also after pre-densification as well. The samples were pre-densified by pure hydrostatic compression (step 1), then
<table><caption>Table 1: Parameter $r_i$ for different loading cases in eq. (4) [59, 61].</caption><thead><tr><th>Meridian section</th><th>ϑ</th><th>r<sub>1</sub></th><th>r<sub>2</sub></th><th>r<sub>3</sub></th></tr></thead><tbody><tr><td>Tensile</td><td>0°</td><td>2/√3</td><td>-1/√3</td><td>-1/√3</td></tr><tr><td>Pure shear</td><td>30°</td><td>1</td><td>0</td><td>-1</td></tr><tr><td>Compression</td><td>60°</td><td>1/√3</td><td>1/√3</td><td>-2/√3</td></tr></tbody></table>

the pressure was relaxation (step 4) to zero.

To ensure that deformation proceeds at a constant meridian angle, we followed a loading paths obeying:

$$ 
\begin{aligned}
\sigma_2 &= (\sigma_1 + p) r_2/r_1 - p, \\
\sigma_3 &= (\sigma_1 + p) r_3/r_1 - p,
\end{aligned}
\qquad (4) 
$$

where $r_1$, $r_2$ and $r_3$ are ratios given in Tab. 1. The term in the brackets is the stress difference applied by the deviatoric deformation compared to the original pressure state. For example if $\vartheta = 30^\circ$: the difference in the first principal stress is positive, while the difference in the third is negative and $\sigma_2$ remains unchanged.

If the material is elastic, a similar relationship for principal strains hold as well. However, when plasticity arises, the axial strains have to be modified separately. Fig. A.2 shows the evolution of both principal stresses and strains with increasing number of deformation steps. A constant strain increment $d\epsilon_1 = 10^{-4}$ (axial displacement $du_1=Ld\epsilon_1$) was used. Initially the other two axial strain increments were set as shown in Tab. 1 through $du_{2,3} = du_1r_{2,3}/r_1$. If the stresses were found to deviate from eq. (4) more than 0.005 GPa, the displacements were corrected in the next step. For example if the second principal stress at step $n$ ($\sigma_2^n$) was larger than requested, the displacement $du_2^{n+1}$ was reduced in the next step according to the feedback equation $du_2^{n+1} = du_1(r_2/r_1 + (\sigma_2 - \sigma_2^n)B)$, where $\sigma_2$ is the target value according to eq. (4) and $B$ is a convergence parameter set to 30 1/GPa.

Using this procedure we found that negligible rotation was induced ($\tau_{xy} = \tau_{xz} = \tau_{yz} < 0.005$ GPa) so that the axial stresses were directly equal to the principal stresses. As a result, to calculate the permanent volumetric strain, the normal strain components were calculated directly from the logarithm of the individual Hencky stretch ratios ($\lambda_i$):

$$ \varepsilon_i = \ln(\lambda_i) = \ln(1 + \varepsilon_i^{\text{SS}}). \qquad (5) $$
In eq. (5) $\epsilon_i$ are the finite and $\epsilon_i^{ss}$ are (small) engineering normal strains ($i = x, y, z$). The volumetric strain ($\epsilon_V$) can be calculated as the sum of the individual axial strain components:

$$ \epsilon_V = \sum_{i=x,y,z} \epsilon_i = \sum_{i=x,y,z} \ln(1 + \epsilon_i^{ss}) = \ln\left(1 + \frac{dV}{V_0}\right), \quad (6) $$

where, $dV$ is the volume change and $V_0$ is the original volume size.

Fig. A.2 shows how the principal strains had to be set in order to achieve the desired relationship between principal stresses for a $\vartheta = 30^\circ$ and $\rho = 0$ GPa ($p = 0$ GPa) deformation case. Using this iteration process we were able to conduct deviatoric (shear) tests at constant hydrostatic stress ($\rho$) and constant meridian angle ($\vartheta$).

### 3. Results

The results are divided into three sections. First we investigate pure hydrostatic compression and the resulting permanent volumetric strain as a function of pressure. Then we focus on the connection between deviatoric deformation and densification. Finally, the two basic deformation schemes (hydrostatic, deviatoric) are combined to map permanent volumetric deformation as a function of stresses to develop a yield surfaces. The dependence of these yield surfaces upon initial structures (composition or pre-densification) is highlighted. Finally we show how these results can lead to a parametrization as a function of the key internal variable: the density.

#### 3.1. Hydrostatic deformation

Simple hydrostatic loading experiments were carried out without shear. Fig. A.3 shows the permanent volumetric strain as a function of maximum hydrostatic pressure. The permanent volume variations were calculated by compressing the samples isotropically until the desired pressure value was reached, then the deformation was reversed in order to relax the pressure. The difference between the initial and the relaxed volume was compared to compute $\epsilon_V^{pl}$.

The response can be divided into three separate stages. In a first stage, at low pressure, the response is quasi-elastic and the volume change is relatively small. In a second stage, densification is roughly linear with pressure although the details are affected by Na content. More Na tends to reduce
<table><caption>Table 2: Parameters of eq. 10 for different compositions. R² is the coefficient of determination of the fit.</caption><thead><tr><th>Composition</th><th>ε<sub>V</sub><sup>pl,max</sup> [-]</th><th>p<sub>c</sub> [GPa]</th><th>m [-]</th><th>R<sup>2</sup></th></tr></thead><tbody><tr><td>NSx5</td><td>-0.31</td><td>14.61</td><td>3.37</td><td>0.9984</td></tr><tr><td>NSx15</td><td>-0.29</td><td>12.01</td><td>2.37</td><td>0.9989</td></tr><tr><td>NSx30</td><td>-0.21</td><td>12.97</td><td>2.12</td><td>0.9994</td></tr></tbody></table>

the densification threshold: NSx5 begins to densify at $p = 3$ GPa while NSx30 has almost no elastic resistance and plasticity starts at a very early stage ($p = 0.5$ GPa). In a third stage, the linear regime gradually bends over and saturates to a maximum value ($ε_V^{pl,max}$). The maximum permanent volumetric strain decreases with Na content.

More precisely the results can be adequately fitted using the following sigmoidal curve:

$$ε_V^{pl} = ε_V^{pl,max} \left( 1 - \frac{1}{1 + (p/p_c)^m} \right) \quad (7)$$

where $p_c$ and $m$ are material constants: $p_c$ is the center of the sigmoidal function close to the inflection point, where the increment of densification is the largest. Coefficient $m$ is the hardening exponent. The larger $m$, the larger the plastic strain increment near $p_c$. The material parameters $ε_V^{pl,max}$, $p_c$, $m$ are shown in Tab. 2.

Qualitatively, the evolution of the plastic properties with increasing Na content (reduction of the yield threshold and of the permanent volumetric strain at saturation) is in agreement with experimental results. Note however that the permanent volumetric strain at saturation is significantly larger than expected. The maximum volume loss for NSx30 is found at 16 %. However, experiments measured a much lower value for window glass around 6 % [6, 7]. This discrepancy is probably an artifact of the potential. Nevertheless, as usual with atomistic simulations our aim is not to predict material properties quantitatively, but to observe the atomistic mechanisms and the nature of the mechanical response. Then, the exact values of the response parameters should be determined using real life experiments.
## 3.2. Elementary shear deformation

In Fig. A.4 the deviatoric stress and the permanent volumetric strain of NSx5 are shown as a function of applied deviatoric strain<fnref n="2" /> ($\epsilon_s$) at constant hydrostatic tension $p = -2$ GPa. For this small negative pressure, pristine NSx5 is elastic up to $\epsilon_s = 0.15$ where it begins to densify gradually. After reaching its shear strength (peak stress) a softening stage appears, with densification. In this regime the plastic volumetric strain increases linearly with applied deviatoric strain. For large shear, a plastic plateau is observed, defining the flow stress. Under these conditions, densification even exceeds the maximum value reached with pure hydrostatic pressure. A possible reason for the saturation upon hydrostatic compression is that silica forms a stable load bearing network with 6 coordinated silicons at high pressure, whereas shear flow breaks this strong structures and allows the material to densify further [62]. Some samples were also partially densified using hydrostatic pressure before the same type of shear experiment was carried out at $p = -2$ GPa (Fig. A.4). We find that with pre-densification the material exhibits reduced strength, earlier plasticity and tends to expand at the initial deformation stage. Then we find a gradual transition, until the samples saturates at the same increased density $\epsilon_V^{pl,max}$ and the same flow stress on the plastic plateau. Thus, at large strains shear erases the structural modifications induced by pre-densification.

The softening phenomena shown in Fig. A.4 is generally explained by the localization of internal plasticity (or so-called shear transformation zones) [39]. These localized regions form shear bands which allow the material to lose elastic energy in an increased way, thus reducing the shear strength of the material. Our tests connecting densification and softening is not a contradictory hypothesis, but a result which shows, that in silicates this shear band causes not only softening but permanent volume loss as well. The detailed analysis of the localization of internal plasticity exceeds the topic of this paper and does not affects the final results concerning the macroscopic yield properties.

## 3.3. Determination of yield criteria

The central aim of this study is to define a yield surface for sodium silicate-like materials as a function of composition. In the previous sections we found

<footnote n="2">²The deviatoric strain was calculated from the Hencky strain tensor in a manner similar to deviatoric stress</footnote>
a strong correlation between densification and yield strength. Therefore, to
investigate the pressure and densification dependence of the yield strength
several tests were carried out under combined hydrostatic and deviatoric
loading (the method is explained in section 2).

Fig. A.5 shows permanent volumetric strain as a function of maximum stress state in the pure shear meridian ($\vartheta = 30^\circ$) for each of the three compositions. The squares show the maximum stress that can be reached during deviatoric deformation due to the plastic plateau or even softening (see Fig. A.4). It can be noticed that both shear strength and hydrostatic tensile strength are reduced by the addition of sodium. The elastic domain is reduced, which suggests that sodium weakens the sample and facilitates early plastification. More importantly, for low Na content, a minimum in the deviatoric strength is found in the early stage: for NSx5 the local minimum is at $\rho = -10$ GPa ($p = 5.77$ GPa), for NSx15 is at $\rho = -3$ GPa ($p = 1.73$ GPa). After this local minimum, the deviatoric strength increases again with pressure, similar to a Drucker-Prager model [63]. Interestingly NSx30 does not show any sign of a local minimum and the deviatoric strength increases monotonically with increasing pressure. Hence, for high sodium content, an extended Drucker-Prager model [64] is suitable to account for the yield strength evolution.

The previous representation suffers from several drawbacks. The permanent volumetric strain can be shown only up to the maximum stress value. For a softening material, such as shown in Fig. A.4, the softening stage cannot be accounted for. More significantly, with the local minimum, shear strength curve is concave, which appears to violate Drucker’s postulate [65] and leads to instability. These observations suggest that the material actually evolves during plastic deformation and that parametrization by some internal variable is necessary. From Fig. A.4 it is clear that densification lowers the deviatoric strength. To clarify the situation, we have first pre-densified samples using hydrostatic pressure, then performed the same density mapping procedure through combined pressure-shear loading. In Fig. A.6 the maximum deviatoric stress is shown as a function of relative densification $\gamma = \varepsilon_V^{pl} / \varepsilon_{V,\max}^{pl,max}$ where $\varepsilon_V^{pl}$ is the permanent volumetric strain and $\varepsilon_{V,\max}^{pl,max}$ is the maximum value reachable using hydrostatic compression. In the same figure, isodensification curves taken from Fig. A.5 illustrate the transformation of the initial yield surface into the final one: as permanent volumetric strain increases, the positive yield pressure increases, the deviatoric yield
strength decreases and the yield curve becomes flatter. Most importantly, Fig. A.6 shows that once the dependence upon density has been explicitly taken into account, the yield surfaces for each given density are convex so that Drucker's postulate is indeed satisfied [65].

The evolution of the yield surfaces does not only depend upon pre-densification but also upon composition. With increasing sodium content the difference between initial and final deviatoric yield strength reduces: for NSx5 yield strength is strongly affected by densification while for NSx30 densification has almost no effect on maximum strength. However, in all cases the final yield strength increases linearly at high pressure, but curves down to zero as a power law near the tensile threshold.

We also found that the third stress invariant (meridian angle) affects deviatoric strength. In Fig. A.7 the yield strength is shown at pressure $p = 0$ GPa (in the $\pi$-plane) for the three different compositions. For pristine samples ($\gamma \approx 0\%$, in red) the plots are not circular, which denotes an impact of meridian angle on shear strength. The shear strength is lower under compression dominated shear ($\vartheta = 60^\circ$) than under tension dominated shear ($\vartheta = 0^\circ$). For fully densified samples ($\gamma = 100\%$), The strength curves at $p = 0$ GPa are shown in black. Upon densification, for NSx5, the deviatoric strength remains almost unchanged. This is an accident connected to this specific pressure value $p = 0$ GPa, as can be seen in Fig. A.6. With more sodium, pre-densification not only reduces deviatoric strength but also suppresses anisotropy of the yield surface in the $\pi$-plane. Generally speaking, we find that yield stress in shear is lower for a compression dominated stress state ($\vartheta = 60^\circ$) and higher for a tension dominated stress state ($\vartheta = 0^\circ$). This response can tentatively be related to the impact of compression which is known to lower yield stress in open structure materials. Indeed, for fully densified samples, the yield stress becomes independent of meridian angle, as in the Drucker-Prager scheme. This observation is consistent with the idea that the response of these materials generally transforms from a compression to a deviatoric stress sensitive yield criteria as the material densifies.

## 4. Constitutive model

Using microscopic experiments, such as micro-pillar compression and indentation, Kermouche et al. [19] showed that silica can be modeled appropriately using an elliptic yield surface (i.e. a prolate spheroidal surface in HW space). Permanent densification gradually transforms this initial yield
function into a pressure insensitive von Mises criterion [66] (i.e. a cylinder in HW space). Molecular dynamics simulation on silica-like materials have confirmed the overall curvature of the yield surface [67].

Our present atomistic results confirm this hypothesis. They also suggest a more appropriate model for amorphous silicates. Our simulations show that the yield function transforms gradually from an initial elliptic to a Drucker-Prager-like shape as densification increases. As a result we can propose a density dependent yield function to capture this evolution as a function of permanent volumetric strain (Fig. A.8a).

The yield function (referred as DP-cap in the following) is compiled from two parts: an extended Drucker-Prager-like model and an elliptic cap. The tension side, from $\rho^{y,+}$ to $\rho^{int}$ is described with the extended Drucker-Prager function [64]:

$$F_{PD} = \left( \frac{s}{c(\gamma)} \right)^{b(\gamma)} + \frac{\rho - \rho^{y,+}(\gamma)}{a(\gamma)}. \qquad (8)$$

In eq. (8) $1/a$ is similar to the slope found in the original Drucker-Prager model [63], $b$ and $c$ are empirical material constants and $\rho^{y,+}$ is the hydrostatic tensile strength. All material constants are linear functions of relative permanent volume change ($\gamma = \epsilon_V^{pl} / \epsilon_V^{pl,max}$, where $\epsilon_V^{pl}$ is the volumetric plastic strain).

The elliptic cap is defined for the compression side between $\rho^{int}$ and $\rho^{y,-}$ as follows:

$$F_{Ell} = \left( \frac{\rho - h}{d} \right)^2 + \left( \frac{s}{e} \right)^2 - 1, \qquad (9)$$

where $d$, $e$ are the radii and $h$ is the center of the ellipse on the hydrostatic axis. The cap is always fitted in connection to the power function by the following three criteria: 1) the values of both functions (8) and (9) have to be equal at $\rho_{int}$; 2) the connection needs to be smooth, therefore the derivatives of $s$ respect to $\rho$ for both functions have to be continuous at $P(\rho^{int}, s^{int})$ and finally $\rho^{y,-} = h - d$, where $\rho^{y,-}$ is the compressive yield pressure calculated as the inverse of equation (7):

$$\rho^{y,-} = -\sqrt{3}p_c \left( \frac{\gamma}{1-\gamma} \right)^{1/m}. \qquad (10)$$
<table><caption>Table 3: Material constants for the DP-cap model as a function of composition and relative permanent volumetric strain. The curves were fitted to the results presented in Fig. A.6.</caption><thead><tr><th>Composition</th><th colspan="2">NSx5</th><th colspan="2">NSx15</th><th colspan="2">NSx30</th></tr><tr><th>γ [%]</th><th>0</th><th>100</th><th>0</th><th>100</th><th>0</th><th>100</th></tr></thead><tbody><tr><td>ρ<sup>y,+</sup> [GPa]</td><td>22.79</td><td>16.94</td><td>15.31</td><td>13.16</td><td>2.51</td><td>6.00</td></tr><tr><td>ρ<sup>y,-</sup> [GPa]</td><td colspan="6">Shown by eq. (10) and Tab. 2</td></tr><tr><td>ρ<sup>int</sup> [GPa]</td><td>20.0 + ρ<sup>y,-</sup></td><td>15.0 + ρ<sup>y,-</sup></td><td>6.5 + ρ<sup>y,-</sup></td><td></td><td></td><td></td></tr><tr><td>a [GPa]</td><td colspan="2">13.0</td><td colspan="2">20.0</td><td colspan="2">24.0</td></tr><tr><td>b [-]</td><td>1.5</td><td>4.4</td><td>1.5</td><td>3.7</td><td>1.5</td><td>3.4</td></tr><tr><td>c [GPa]</td><td>14.0</td><td>7.0</td><td>20.0</td><td>7.0</td><td>24.0</td><td>6.4</td></tr></tbody></table>

Of course, although the results of the numerical simulations were best fitted with function (10), this equation can be replaced with any other sigmoidal curve used in the literature [6, 20] if it fits experimental results better. From the three independent criteria, the three parameters of the ellipse can be calculated. Computational details concerning the relationship between the elliptic cap and the power function can be found in Appendix A.

The DP-cap yield function changes only with relative densification, which can result both from hydrostatic and deviatoric loading. It handles hydrostatic hardening as well as shear softening at constant pressure (Fig. A.4).

Fig. A.8b shows the DP-cap model fitted on the results obtained for NSx5. All material constants are summarized in Tab. 3.

The advantage of the DP-cap yield function is that it is capable of modeling not only one directional but cyclic loading as well, using density as an internal variable. As to the tensile side of the ellipse, for samples with small sodium content, negative volume variation can be applied in a very restricted way only. The sample loses its load bearing capacity very early and then a macroscopic crack appears, from which point tensile side of the strength curve was defined. Of course it is doubtful that this model can describe rupture adequately. Therefore, to give a physical meaning to parameters (b and c) which affect the tensile side of the power law is difficult. To provide a better brittle response to the model, a different method (e.g. phase-field method, XFEM or a hybrid discrete and finite element method) should be used in parallel to this yield criterion. Note also that, to a lesser extent, the samples are affected by irreversible shear strain and the meridian angle as well. Therefore, both deviatoric strain hardening and the Lode parameter should be taken into account for a better representation of the atomistic results.
5. Conclusion and remarks on future work

Molecular statics simulations were performed to investigate the mechan-
ical response of silicate glasses with increasing network depolymerization.
Multiaxial loading conditions were used, where the samples were first de-
formed hydrostatically, then sheared. After relaxing the pressure and shear
stress to zero, permanent volumetric strains were computed from initial and
final volumes. The results show that densification not only affects hydro-
static hardening, but also has significant impact on shear yield strength as
well. Sizable loss of free volume is also registered for large amplitude shear
deformations. Under these conditions, densification can be interpreted as a
microstructural change caused by the collapse of large rings upon shear [62].

We find that shear strength is larger for more polymerized samples (i.e. with low sodium content). Larger amounts of sodium facilitate early plasticization, with less softening and more ductility. This result is consistent with the observation that sodium rich silicates undergo more shear flow and less densification. Low sodium content samples enter a softening regime controlled by volumetric strain while sodium rich silicates can reach their shear strength at zero pressure without densification. However, a similar transition can also be achieved in low sodium samples by pre-densification: from an elliptic model, the response transforms into a Drucker-Prager behavior when density increases, confirming the tendency observed in experiments [15].

Based on our numerical results, we have proposed a density dependent
yield function. Considering the complex hardening of the material, density is
introduced as an internal variable taking into account the permanent volume
changes for both hydrostatic and shear deformation. The yield function is
initially elliptic and gradually evolves into a Drucker-Prager criterion when
densification saturates.

In principle, using numerical experiments to derive constitutive equa-
tions, as presented in this paper, has two advantages: 1) it is easier and
much more cost efficient to conduct experiments virtually; 2) for silicate
glasses, it is difficult/impossible to conduct multiaxial tests experimentally.
For instance, yield stress as a function of meridian angle can certainly not
be measured in the present state of experimental development. Therefore,
simulations should be a powerful tool to derive shapes of yield functions and
identify essential material parameters. However, the technique is of course
still very limited by the problem of definition and calibration of the interac-
tion potentials. For example, in our results, although the potential provides
accurate static structures, the density at saturation after deformation cer-
tainly overestimates expected experimental values. This is why, at present,
conclusions can be derived only in a qualitative way. To actually implement
such a constitutive model, it is necessary to determine material constants by
comparing the simulated data with a variety of experimental results, under
as many different loading conditions as possible.

Acknowledgment

This work was supported by the French Research National Agency pro-
gram MultiSil (ANR-13-BS09-0012).

References

[1] D. M. Marsh, Plastic flow in glass, Proceedings of the Royal Society of London, Ser. A 279 (1378) (1964) 420 – 435.

[2] J. D. Mackenzie, High-pressure effects on oxide glasses: I. densification in rigid state, Journal of the American Ceramic Society 46 (10) (1963) 461–470.

[3] D. Vandembroucq, T. Deschamps, C. Coussa, A. Perriot, E. Barthel, B. Champagnon, C. Martinet, Density hardening plasticity and mechanical ageing of silica glass under pressure: a Raman spectroscopic study, Journal of Physics: Condensed Matter 20 (48) (2008) 485221.

[4] T. Deschamps, A. Kassir-Bodon, C. Sonneville, J. Margueritat, C. Martinet, D. de Ligny, A. Mermet, B. Champagnon, Permanent densification of compressed silica glass: a Raman-density calibration curve, Journal of Physics: Condensed Matter 25 (2) (2013) 025402.

[5] T. Rouxel, Driving force for indentation cracking in glass: composition, pressure and temperature dependence, Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 373 (2038) (2015).

[6] H. Ji, V. Keryvin, T. Rouxel, T. Hammouda, Densification of window glass under very high pressure and its relevance to vickers indentation, Scripta Materialia 55 (12) (2006) 1159–1162.
[7] T. Deschamps, C. Martinet, J. L. Bruneel, B. Champagnon, Soda-lime silicate glass under hydrostatic pressure and indentation: a micro-Raman study, Journal of Physics: Condensed Matter 23 (3) (2011) 035402.

[8] F. M. Ernsberger, Mechanical properties of glass, Journal of Non-Crystalline Solids 25 (13) (1977) 293-321.

[9] F. M. Ernsberger, Glass Science and Technology: Elasticity and Strength in Glasses, Academic Press, 1980.

[10] S. Yoshida, H. Sawasato, T. Sugawara, Y. Miura, J. Matsuoka, Effects of indenter geometry on indentation-induced densification of soda-lime glass, Journal of Materials Research 25 (2010) 2203-2211.

[11] K. W. Peter, Densification and flow phenomena of glass in indentation experiments, Journal of Non-Crystalline Solids 5 (2) (1970) 103-115.

[12] F. M. Ernsberger, Role of densification in deformation of glasses under point loading, Journal of the American Ceramic Society 51 (10) (1968) 545-547.

[13] C. R. Kurkjian, G. W. Kammlott, M. M. Chaudhri, Indentation behavior of soda-lime silica glass, fused silica, and single-crystal quartz at liquid nitrogen temperature, Journal of the American Ceramic Society 78 (3) (1995) 737-744.

[14] T. M. Gross, M. Tomozawa, Fictive temperature-independent density and minimum indentation size effect in calcium aluminosilicate glass, Journal of Applied Physics 104 (6) (2008) 1-11.

[15] T. Rouxel, H. Ji, J. P. Guin, F. Augereau, B. Rufflé, Indentation deformation mechanism in glass: Densification versus shear flow, Journal of Applied Physics 107 (9) (2010) 094903.

[16] K. Aakermann, K. Januchta, J. Pedersen, M. Svenson, S. Rzoska, M. Bockowski, J. Mauro, M. Guerette, L. Huang, M. Smedskjaer, Indentation deformation mechanism of isostatically compressed mixed alkali aluminosilicate glasses, Journal of Non-Crystalline Solids 426 (2015) 175-183.
[17] M. Imaoka, I. Yasui, Finite element analysis of indentation on glass, Journal of Non-Crystalline Solids 22 (2) (1976) 315–329.

[18] J. C. Lambropoulos, S. Xu, T. Fang, Constitutive law for the densification of fused silica, with applications in polishing and microgrinding, Journal of the American Ceramic Society 79 (6) (1996) 1441–1452.

[19] G. Kermouche, E. Barthel, D. Vandembroucq, P. Dubujet, Mechanical modelling of indentation-induced densification in amorphous silica, Acta Materialia 56 (13) (2008) 3222–3228.

[20] V. Keryvin, J.-X. Meng, S. Gicquel, J.-P. Guin, L. Charleux, J.-C. Sanglebœuf, P. Pilvin, T. Rouxel, G. L. Quilliec, Constitutive modeling of the densification process in silica glass under hydrostatic compression, Acta Materialia 62 (2014) 250–257.

[21] R. F. Cook, G. M. Pharr, Direct observation and analysis of indentation cracking in glasses and ceramics, Journal of the American Ceramic Society 73 (4) (1990) 787–817.

[22] S. Yoshida, S. Iwata, T. Sugawara, Y. Miura, J. Matsuoka, A. Errapart, C. R. Kurkjian, Elastic and residual stresses around ball indentations on glasses using a micro-photoelastic technique, Journal of Non-Crystalline Solids 358 (24) (2012) 3465–3472.

[23] A. Perriot, D. Vandembroucq, E. Barthel, V. Martinez, L. Grosvalet, C. Martinet, B. Champagnon, Raman microspectroscopic characterization of amorphous silica plastic behavior, Journal of the American Ceramic Society 89 (2) (2006) 596–601.

[24] A. Perriot, E. Barthel, G. Kermouche, G. Quérel, D. Vandembroucq, On the plastic deformation of soda-lime glass: a Cr³⁺ luminescence study of densification, Philosophical Magazine 91 (7-9) (2011) 1245–1255.

[25] Y.-F. Niu, K. Han, J.-R. Guin, Locally enhanced dissolution rate as a probe for nanocontact-induced densification in oxide glasses, Langmuir 28 (29) (2012) 10733–10740.

[26] R. Lacroix, G. Kermouche, J. Teisseire, E. Barthel, Plastic deformation and residual stresses in amorphous silica pillars under uniaxial loading, Acta Materialia 60 (15) (2012) 5555–5566.
[27] D. Wakabayashi, N. Funamori, T. Sato, Enhanced plasticity of silica glass at high pressure, Physical Review B 91 (2015) 014106.

[28] D. Wakabayashi, N. Funamori, T. Sato, T. Taniguchi, Compression behavior of densified $SiO_2$ glass, Physical Review B 84 (2011) 144103.

[29] C. A. Schuh, A. C. Lund, Atomistic basis for the plastic yield criterion of metallic glass, Nature Materials 2 (2003) 449–452.

[30] A. C. Lund, C. A. Schuh, Yield surface of a simulated metallic glass, Acta Materialia 51 (18) (2003) 5399–5411.

[31] A. C. Lund, C. A. Schuh, The Mohr-Coulomb criterion from unit shear processes in metallic glass, Intermetallics 12 (1011) (2004) 1159–1165.

[32] A. C. Lund, C. A. Schuh, Strength asymmetry in nanocrystalline metals under multiaxial loading, Acta Materialia 53 (11) (2005) 3193–3205.

[33] F. Shimizu, S. Ogata, J. Li, Yield point of metallic glass, Acta Materialia 54 (16) (2006) 4293–4298.

[34] H. van Swygenhoven, M. Spaczer, A. Caro, D. Farkas, Competing plastic deformation mechanisms in nanophase metals, Physical Review B 60 (1999) 22–25.

[35] A. C. Lund, T. G. Nieh, C. A. Schuh, Tension/compression strength asymmetry in a simulated nanocrystalline metal, Physical Review B 69 (2004) 012101.

[36] P. H. Mott, A. S. Argon, U. W. Suter, Atomistic modelling of plastic deformation of glassy polymers, Philosophical Magazine A 67 (4) (1993) 931–978.

[37] J. Rottler, M. O. Robbins, Yield conditions for deformation of amorphous polymer glasses, Physical Review E 642 (2001) 051801.

[38] M. L. Falk, C. E. Maloney, Simulating the mechanical response of amorphous solids using atomistic methods, The European Physical Journal B 75 (4) (2010) 405–413.
[39] D. Rodney, A. Tanguy, D. Vandembroucq, Modeling the mechanics of amorphous solids at different length scale and time scale, Modelling and Simulation in Materials Science and Engineering 19 (8) (2011) 083001.

[40] Y. Xu, J. Chen, H. Li, Finite hyperelastic-plastic constitutive equations for atomistic simulation of dynamic ductile fracture, International Journal of Plasticity 59 (2014) 15–29.

[41] D. Theodorou and, U. W. Suter, Atomistic modeling of mechanical properties of polymeric glasses, Macromolecules 19 (1) (1986) 139–154.

[42] S. Plimpton, Fast parallel algorithms for short-range molecular dynamics, Journal of Computational Physics 117 (1995) 1–19.

[43] B. W. H. van Beest, G. J. Kramer, R. A. van Santen, Force fields for silicas and aluminophosphates based on ab initio calculations, Physical Review Letters 64 (1990) 1955–1958.

[44] X. Yuan, A. N. Cormack, Local structures of MD-modeled vitreous silica and sodium silicate glasses, Journal of Non-Crystalline Solids 283 (13) (2001) 69–87.

[45] B. Mantisi, A. Tanguy, G. Kermouche, E. Barthel, Atomistic response of a model silica glass under shear and pressure, The European Physical Journal B 85.

[46] M. Fábián, P. Jóvári, E. Sváb, G. Mészáros, T. Proffen, E. Veress, Network structure of 0.7SiO$_{2}$-0.3Na$_{2}$O glass from neutron and X-ray diffraction and rmc modelling, Journal of Physics: Condensed Matter 19 (33) (2007) 335209.

[47] Q. Zhao, M. Guerette, G. Scannell, L. Huang, In-situ high temperature Raman and Brillouin light scattering studies of sodium silicate glasses, Journal of Non-Crystalline Solids 358 (2012) 3418–3426.

[48] J. F. Emerson, P. E. Stallworth, P. J. Bray, High-field $^{29}$Si NMR studies of alkali silicate glasses, Journal of Non-Crystalline Solids 113 (1989) 253–259.

[49] H. Maekawa, T. Maekawa, K. Kawamura, T. Yokokawa, The structural groups of alkali silicate glasses determined from $^{29}$Si MAS-NMR, Journal of Non-Crystalline Solids 127 (1) (1991) 53–64.
[50] T. Charpentier, S. Ispas, M. Profeta, F. Mauri, C. J. Pickard, First-principles calculation of $^{17}$O, $^{29}$Si, and $^{23}$Na NMR spectra of sodium silicate crystals and glasses, The Journal of Physical Chemistry B 108 (13) (2004) 4147-4161.

[51] G. Molnár, P. Ganster, J. Török, A. Tanguy, Sodium effect on static mechanical behavior of md-modeled sodium silicate glasses, Journal of Non-Crystalline Solids 440 (2016) 12-25.

[52] B. P. Haigh, The strain energy function and the elastic limit, Engineering 109 (1920) 158-160.

[53] H. M. Westergaard, On the resistance of ductile materials to combined stresses in two or three directions perpendicular to one another, Journal of The Franklin Institute 189 (1920) 627-640.

[54] W.-F. Chan, D.-J. Han, Plasticity for Structural Engineers, Springer-Verlag, 1988.

[55] G. A. Holzapfel, Nonlinear Solid Mechanics, Wiley, 2001.

[56] M.-H. Yu, Generalized Plasticity, Springer Science and Business Media, 2006.

[57] D. Bigoni, Nonlinear Solid Mechanics: Bifurcation Theory and Material Instability, Cambridge University Press, 2012.

[58] A. Puzrin, Constitutive Modelling in Geomechanics, Springer-Verlag, 2012.

[59] K. Danas, P. P. Castañeda, Influence of the lode parameter and the stress triaxiality on the failure of elasto-plastic porous materials, International Journal of Solids and Structures 49 (1112) (2012) 1325-1342.

[60] A. P. Thompson, S. J. Plimpton, W. Mattson, General formulation of pressure and stress tensor for arbitrary many-body interaction potentials under periodic boundary conditions, The Journal of Chemical Physics 131 (15).

[61] J. Hutchinson, V. Tvergaard, Comment on "Influence of the Lode pa- rameter and the stress triaxiality on the failure of elasto-plastic porous
material” by K. Danas and P. Ponte Castañeda, International Journal of Solids and Structures 49 (23) (2012) 3484–3485.

[62] N. S. Shcheblanov, B. Mantisi, P. Umari, A. Tanguy, Detailed analysis of plastic shear in the raman spectra of SiO₂ glass, Journal of Non-Crystalline Solids 2428 (2015) 6–19.

[63] D. C. Drucker, W. Prager, Soil mechanics and plastic analysis for limit design, Quarterly of Applied Mathematics 10 (2) (1952) 157–165.

[64] ABAQUS, ABAQUS Documentation, Dassault Systemes, Providence, RI, USA, 2011.

[65] D. C. Drucker, A definition of a stable inelastic material, ASME Journal of Applied Mechanics 26 (1959) 1011–195.

[66] R. von Mises, Mechanik der festen Körper im plastisch deformablen Zustand, Königliche Gesellschaft der Wissenschaften zu Göttingen 1 (1913) 582–592.

[67] B. Mantisi, G. Kermouche, E. Barthel, A. Tanguy, Impact of pressure on plastic yield in amorphous solids with open structure, Physical Review E 93 (2016) 033001.
## Appendix A. Cap fitting

To fit the three parameters of the elliptic cap to the power function, three independent criteria is used:

* The deviatoric stress values of both elliptic and power function have to be equal at the intersection point ($\rho_{int}$): $s_{Ell} = s_{PD}$.

* The connection needs to be smooth, therefore the derivatives of s respect to $\rho$ for both functions again have to be the same at $P(\rho^{int}, s^{int})$.

* Finally $\rho^{y,-} = h - d$, where $\rho^{y,-}$ is the compressive yield pressure calculated as the inverse of equation (7).

To satisfy the first criteria the deviatoric stress value ($s^{int}$) at $\rho^{int}$ can be calculated from equation (8). Then both $\rho^{int}$ and $s^{int}$ has to be substituted into equation (9):

$$ \left( \frac{\rho^{\text{int}} - h}{d} \right)^2 + \left( \frac{s^{\text{int}}}{e} \right)^2 = 1, \qquad (A.1) $$

where $\rho^{int}$ and $s^{int}$ are known and $e, d$ and $h$ are the unknowns.

The derivative of the power function (8) at the intersection point can be expressed as:

$$ \frac{ds_{PD}}{d\rho} = -\frac{c}{ab} \left( \frac{\rho^{y,+} - \rho^{\text{int}}}{a} \right)^{1/b-1}. \qquad (A.2) $$

Using implicit differentiation the same slope can be expressed for the elliptic cap as well:

$$ \frac{ds_{Ell}}{d\rho} = -\frac{e^2 \rho^{\text{int}} - h}{d^2 s^{\text{int}}}. \qquad (A.3) $$

Knowing that (A.2) and (A.3) has to be equal, the second equation can be expressed. Using the third condition ($\rho^{y,-} = h - d$), the three parameters of the ellipse ($e, d$ and $h$) can be calculated from three independent equations.
<figure><img src="image_4.png" /><figcaption>Figure A.1: To calculate permanent volumetric strains at different stress states, the samples were first deformed hydrostatically (step 1), then pure deviatoric stress was applied following eq. (4) (step 2) and finally deviatoric stress (step 3) and pressure (step 4) were relaxed to zero. The final shape of the simulation box was compared to the original one to compute permanent volumetric strain $\epsilon_V^{pl}$ as a function of applied stresses $\rho, s$. The figure also shows both initial (quasi-elastic region) and final (fully densified) yield strength. The thin lines show the evolution of the yield surface during gradual densification upon hydrostatic and shear deformation.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure A.2: Principal stresses and principal strains as a function of load step for deformation in the pure shear meridian ($\vartheta = 30^\circ$). A constant strain increment $d\epsilon_1 = 10^{-4}$ was used and the other two principal strain steps ($d\epsilon_2$ and $d\epsilon_3$) were iteratively adjusted to match the desired relations between principal stresses (eq. (4)).</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure A.3: Permanent volumetric strain ($\epsilon_V^{pl}$) as a function of applied hydrostatic pressure for different compositions. The results are fitted using a sigmoidal fit (see eq. (7)). The parameters of the fit can be found in Tab. 2.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure A.4: Top: deviatoric stress ($k$) as a function of applied deviatoric strain ($\epsilon_s$) under constant hydrostatic tension $p = -2$ GPa for NSx5. The results are shown for three levels of pre-densification defined by the initial permanent volumetric strain ($\epsilon_V^{pl,ini}$). Bottom: permanent volumetric strain $\epsilon_V^{pl}$. The results demonstrate that shear at constant pressure induces densification at large shear deformations.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure A.5: Residual volumetric strain as a function of applied pressure (or hydrostatic stress) and deviatoric stress for three different compositions (see section 2). The open symbols shown the maximum yield strength reached. In the region in red permanent volume increase for NSx30 was found. Other compositions initially did not show any permanent plastic extension. The quasi-elastic regime is found in region in white.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure A.6: Interpolated curves using the results presented in Fig. A.5. All yield surfaces are shown as a function of relative densification ($\gamma = \epsilon_V^{pl}/\epsilon_V^{pl,max}$) at $\vartheta = 30^{\circ}$ (pure shear meridian). The dotted lines show intermediate stages of volumetric densification. The solid lines with different symbols represent initially densified samples (see Fig. A.4).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure A.7: Maximum yield surface shown in the π-plane ($p = 0$ GPa) as a function of densification for different compositions. Red curve shows the yield strength for pristine samples ($\gamma = 0\%$). Black curve shows the yield strength for fully densified samples ($\gamma = 100\%$). Each symbol shows a simulated case.</figcaption></figure>
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>Figure A.8: a) Schematic illustration of the proposed yield function which combines extended Drucker-Prager model (see eq. (8) – in black) and elliptic cap (see eq. (9) – in red). The connecting point is shown by $\rho^{int}$ and $s^{int}$. The positive and the negative yield point on the hydrostatic stress axis are shown by $\rho^{y,+}$ and $\rho^{y,-}$. b) Fitted yield functions for NSx5 for the initial, undensified state ($\gamma = 0\%$; model: solid line; numerical results: squares) and final, fully densified state ($\gamma = 100\%$; model: dashed line; numerical results: triangles). Dotted lines show the gradual transformation between the two.</figcaption></figure>