text:
# Energy distribution analysis of the wavepacket simulations of CH4 and CD4 scattering

## 1 Introduction

The dissociation of methane on transition metals is an important reaction in catalysis. It is the rate limiting step in steam reforming to produce syngas. It is also prototypical for C–H activation in other processes. A large number of molecular beam experiments in which the dissociation energy was measured as a function of translational energy have already been done on this system. These experiments have contributed much to our understanding of the mechanism of the dissociation. Some of them observed that vibrationally hot $`\mathrm{CH}_4`$ dissociates more readily than cold $`\mathrm{CH}_4`$, with the energy in the internal vibrations being about as effective as the translational energy in inducing dissociation. A more detailed assessment of the importance of the internal vibrations could not be made, because of the large number of internal vibrations. A recent molecular beam experiment with laser excitation of the $`\nu _3`$ mode succeeded in measuring a dramatical enhancement of the dissociation on a Ni(100) surface, but it is still much too low to account for the vibrational activation observed in previous studies and indicates that other vibrationally excited modes contribute significantly to the reactivity of thermal samples.

Wavepacket simulations are being used more and more to study the dynamics of this kind of molecule surface reactions. The published wavepacket simulations on the methane dissociation reaction on transition metals have treated the methane molecule always as a diatomic up to now. Besides the C–H bond and molecule surface distance, a combination of other coordinates were included, like (multiple) rotations and some lattice motion. None of them have looked at the role of the internal vibrations. Various theoretical studies have obtained reaction pathways and barriers for dissociation by DFT calculations, but they cannot explain the role of the vibrational modes in the reaction dynamics either.

A nice way to study reaction dynamics is the use of isotopes. The most recent wavepacket simulation on the dissociation probability of CH4 and CD4 showed a semiquantitative agreement with the molecular beam experiments of Ref., except for the isotope effect and the extracted vibrational efficacy. The molecular beam study with laser excitation of the $`\nu _3`$ asymmetrical stretch mode shows that the incorrect vibrational efficacy is caused by the assumptions in the fit procedure that both stretch modes behaves identical. One of the possible explanation of the incorrect isotope effect can be the role played by the non-included intramolecular vibrations.

In a previous paper we reported on wavepacket simulations that we have done to determine which and to what extent internal vibrations are important for the dissociation of CH4. We were not able yet to simulate the dissociation including all internal vibrations. Instead we simulated the scattering of methane, for which all internal vibrations can be included, and used the results to deduce consequences for the dissociation. We used model potential energy surfaces (PESs) that have been developed with Ni(111) in mind, but our results should hold for other surfaces as well. At a translational energy up to 96 kJ/mol we found that the scattering is almost completely elastic. Vibrational excitations when the molecule hits the surface and the corresponding deformation depend on generic features of the potential energy surface. In particular, our simulations indicate that for methane to dissociate the interaction of the molecule with the surface should lead to an elongated equilibrium C–H bond length close to the surface.

We have been using the multiconfigurational time-dependent Hartree (MCTDH) method for our wavepacket simulation, because it can deal with a large number of degrees of freedom and with large grids. This method has been applied successfully to gas phase reactions and reactions at surfaces.

In this paper we report wavepacket simulations of CD4 scattering including all internal vibrations for fixed orientations, performed on the same model PESs as in our previous paper. Translational motion parallel to the surface and all rotational motion was neglected. No degrees of freedom of the surface were included. Experiments show that coupling with these degrees of freedom is dependent on the metal surface. For example, the observed surface temperature effect are small on Ni, but quite large on Pt. As we are only interested in the role of internal vibrations, we have not included degrees of freedom of the surface to keep the simulations as simple as possible. We will discuss the vibrational excitation and the deformation of the CD4 molecule when it hits the surface and compare it with CH4. Later on we will look at the energy distribution of the kinetic energy per mode and the potential energy in some terms of the PES with the elongated equilibrium bond length close to the surface for both isotopes. The transfer of translational kinetic energy towards vibrational kinetic energy gives an indication about the dissociation probability, since vibrational kinetic energy helps in overcoming the dissociation barrier. It gives a better idea too about which modes are essential to include in a more accurate wavepacket simulation of methane dissociation. After that we will discuss the implications of this for the dissociation and give a summary with some general conclusions.

## 2 Computational details

### 2.1 The Potential Energy Surfaces

We used for the scattering of CD4 the same model PESs as we did for CH4. Since we expressed the PES in mass-weighted coordinated the parameters in the PESs for CD4 differs from CH4. We will now give an overview of our model PESs and the corresponding parameters for CD4. The parameters of CH4 for these PESs were already given in Ref. , where also more detailed information about our assumptions and contour plots of some cross-section of the model PESs can be found.

The PESs we used can all be written as

$$V_{\mathrm{total}}=V_{\mathrm{intra}}+V_{\mathrm{surf}},$$
(1)
where $`V_{\mathrm{intra}}`$ is the intramolecular PES and $`V_{\mathrm{surf}}`$ is the repulsive interaction with the surface. For the $`V_{\mathrm{intra}}`$ we looked at four different types of PESs. The $`V_{\mathrm{intra}}`$ include also for two types changes in the intramolecular potential due to interactions with the surface.

#### 2.1.1 A harmonic potential

The first one is completely harmonic. We have used normal mode coordinates for the internal vibrations, because these are coupled only very weakly. In the harmonic approximation this coupling is even absent so that we can write $`V_{\mathrm{intra}}`$ as

$$V_{\mathrm{intra}}=V_{\mathrm{harm}}=\frac{1}{2}\underset{i=2}{\overset{10}{}}k_iX_i^2,$$
(2)
the summation is over the internal vibrations, $`X_i`$’s are mass-weighted displacement coordinates and $`k_i`$ are mass-weighted force constants. (see Table 1 for definitions and values); ($`X_1`$ is the mass-weighted overall translation along the surface normal). The force constants have been obtained by fitting them on the experimental vibrational frequencies of CH4 and CD4.

We have assumed that the repulsive interaction with the surface is only through the deuterium atoms that point towards the surface. We take the $`z`$-axis as the surface normal. In this case the surface PES is given by

$$V_{\mathrm{surf}}=\frac{A}{N_D}\underset{i=1}{\overset{N_D}{}}e^{\alpha z_i},$$
(3)
where $`N_D`$ is the number of deuteriums that points towards the surface, $`\alpha `$=1.0726 atomic units and $`A`$=6.4127 Hartree. These parameters are chosen to give the same repulsion as the PES that has been used in an MCTDH wavepacket simulation of CH4 dissociation.

If we write $`V_{\mathrm{surf}}`$ in terms of normal mode coordinates, then we obtain for one deuterium pointing towards the surface

$$V_{\mathrm{surf}}=Ae^{\alpha _1X_1}e^{\alpha _2X_2}e^{\alpha _3X_3}e^{\alpha _4X_4},$$
(4)
where $`A`$ as above, and $`\alpha `$’s as given in Table 2. $`X_2`$, $`X_3`$ and $`X_4`$ correspond all to $`a_1`$ modes of the C3v symmetry (see Fig. 1). There is no coupling between the modes $`X_5`$ to $`X_{10}`$ in the $`V_{\mathrm{surf}}`$ part of the PES, which are all $`e`$ modes of the C3v symmetry.

For two deuteriums we obtain

$`V_{\mathrm{surf}}=A`$ $`e^{\alpha _1X_1}e^{\alpha _2X_2}e^{\alpha _3X_3}e^{\alpha _4X_4}e^{\alpha _5X_5}`$ (5)
$`\times {\displaystyle \frac{1}{2}}[`$ $`e^{\beta _3X_7}e^{\beta _3X_8}e^{\beta _5X_9}e^{\beta _5X_{10}}`$
$`+`$ $`e^{\beta _3X_7}e^{\beta _3X_8}e^{\beta _5X_9}e^{\beta _5X_{10}}],`$
with $`A`$ again as above, $`\alpha `$’s and $`\beta `$’s as given in Table 2. $`X_2`$, $`X_3`$, $`X_4`$ and $`X_5`$ correspond all to $`a_1`$ modes of C2v. $`X_7`$, $`X_8`$, $`X_9`$ and $`X_{10}`$ correspond to $`b_1`$ and $`b_2`$ modes of C2v. $`X_6`$ corresponds to the $`a_2`$ mode of C2v and has no coupling with the other modes in $`V_{\mathrm{surf}}`$.

For three deuteriums we obtain

$`V_{\mathrm{surf}}=A`$ $`e^{\alpha _1X_1}e^{\alpha _2X_2}e^{\alpha _3X_3}e^{\alpha _4X_4}`$ (6)
$`\times {\displaystyle \frac{1}{3}}[`$ $`e^{\beta _1X_5}e^{\beta _2X_6}e^{\beta _3X_7}e^{\beta _4X_8}e^{\beta _5X_9}e^{\beta _6X_{10}}`$
$`+`$ $`e^{\beta _1X_5}e^{\beta _2X_6}e^{\beta _3X_7}e^{\beta _4X_8}e^{\beta _5X_9}e^{\beta _6X_{10}}`$
$`+`$ $`e^{2\beta _1X_5}e^{2\beta _3X_7}e^{2\beta _5X_9}],`$
with $`A`$ again as above, $`\alpha `$’s and $`\beta `$’s as given in Table 2. $`X_2`$, $`X_3`$ and $`X_4`$ corresponds to $`a_1`$ modes in the C3v symmetry (see Fig. 1). Because these last six coordinates correspond to degenerate $`e`$ modes of the C3v symmetry, the $`\beta `$ parameters are not unique.

#### 2.1.2 An anharmonic intramolecular potential

Even though we do not try to describe the dissociation of methane in this and our previous paper, we do want to determine which internal vibration might be important for this dissociation. The PES should at least allow the molecule to partially distort as when dissociating. The harmonic PES does not do this. A number of changes have therefor been made. The first is that we have describe the C–D bond by a Morse PES.

$$V_{\mathrm{Morse}}=D_e\underset{i=1}{\overset{4}{}}\left[1e^{\gamma \mathrm{\Delta }r_i}\right]^2,$$
(7)
where $`D_e=0.1828`$ Hartree (the dissociation energy of methane in the gas-phase) and $`\mathrm{\Delta }r_i`$ the change in bond length from the equilibrium distance. $`\gamma `$ was calculated by equating the second derivatives along one bond of the harmonic and the Morse PES. If we transform Eq. (7) back into normal mode coordinates, we obtain

$$V_{\mathrm{Morse}}=D_e\underset{i=1}{\overset{4}{}}\left[1e^{\gamma _{i2}X_2}e^{\gamma _{i3}X_3}e^{\gamma _{i4}X_4}e^{\gamma _{i7}X_7}e^{\gamma _{i8}X_8}e^{\gamma _{i9}X_9}e^{\gamma _{i,10}X_{10}}\right]^2,$$
(8)
with $`D_e`$ as above. $`\gamma `$’s are given in Tables 3 and 4. Note that, although we have only changed the PES of the bond lengths, the $`\nu _4`$ umbrella modes are also affected. This is because these modes are not only bending, but also contain some changes of bond length.

The new intramolecular PES now becomes

$$V_{\mathrm{intra}}=V_{\mathrm{harm}}+V_{\mathrm{Morse}}V_{\mathrm{corr}},$$
(9)
where $`V_{\mathrm{harm}}`$ is as given in Eq. (2) and $`V_{\mathrm{corr}}`$ is the quadratic part of $`V_{\mathrm{Morse}}`$, which is already in $`V_{\mathrm{harm}}`$.

#### 2.1.3 Intramolecular potential with weakening C–D bonds

When the methane molecule approach the surface the overlap of substrate orbitals and anti-bonding orbitals of the molecule weakens the C–D bonds. We want to include this effect for the C–D bonds of the deuteriums pointing towards the surface. We have redefined the $`V_{\mathrm{Morse}}`$ given in Eq. (8) and also replace it in Eq. (9). A sigmoidal function is used to switch from the gas phase C–D bond to a bond close to the surface. We have used the following, somewhat arbitrary, approximations. (i) The point of inflection should be at a reasonable distance from the surface. It is set to the turnaround point for a rigid methane molecule with translation energy 93.2 kJ/mol plus twice the fall-off distance of the interaction with the surface. (ii) The depth of the PES of the C–D bond is 480 kJ/mol in the gas phase, but only 93.2 kJ/mol near the surface. The value 93.2 kJ/mol corresponds to the height of the activation barrier used in our dissociation. (iii) The exponential factor is the same as for the interaction with the surface.

If we transform to normal-mode coordinates for the particular orientations, we then obtain

$$V_{\mathrm{weak}}=D_e\underset{i=1}{\overset{4}{}}W_i\left[1e^{\gamma _{i2}X_2}e^{\gamma _{i3}X_3}e^{\gamma _{i4}X_4}e^{\gamma _{i7}X_7}e^{\gamma _{i8}X_8}e^{\gamma _{i9}X_9}e^{\gamma _{i,10}X_{10}}\right]^2,$$
(10)
where $`W_i=1`$ for non-interacting bonds and

$$W_i=\frac{1+\mathrm{\Omega }e^{\alpha _1X_1+\omega }}{1+e^{\alpha _1X_1+\omega }}$$
(11)
for the interacting bonds pointing towards the surface. $`\alpha _1`$ is as given in Table 2, $`\gamma `$’s are given in Tables 3 and 4, $`\mathrm{\Omega }=1.94210^1`$ and $`\omega =7.197`$.

#### 2.1.4 Intramolecular potential with elongation of the C–D bonds

A weakened bond generally has not only a reduced bond strength, but also an increased bond length. We include the effect of the elongation of the C–D bond length of the deuteriums that point towards the surface due to interactions with the surface. We have redefined the $`V_{\mathrm{Morse}}`$ given in Eq. (8) and also replace it in Eq. (9) for this type of PES. We have used the following approximations: (i) The transition state, as determined by Ref.  and , has a C–H bond that is 0.54 Å longer than normal. This elongation should occur at the turn around point for a rigid methane molecule with a translation energy of 93.2 kJ/mol. (ii) The exponential factor is again the same as for the interaction with the surface.

If we transform to normal-mode coordinates for the particular orientations, then we obtain

$$V_{\mathrm{shift}}=D_e\underset{i=1}{\overset{4}{}}\left[1e^{\gamma _{i2}X_2}e^{\gamma _{i3}X_3}e^{\gamma _{i4}X_4}e^{\gamma _{i7}X_7}e^{\gamma _{i8}X_8}e^{\gamma _{i9}X_9}e^{\gamma _{i,10}X_{10}}\mathrm{exp}[S_ie^{\alpha _1X_1}]\right]^2,$$
(12)
where $`\alpha _1`$ is as given in Table 2, $`\gamma `$’s are given in Tables 3 and 4. For orientation with one deuterium towards the surface we obtain; $`S_1=2.94210^2`$ and $`S_2=S_3=S_4=0`$, with two deuteriums; $`S_1=S_2=0`$ and $`S_3=S_4=1.69810^2`$, and with three deuteriums; $`S_1=0`$ and $`S_2=S_3=S_4=2.94210^2`$.

### 2.2 Initial States

The exact wave-function of a $`D`$-dimensional system, is expressed in the MCTDH approximation by the form

$$\mathrm{\Psi }_{\mathrm{MCTDH}}(q_1,\mathrm{},q_D;t)=\underset{n_1\mathrm{}n_D}{}c_{n_1\mathrm{}n_D}(t)\psi _{n_1}^{(1)}(q_1;t)\mathrm{}\psi _{n_D}^{(D)}(q_D;t).$$
(13)
All initial states in the simulations start with the vibrational ground state. The initial translational part $`\psi ^{(\mathrm{tr})}`$ is represented by a Gaussian wave-packet,

$$\psi ^{(\mathrm{tr})}(X_1)=(2\pi \sigma ^2)^{1/4}\mathrm{exp}\left[\frac{(X_1X_0)^2}{4\sigma ^2}+iP_1X_1\right],$$
(14)
where $`\sigma `$ is the width of the wave-packet (we used $`\sigma =320.248`$ atomic units), $`X_0`$ is the initial position (we used $`X_0=11\sigma `$, which is far enough from the surface to observe no repulsion) and $`P_1`$ is the initial momentum. Since we used mass-weighted coordinates the Gaussian wavepacket are identical for CD4 and CH4. We performed simulations in the energy range of 32 to 96 kJ/mol. We here present only the results of 96 kJ/mol (equivalent to $`P_1=0.2704`$ atomic units), because they showed the most obvious excitation probabilities for $`V_{\mathrm{Morse}}`$. We used seven natural single-particle states, 512 grid points and a grid-length of $`15\sigma `$ for the translational coordinate. With this grid-width we can perform simulation with a translational energy up to 144 kJ/mol.

Gauss-Hermite discrete-variable representations (DVR) were used to represent the wavepackets of the vibrational modes. We used for all simulations of CD4 the same number of DVR points as for CH4, which was 5 DVR points for the $`\nu _2`$ bending modes and 8 DVR points for the $`\nu _4`$ umbrella, $`\nu _3`$ asymmetrical stretch, and $`\nu _1`$ symmetrical stretch mode for an numerical exact integration, except for the simulations with $`V_{\mathrm{shift}}`$, where we used 16 DVR points for the $`\nu _1`$ symmetrical stretch mode, because of the change in the equilibrium position.

Also the same configurational basis was used for both isotopes. We did the simulation with one bond pointing towards the surface in eight dimensions, because the $`\nu _2`$ bending modes $`X_5`$ and $`X_6`$ do not couple with the other modes. We needed four natural single-particle states for modes $`X_2`$, $`X_3`$ and $`X_4`$, and just one for the others. So the number of configurations was $`7^14^31^4=448`$. The simulation with two bonds pointing towards the surface was performed in nine dimensions. One of the $`\nu _2`$ bending mode ($`X_6`$) does not couple with the other modes, but for the other mode $`X_5`$ we needed four natural single-particle states. The number of configurations was $`7^14^41^4=1792`$, because we needed the same number of natural single-particle states as mentioned above for the other modes. We needed ten dimensions to perform the simulation with three bonds pointing towards the surface. We used here one natural single-particle state for the modes $`X_5`$ to $`X_{10}`$ and four natural single-particle states for $`X_2`$ to $`X_4`$, which gave us $`7^14^31^6=448`$ configurations.

## 3 Results and Discussion

### 3.1 Excitation probabilities and structure deformation of CD4

The scattering probabilities for CD4 are predominantly elastic, as we also found in our previous simulations of CH4 scattering. The elastic scattering probability is larger than 0.99 for all orientation of the PESs with $`V_{\mathrm{Morse}}`$ and $`V_{\mathrm{weak}}`$ at a translational energy of 96 kJ/mol. For the PES with $`V_{\mathrm{shift}}`$ we observe an elastic scattering probability of 0.981 for the orientation with one, 0.955 with two and 0.892 with three deuteriums pointing towards the surface. This is lower than we have found for CH4, which is 0.956 for the orientation with three hydrogens pointing towards the surface and larger than 0.99 for the others. The higher inelastic scattering probabilities of CD4 was expected, because the force constants $`k_i`$ of CD4 are decreased up to 50% with respect to those of CH4 and the translational surface repulsion fall-off differs only little.

When we look at the excitation probabilities at the surface for the PES with $`V_{\mathrm{Morse}}`$ and $`V_{\mathrm{weak}}`$, then we observe generally an increase for CD4 compared with CH4, except for the $`\nu _4`$ umbrella mode in the orientation with two bond pointing towards the surface. Relevant differences in the structure deformations are observed only in the bond angles, which are increased for CD4 in the orientations with one and three bonds pointing towards the surface. The bond angle deformation of the angle between the bonds pointing towards the surface in the orientation with two bonds pointing towards the surface is decreased for CD4. We observe again that the PES with $`V_{\mathrm{weak}}`$ gives larger structure deformations than the PES with $`V_{\mathrm{Morse}}`$, but the differences are smaller for CD4 than CH4.

For the PES with $`V_{\mathrm{shift}}`$ we do not observe this effect on the bond angle deformation. The bond angle deformation for the orientation with two and three deuteriums pointing towards the surface is the same as for CH4 and it is just $`0.1^{}`$ less for the bond angle at the surface side in the orientation with one deuterium pointing towards the surface. The excitation probabilities (see Table 5) for the $`\nu _2`$ bending and $`\nu _4`$ umbrella modes become higher for all orientations for CD4, which is necessary for getting the same bond angle deformations as CH4.

The changes in the bond distances for the orientations with one and two bonds pointing towards the surface is for CD4 almost the same as for CH4. For the orientation with three bonds pointing towards the surface, we found that the maximum bond lengthening of the bonds on the surface side was $`0.032`$ Å less for CD4 than CH4. We also found that the bond shortening of the bond pointing away from the surface is $`0.010`$ Å more for CD4. These are only minimal differences, which also only suggest that the bond deformation for CD4 has been influenced slightly more by the $`\nu _3`$ asymmetrical stretch mode than the $`\nu _1`$ symmetrical stretch mode. The observed excitation probabilities for these modes do not contradict this, but are not reliable enough for hard conclusions because of their high magnitude. It is also not clear, beside of this problem, what they really represent. Is the excitation caused by a different equilibrium position of the PES at the surface in a mode or is it caused by extra energy in this mode? To answer these questions we decided to do an energy distribution analysis during the scattering for both isotopes.

### 3.2 Energy distribution in CH4 and CD4

The energy distribution analysis is performed by calculating the expectation values of the important term of the Hamiltonian $`H`$ for the wave-function $`\mathrm{\Psi }(t)`$ at a certain time $`t`$ during the scattering of CD4 and CH4 for all presented orientations in this and our previous paper. We will present here only the results of the PES with $`V_{\mathrm{shift}}`$, because it is the only model PES for which the energy distribution analysis is relevant for the discussion of the dissociation hypotheses later on.

We can obtain good information about the energy distribution per mode by looking at the kinetic energy expectation values $`\mathrm{\Psi }(t)|T_j|\mathrm{\Psi }(t)`$ per mode $`j`$ (see Table 6), because the kinetic energy operators $`T_j`$ have no cross terms like the PESs have. When we discuss the kinetic energy of a mode we normally refer to the $`a_1`$ mode of the C3v or C2v symmetry, because in these modes we have observed the highest excitation probabilities and the change in kinetic energy in the other modes is generally small.

By looking at the expectation values of some terms of the PES $`\mathrm{\Psi }(t)|V_{\mathrm{term}}|\mathrm{\Psi }(t)`$ (see Table 7), we obtain information about how the kinetics of the scattering is driven by the PES. The $`V_{\mathrm{surf}}`$ PES \[see Eqs. (4), (5) and (6)\] is the surface hydrogen/deuterium repulsion for a given orientation. $`V_{\mathrm{harm}}(\nu _2)`$ and $`V_{\mathrm{harm}}(\nu _4)`$ \[see Eq. (2)\] are the pure harmonic terms of the intramolecular PES of the $`a_1`$ modes in the C3v and C2v symmetry corresponding to a $`\nu _2`$ bending and $`\nu _4`$ umbrella modes, respectively. The pure harmonic correction terms of $`V_{\mathrm{corr}}`$ \[see Eq. (9)\] are included in them. $`V_{\mathrm{bond}}(R_{\mathrm{up}})`$ and $`V_{\mathrm{bond}}(R_{\mathrm{down}})`$ are the potential energy in a single C–H or C–D bond pointing respectively towards and away from the surface, and they give the expectation value of one bond term of $`V_{\mathrm{shift}}`$ \[see Eq. (12)\]. All given expectation values are the maximum deviation of the initial values, which effectively means the values at the moment the molecule hits the surface.

The largest changes in expectation values are, of course, in the kinetic energy of the translational mode. The translational kinetic energy does not become zero as we should expect in classical dynamics. The loss of translational kinetic energy is primary absorbed by the $`V_{\mathrm{surf}}`$ terms of the PESs. The expectation values of the $`V_{\mathrm{surf}}`$ terms show the ability of the hydrogens or deuteriums to come close to the metal surface, since in real space their exponential fall-offs are the same for both isotopes. For a rigid molecule the sum of the translational kinetic energy and $`V_{\mathrm{surf}}`$ should be constant, so all deviations of this sum have to be found back in the intramolecular kinetic energy and other PES terms.

We observe that both the minimum in the translational kinetic energy and the maximum in the $`V_{\mathrm{surf}}`$ terms were higher for CH4 than CD4, so we have to find more increase in energy in the intramolecular modes and PES terms for CD4 than CH4. We indeed do so and that can be one of the reasons we found higher inelastic scatter probabilities for CD4 for the PES with $`V_{\mathrm{shift}}`$ .

For the orientations with one and two bonds pointing towards the surface we observe a large increase of the kinetic energy in the $`\nu _3`$ asymmetrical stretch mode. If we compare this with the excitation probabilities, we find that the kinetic energy analysis gives indeed a different view on the dynamics. For the orientation with two bond pointing towards the surface we have found for both isotopes very high excitation probabilities in the $`\nu _1`$ and $`\nu _3`$ stretch modes. We know now from the kinetic energy distribution that for the $`\nu _1`$ symmetrical stretch mode the high excitation probability is caused by the change of the equilibrium position of the $`\nu _1`$ mode in the PES and that for the $`\nu _3`$ stretch mode probably the PES also has become narrower.

For the orientation with three bonds pointing towards the surface we also obtain an large increase of the kinetic energy of the $`\nu _3`$ asymmetrical stretch mode, but we also find an even larger increase in the kinetic energy of the $`\nu _1`$ symmetrical stretch mode. The total kinetic energy was extremely large, because the kinetic energy of the translational mode becomes also much larger than for the other orientations. Because of this the $`V_{\mathrm{surf}}`$ terms had to be around twice as low as for the other orientations.

All $`V_{\mathrm{bond}}(R_{\mathrm{up}})`$ terms become lower compared to the initial value, especially in the orientation with two bond pointing towards the surface. In the orientation with one bond pointing towards the surface, the $`V_{\mathrm{bond}}(R_{\mathrm{down}})`$ term became higher. This is caused by the repulsion of $`V_{\mathrm{surf}}`$ in the direction of the bond. The increase of this PES term value is higher for CD4 than CH4.

In the orientation with three bond pointing towards the surface we also observe a higher $`V_{\mathrm{bond}}(R_{\mathrm{down}})`$ value, with also the highest increase for CD4. In relation with the somewhat shorter bond distance for the $`R_{\mathrm{down}}`$ of CD4 compared with CH4, which was also a bit lower compared with the other orientations, we know now that the hydrogens and especially the deuterium have problems in following the minimum energy path of the PES with $`V_{\mathrm{shift}}`$ during the scattering dynamics. This leads to higher kinetic energy in the vibrational modes, which results in more inelastic scattering.

The $`V_{\mathrm{harm}}(\nu _2)`$ term increases in respect to the initial value, but not as much as the increase of the $`V_{\mathrm{harm}}(\nu _4)`$ is for the orientation with two bonds pointing towards the surface. The values of $`V_{\mathrm{harm}}(\nu _4)`$ for CD4 are even higher than for CH4. We observe also a larger increase of the kinetic energy in the $`\nu _4`$ umbrella mode for CD4 than for CH4. So although there is somewhat more energy transfer to the vibrational modes for CD4 than CH4, this extra vibrational energy is absorbed especially in the $`\nu _4`$ umbrella mode of CD4.

### 3.3 Dissociation hypotheses

We like now to discuss some possible implications of the scattering simulation for the isotope effect on the dissociation of methane. In our previous paper we have already drawn some conclusions about the possible reaction mechanism and which potential type would be necessary for dissociation. We found the direct breaking of a single C–H bond in the initial collision more reasonable than the splats model with single bond breaking after an intermediary Ni–C bond formation as suggested by Ref. , because the bond angle deformations seems to small to allow a Ni–C to form. From the simulations with CD4 we can draw the same conclusions. The PES with $`V_{\mathrm{shift}}`$ gives the same angle deformations for both isotopes, which is not sufficient for the splats model. The other potentials give higher bond angle deformations for the orientation with three deuteriums pointing towards the surface. If the Ni–C bond formation would go along this reaction path, then the dissociation of CD4 should be even more preferable than CH4, which is not the case. So we only have to discuss the implication of the scattering simulation for the dissociation probabilities of CH4 and CD4 for a direct breaking of a single bond reaction mechanism. This reaction mechanism can be influenced by what we will call a direct or an indirect effect.

A direct effects is the expected changes in the dissociation probability between CH4 and CD4 for a given orientation. Since we expect that we need for dissociation a PES with an elongation of the bonds pointing towards the surface, we only have to look at the isotope effect in the simulation for the PES with $`V_{\mathrm{shift}}`$ for different orientations to discuss some direct effect. It is clear from our simulations that the bond lengthening of CD4 is smaller than CH4 for the orientation with three bonds pointing towards the surface. If this orientation has a high contribution to the dissociation of methane, then this can be the reason of the higher dissociation probability of CH4. In this case our simulations also explain why Ref. did not observe a high enough isotope effect in the dissociation probability of their simulation with CH4 and CD4 modelled by a diatomic, because we do not observe a change in bond lengthening between the isotopes for the orientation with one bond pointing towards the surface.

The orientation with three bonds pointing towards the surface is also the orientation with the highest increase of the total vibrational kinetic energy for the PES with $`V_{\mathrm{shift}}`$, because the energy distribution analysis shows besides an high increase of the kinetic energy in the $`\nu _3`$ asymmetrical stretch mode also an high increase in the $`\nu _1`$ symmetrical stretch mode. Since vibrational kinetic energy can be used effectively to overcome the dissociation barrier, the orientation with three bonds indicates to be a more preferable orientation for dissociation. Moreover the relative difference in kinetic energy between both isotopes is for the $`\nu _1`$ stretch mode larger than for the $`\nu _3`$ stretch mode. If the kinetic energy in the $`\nu _1`$ stretch mode contributes significantly to overcoming the dissociation barrier, then it is another explanation for the low isotope effect in Ref..

An indirect effect is the expected changes in the dissociation probability between CH4 and CD4 through changes in the orientations distribution caused by the isotope effect in the vibrational modes. This can be the case if the favourable orientation for dissociation is not near the orientation with three bonds pointing towards the surface, but more in a region where one or two bonds pointing towards the surface. These orientations do not show a large difference in deformation for the PES with $`V_{\mathrm{shift}}`$. We can not draw immediate conclusion about the indirect effect from our simulations, since we did not include rotational motion, but our simulation show that an indirect isotope effect can exist. For the PES with $`V_{\mathrm{Morse}}`$ in the orientation with three bonds pointing towards the surface, we observe that CD4 is able to come closer to the surface than CH4. So this rotational orientation should be more preferable for CD4 than for CH4. On the other hand, if the PES is in this orientation more like $`V_{\mathrm{shift}}`$ the dissociation probability in other orientation can be decreased for CD4 through higher probability in inelastic scattering channels.

So for both effects the behaviour of the orientation with three bonds pointing towards the surface seems to be essential for a reasonable description of the dissociation mechanism of methane. A wavepacket simulation of methane scattering including one or more rotational degrees of freedom and the vibrational stretch modes will be a good starting model to study the direct and indirect effects, since most of the kinetic energy changes are observed in the stretch modes and so the bending and umbrella modes are only relevant with accurate PESs. Eventually dissociation paths can be introduced in the PES along one or more bonds.

Beside of our descriptions of the possible isotopes effect for the dissociation extracted of the scatter simulations we have to keep in mind that also a tunneling mechanism can be highly responsible for the higher observed isotope effect in the experiment and that a different dissociation barrier in the simulations can enhance this effect of tunneling.

## 4 Conclusions

The scattering is in all cases predominantly elastic. However, the observed inelastic scattering is higher for CD4 compared with previous simulation on CH4 for the PES with an elongated equilibrium bond length close to the surface. When the molecule hits the surface, we observe in general a higher vibrational excitation for CD4 than CH4. The PES with an elongated equilibrium bond length close to the surface gives for both isotopes almost the same deformations, although we observe a somewhat smaller bond lengthening for CD4 in the orientation with three bonds pointing towards the surface. The other model PESs show differences in the bond angle deformations and in the distribution of the excitation probabilities of CD4 and CH4, especially for the PES with only an anharmonic intramolecular potential.

Energy distribution analysis contributes new information on the scattering dynamics. A higher transfer of translational energy towards vibrational kinetic energy at the surface results in higher inelastic scattering. The highest increase of vibrational kinetic energy is found in the $`\nu _3`$ asymmetrical stretch modes for all orientations and also in the $`\nu _1`$ symmetrical stretch mode for the orientation with three bonds pointing towards the surface, when the PES has an elongated equilibrium bond length close to the surface.

Our simulations give an indication that the isotope effect in the methane dissociation is caused mostly by the difference in the scattering behaviour of the molecule in the orientation with three bonds pointing towards the surface. At least multiple vibrational stretch modes should be included for a reasonable description of isotope effect in the methane dissociation reaction.

## Acknowledgments

This research has been financially supported by the Council for Chemical Sciences of the Netherlands Organization for Scientific Research (CW-NWO). This work has been performed under the auspices of NIOK, the Netherlands Institute for Catalysis Research, Lab Report No. TUE-99-5-02.


format: text

text:
# Macroscopic material properties from quasi-static, microscopic simulations of a two-dimensional shear-cell

## 1 Introduction

Macroscopic continuum equations for the description of the behavior of granular media rely on constitutive equations for stress, strain, and other physical quantities describing the state of the system. One possible way of obtaining an observable like the stress is to perform discrete particle simulations herrmann98 ; cundall79  and to average over the “microscopic” quantities in the simulation, in order to obtain an averaged macroscopic quantity. In the literature, slightly different definitions for stress and strain averaging procedures can be found drescher72 ; rothenburg81 ; savage81 ; cundall82 ; goddard86 ; bathurst88 ; kruyt96 ; liao97b ; goldhirsch99 ; kuhl99 .

The outline of this study is as follows. In section 2 the discrete element simulation method is discussed and, in section 3 some averaging methods are introduced and applied to a scalar quantity, namely the volume fraction. Section 4 contains the definitions and averaging strategies for fabric, stress, and elastic strain and in section 5 some material properties are extracted from the results obtained in section 4.

## 2 Modelling discrete particles

The elementary units of granular materials are mesoscopic grains. In order to account for the excluded volume, one can assume that the grains are impenetrable but deform under stress. Since the realistic modelling of the deformations of the particles in the framework of a continuum theory walton93 ; lian96  would be much too complicated, we relate the interaction force to the overlap $`\delta `$ of two particles. Note that the evaluation of the inter particle forces based on the overlap may not be sufficient to account for the nonlinear stress distribution inside the particles. Consequently, our results presented below are of the same quality as this basic assumption.

The force laws used are material dependent, involving properties such as Young’s modulus of elasticity, and have to be validated by comparison with experimental measurements foerster94 ; labous97 ; falcon98 .

When all forces $`\stackrel{}{f}_i`$ acting on the particle $`i`$, either from other particles, from boundaries or from external forces, are known, the problem is reduced to the integration of Newton’s equations of motion for the translational and the rotational degrees of freedom

$$m_i\frac{\mathrm{d}^2}{\mathrm{d}t^2}\stackrel{}{r}_i=\stackrel{}{f}_i,\mathrm{and}I_i\frac{\mathrm{d}^2}{\mathrm{d}t^2}\stackrel{}{\mathrm{\Phi }}_i=\stackrel{}{M}_i.$$
(1)
The mass of particle $`i`$ with diameter $`d_i`$ is $`m_i`$, and its moment of inertia is $`I_i=q_im_i(a_i)^2`$, with the radius $`a_i=d_i/2`$ and the dimensionless shape factor $`q_i`$. The vectors $`\stackrel{}{r}_i`$ and $`\stackrel{}{\mathrm{\Phi }}_i`$ give the position and the orientation angles of particle $`i`$, respectively. In our model attractive forces and the presence of other phases are neglected, we focus on “dry granular media”. Particle-particle interactions are short range and active on contact only. The total force (torque) due to other particles is thus $`\stackrel{}{f}_i=_c\stackrel{}{f}_i^c`$ ($`\stackrel{}{M}_i=_c\stackrel{}{M}_i^c`$), where the sum runs over all contacts of particle $`i`$. The torque $`\stackrel{}{M}_i^c=\stackrel{}{l}_i^c\times \stackrel{}{f}_i^c`$ is related to the force $`\stackrel{}{f}_i^c`$ via the cross product with the branch vector $`\stackrel{}{l}_i^c`$ from the particle center to the point of contact $`c`$. Eq. (1) consists of six scalar equations in three dimensions and reduces to three equations in two dimensions (2D), two for the linear and one for the rotational degree of freedom. In the following the force laws for $`\stackrel{}{f}_i^c`$ accounting for excluded volume, dissipation, and friction are introduced.

### 2.1 Force laws

The particles $`i`$ and $`j`$ interact only when they are in contact so that their overlap $`\delta =\frac{1}{2}(d_i+d_j)(\stackrel{}{r}_i\stackrel{}{r}_j)\stackrel{}{n}`$ is positive, with the unit vector $`\stackrel{}{n}=(\stackrel{}{r}_i\stackrel{}{r}_j)/|\stackrel{}{r}_i\stackrel{}{r}_j|`$ that points from $`j`$ to $`i`$. The symbol ‘$``$’ denotes the scalar product of vectors or, more generally, the order-reduction by one for each of two tensors.

The first contribution to the force acting on particle $`i`$ from $`j`$ is an elastic repulsive force

$$\stackrel{}{f}_{\mathrm{n},\mathrm{el}}=k_n\delta \stackrel{}{n},$$
(2)
where $`k_n`$ is proportional to the material’s modulus of elasticity with units $`[N/m]`$. Since we are interested in disks rather than spheres, we use a linear spring that follows Hooke’s law, whereas in the case of elastic spheres, the Hertz contact law would be more appropriate hertz82 ; landau75 .

The second contribution, a viscous dissipation, is given by the damping force in the normal direction

$$\stackrel{}{f}_{\mathrm{n},\mathrm{diss}}=\gamma _n\dot{\delta }\stackrel{}{n},$$
(3)
where $`\gamma _n`$ is a phenomenological normal viscous dissipation coefficient with units \[kg s-1\] and $`\dot{\delta }=\stackrel{}{v}_{ij}\stackrel{}{n}`$ the relative velocity in the normal direction $`\stackrel{}{v}_{ij}=\stackrel{}{v}_i\stackrel{}{v}_j`$.

The third contribution to the contact force – accounting for tangential friction – can be chosen in the simplest case, according to Coulomb, as $`\stackrel{}{f}_{\mathrm{t},\mathrm{friction}}\mu |\stackrel{}{f_n}|\stackrel{}{t}`$, where $`\mu `$ is the friction coefficient and $`\stackrel{}{t}=\dot{\stackrel{}{\xi }}/|\dot{\stackrel{}{\xi }}|`$ is the tangential unit-vector parallel to the tangential component of the relative velocity $`\dot{\stackrel{}{\xi }}=\stackrel{}{v}_{ij}(\stackrel{}{v}_{ij}\stackrel{}{n})\stackrel{}{n}`$. Because this non-smooth ansatz leads to numerical problems for small $`\dot{\stackrel{}{\xi }}`$ a regularizing viscous force $`\stackrel{}{f}_{\mathrm{t},\mathrm{viscous}}=\gamma _t\dot{\stackrel{}{\xi }}`$ is added. The two forces are combined by taking the minimum value

$$\stackrel{}{f}_\mathrm{t}=\mathrm{min}(|\gamma _t\dot{\stackrel{}{\xi }}|,|\mu \stackrel{}{f_n}|)\stackrel{}{t}.$$
(4)
The effect of a more realistic tangential force law according to Cundall and Strack cundall83  will be reported elsewhere latzel99a . Due to the boundary conditions introduced below, it is also necessary to account for the friction with the bottom

$$\stackrel{}{f}_\mathrm{b}=\mu _\mathrm{b}mg\widehat{\stackrel{}{v}},$$
(5)
with the unit vector in the direction of the particles velocity $`\widehat{\stackrel{}{v}}=\stackrel{}{v}/|\stackrel{}{v}|`$. The effect of a bottom friction was discussed also in schollmann99 . In summary, we combine the forces and use

$$\stackrel{}{f}_i=\underset{c}{}(\stackrel{}{f}_{\mathrm{n},\mathrm{el}}+\stackrel{}{f}_{\mathrm{n},\mathrm{diss}}+\stackrel{}{f}_\mathrm{t})+\stackrel{}{f}_\mathrm{b}$$
(6)
for the forces acting on particle $`i`$ at its contact $`c`$ with particle $`j`$.

### 2.2 Model system

In the simulations presented in this study, a two-dimensional Couette shear-cell is used, filled with a bidisperse packing of disks, as sketched in Fig. 1. The system undergoes slow shearing introduced by turning the inner ring. The properties of the particles, used for the force laws of subsection 2.1, are summarized in table 1.

The boundary conditions are based on an experimental realization howell99 ; veje98b ; veje99 . For more details on other simulations, see veje98b ; schollmann99 .

In the simulations different global volume fractions

$$\overline{\nu }=\frac{1}{V_{\mathrm{tot}}}\underset{p=1}{\overset{N}{}}V^p$$
(7)
of the shear-cell are examined. The sum in Eq. (7) runs over all particles $`p`$ with volume $`V^p`$ in the cell with $`V_{\mathrm{tot}}=\pi (R_o^2R_i^2)`$. In this study $`\overline{\nu }`$ is varied between $`0.8084`$ and $`0.8194`$, by varying the particle number, see table 2. For the calculation of the global volume fraction, the small particles glued to the wall are counted with half their volume only, and thus contribute with $`\overline{\nu }_{\mathrm{wall}}=0.0047`$ to $`\overline{\nu }`$.

Note that we use the phrases “volume” and “volume fraction” even if, strictly speaking, the unfamiliar terms “disk area” and “area fraction” could be used. The reasons for this choice are: (i) The methods discussed in this study are straightforwardly generalized to three dimensions and (ii) the particles are three dimensional objects with height $`h`$ anyway, so that the use of the word “volume” is justified.

### 2.3 Initial conditions and steady state

The simulations are started in a dilute state with an extended outer ring $`R_o(t=0)>R_o=0.2524`$ m, and the inner ring already rotates counterclockwise with constant angular frequency $`\mathrm{\Omega }=2\pi /T_i=0.1`$ s-1 and period $`T_i=62.83`$ s. The extended outer ring is used in order to allow for a random, dilute initial configuration. The desired density is then reached by reducing the volume. The radius of the outer ring is reduced within about two seconds to reach its desired value $`R_o`$. Afterwards, the outer ring is kept fixed and the inner ring continues to rotate until at $`t=t_{\mathrm{max}}`$ the simulation ends. If not explicitly mentioned, averages are performed after about one rotation at $`t=60`$ s (to get rid of the arbitrary initial configuration), and during about one rotation, until $`t=119`$ s.

## 3 From the micro- to a macro-description

In the previous section, the microscopic point of view was introduced, as used for the discrete element method. Particles are viewed as independent entities which interact when they come in contact. In this framework, the knowledge of the forces at each contact is sufficient to model the dynamics and the statics of the system. Tensorial quantities like the stress or the deformation gradient are not necessary for a discrete modelling. However, subject of current research is to establish a correspondence to continuum theories by computing tensorial quantities like the stress $`𝝈`$, the strain $`𝜺`$, as well as scalar material properties like, e.g., the bulk and shear moduli goddard86 ; kruyt96 ; liao97b . In the course of this process, we first discuss averaging strategies using the material density as an example.

### 3.1 Averaging strategies

Most of the measurable quantities in granular materials vary strongly on short distances. Thus, computing averages necessitates dealing with or smearing out the fluctuations. In order to suppress the fluctuations, we perform averages in both time and space. This is possible due to the chosen boundary condition. The system can run for long time in a quasi-steady state and, due to the cylindrical symmetry, points at a certain distance $`R`$ from the origin are equivalent to each other. Therefore, averages are taken over many snapshots in time with time steps $`\mathrm{\Delta }t`$ and on rings of material at a distance $`r=RR_i`$ from the inner ring. The width of the averaging rings is $`\mathrm{\Delta }r`$, so that the averaging volume of one ring is $`V_r=2\pi r\mathrm{\Delta }r`$. For the sake of simplicity (and since the procedure is not restricted to cylindrical symmetry), the averaging volume is denoted by $`V=V_r`$ in the following. The rings are numbered from $`s=0`$ to $`B1`$, with $`B=(R_oR_i)/\mathrm{\Delta }r`$, and ring $`s`$ reaches from $`r_s=r\mathrm{\Delta }r/2`$ to $`r_{s+1}=r+\mathrm{\Delta }r/2`$. The averaging over many snapshots is somehow equivalent to an ensemble average. However, we remark that different snapshots are not necessarily independent of each other as discussed in subsection 3.4. Also the duration of the simulation maybe not long enough to explore a representative part of the phase space.

The cylindrical symmetry is accounted for by a rotation of all directed quantities like vectors, depending on the cartesian position $`\stackrel{}{r}_i=(x_i,y_i)`$ of the corresponding particle $`i`$. The orientation of particle $`i`$ is $`\varphi _i=\mathrm{arctan}(y_i/x_i)`$ for $`x_i>0`$ and periodically continued for $`x_i<0`$ so that $`\varphi _i`$ can be found in the interval $`[\pi ,\pi ]`$. The vector $`\stackrel{}{n}^c`$ that corresponds to contact $`c`$ of particle $`i`$ is then rotated about the angle $`\varphi _i`$ from its cartesian orientation before being used for an average. Note that this does not correspond to a transformation into orthonormal cylindrical coordinates.

Finally, we should remark that the most drastic assumption used for our averaging procedure is the fact, that all quantities are smeared out over one particle. Since it cannot be the goal to solve for the stress field inside one particle, we assume that a measured quantity is constant inside the particle. This is almost true for the density, but not e.g. for the stress. However, since we average over all positions with similar distance from the origin, i.e. averages are performed over particles with different positions inside a ring, details of the position dependency inside the particles will be smeared out anyway. An alternative approach was recently proposed by I. Goldhirsch goldhirsch99  who smeared out the averaging quantities along the lines connecting the centers of the particles and weighed the contribution according to the fraction of this line within the averaging volume.

### 3.2 Volume fractions

The first quantity to measure is the local volume fraction

$$\nu =\nu (r)=\frac{1}{V}\underset{pV}{}w_V^pV^p$$
(8)
with the particle volume $`V^p`$. $`w_V^p`$ is the weight of the particle’s contribution to the average.

This formalism can be extended to obtain the mean value of a quantity $`Q`$ in the following way:

$$Q=Q^p=\frac{1}{V}\underset{pV}{}w_V^pV^pQ^p,$$
(9)
with the pre-averaged particle quantity

$$Q^p=\underset{c=1}{\overset{𝒞^p}{}}Q^c.$$
(10)
with the quantity $`Q^c`$ attributed to contact $`c`$ of particle $`p`$. In the following the brackets $`Q`$ denote the average over a quantity $`Q`$ of the particles in an averaging volume $`V`$.

The simplest choice for $`w_V^p`$ is to use $`w_V^p=1`$, if the center of the particle lies inside the ring, and $`w_V^p=0`$ otherwise. This method will be referred to as particle-center averaging in the following.

A more complicated way to account for those parts of the particles which lie partially inside a ring is to use only the fraction of the particle volume that is covered by the averaging volume. Since an exact calculation of the area of a disk that lies in an arbitrary ring is rather complicated, we assume that the boundaries of $`V`$ are locally straight, i.e. we cut the particle in slices, as shown in Fig. 2. This method will be referred to as slicing in the following. The error introduced by using straight cuts is well below one per-cent in all situations considered here.

The volume $`V_s^p=w_V^pV^p`$ of a particle $`p`$ which partially lies between $`r_s`$ and $`r_{s+1}`$ is calculated by subtracting the external volumes $`V_s^{}`$ and $`V_s^+`$ from the particle volume $`V^p=\pi (d/2)^2`$ so that

$`V_s^p`$ $`=`$ $`V^pV_s^{}V_s^+`$ (11)
$`=`$ $`(d/2)^2[\pi \varphi _s+\mathrm{sin}(\varphi _s)\mathrm{cos}(\varphi _s)`$
$`\varphi _{s+1}+\mathrm{sin}(\varphi _{s+1})\mathrm{cos}(\varphi _{s+1})]`$
with $`\varphi _s=\mathrm{arccos}(2(r^pr_s)/d)`$ and $`\varphi _{s+1}=\mathrm{arccos}(2(r_{s+1}r^p)/d)`$. The term $`(d/2)^2\varphi `$ is the area of the segment of the circle with angle $`2\varphi `$, and the term $`(d/2)^2\mathrm{sin}(\varphi )\mathrm{cos}(\varphi )`$ is the area of the triangle belonging to the segment. In Fig. 2 the case $`s=1`$ is highlighted, and the boundaries between $`V_s^{}`$, $`V_s^p`$, and $`V_s^+`$ are indicated as thick solid lines. The two outermost slices $`V_0^p=V_1^{}`$ and $`V_s^p=V_{s1}^+`$ have to be calculated separately.

When the two averaging strategies (particle-center and slicing) are compared for different widths $`\mathrm{\Delta }r`$ of the averaging rings, so that the number of intervals is $`B=20`$, 40, or 60. For $`\mathrm{\Delta }rd_{\mathrm{small}}`$ ($`B=20`$), the results are almost independent of the averaging procedure. For finer binning $`B_s30`$, the slicing procedure converges on a master curve with weak oscillations close to the walls which are due to the wall-induced layering of the particles. The results of the particle-center averaging strongly fluctuate for $`B_c24`$. These oscillations come from ordered layers of the particles close to the walls so that the slicing method reflects the real density distribution for fine enough binning. The particle-center method, on the other hand, leads to peaks, where the centers of the particles in a layer are situated and to much smaller densities where few particle centers are found; the particle-center density is obtained rather than the material density.

### 3.3 Representative elementary volume REV

An important question is, how does the result of an averaging procedure depend on the size of the averaging volume $`V`$. We combine time- and space averaging, i.e. we average over many snapshots and over rings of width $`\mathrm{\Delta }r`$, so that the remaining “size” of the averaging volume is the width of the rings $`\mathrm{\Delta }r`$. In Fig. 3 data for $`\nu `$ at fixed position $`R=0.12`$, 0.13, 0.14, and 0.20 m, but obtained with different width $`\mathrm{\Delta }r`$, are presented. The positions correspond to $`r/d_{\mathrm{small}}2.2`$, 3.6, 4.9, and 13, when made dimensionless with the diameter of the small particles. Both the particle-center method (open symbols) and the slicing method (solid symbols) are almost identical for $`\mathrm{\Delta }r/d_{\mathrm{small}}2`$, for the larger $`\mathrm{\Delta }r`$ the averaging volume can partially lie outside of the system. For very small $`\mathrm{\Delta }r/d_{\mathrm{small}}0.1`$ the different methods lead to strongly differing results, however, the values in the limit $`\mathrm{\Delta }r0`$ are consistent, i.e. independent of $`\mathrm{\Delta }r`$ besides statistical fluctuations. In the intermediate regime $`0.1`$ $`<\mathrm{\Delta }r/d_{\mathrm{small}}<2`$, the particle-center method strongly varies, while the slicing method shows a comparatively smooth variation.

Interestingly, all methods seem to collapse at $`\mathrm{\Delta }r_{\mathrm{REV}}`$ $`d_{\mathrm{small}}`$ (and twice this value), nearly to the size of the majority of the particles. For the examined situations, we observe that the particle-center and the slicing method lead to similar results for $`0.97\mathrm{\Delta }r_{\mathrm{REV}}/d_{\mathrm{small}}1.03`$. This indicates that the systems (and measurements of system quantities) are sensitive to a typical length scale, which is here somewhat smaller than the mean particle size. When using this special $`\mathrm{\Delta }r_{\mathrm{REV}}`$ value, one has $`B=20`$ or $`B=21`$ binning intervals. The open question of this being a typical length scale that also occurs in systems with a broader size spectrum, cannot be answered with our setup, due to the given particle-size ratio.

### 3.4 Time averaging

In order to understand the fluctuations in the system over time, and to test whether subsequent snapshots can be assumed to be independent, the volume fraction $`\nu `$ is displayed for snapshots at different averaging times in Fig. 4.

Changes in density are very weak and mostly occur in the dilated shear zone for small $`r`$. From one snapshot to the next, we frequently find, that the configuration in the outer part of the shear cell has not changed, whereas a new configuration is found in the inner part. Only after rather long times does the density change also in the outer part. Thus, simulation results in the outer part are subject to stronger fluctuations because the average is performed over less independent configurations than in the inner part.

## 4 Macroscopic tensorial quantities

In this section, the averaged, macroscopic tensorial quantities in our model system are presented. The fabric tensor describes the contact network, the stress tensor describes, in this study, the stress due to normal forces, and the deformation gradient is a measure for the corresponding elastic, reversible deformations. A more detailed analysis, where the tangential forces are also taken into account, is in progress latzel99a  and some details in that direction can be found in Refs. kruyt96 ; liao97b ; kuhl99 . However, we checked the influence of the tangential forces in our system and found their effect to be always smaller than ten per-cent.

### 4.1 Micro-mechanical fabric tensor

In assemblies of grains, the forces are transmitted from one particle to the next only at the contacts of the particles. In the general case of non-spherical particles, a packing network is characterized by the vectors connecting the centers of the particles and by the particle-contact vectors. Furthermore, the local geometry of each contact is important goddard98 ; cowin88 , see Fig. 5. In our case, with spherical particles, the situation is simpler with respect to both the spherical contact geometry and the fabric.

#### 4.1.1 Fabric tensor for one particle

One quantity that describes the contact configuration of one particle to some extent, is the second order fabric tensor goddard98 ; cowin88

$$𝐅^p=\underset{c=1}{\overset{𝒞^p}{}}\stackrel{}{n}^c\stackrel{}{n}^c,$$
(12)
where $`\stackrel{}{n}^c`$ is the unit normal vector at contact $`c`$ of particle $`p`$ with $`𝒞^p`$ contacts. The symbol $``$ denotes the dyadic product in this study. Other definitions of the fabric use the so-called branch vector $`\stackrel{}{l}^{pc}`$ from the center of particle $`p`$ to its contact $`c`$, however, the unit normal and the unit branch vector are related by $`a_p\stackrel{}{n}=\stackrel{}{l}^{pc}`$ in the case of spheres or disks, so that the definition

$$𝐅^p=\frac{1}{a_p^2}\underset{c=1}{\overset{𝒞^p}{}}\stackrel{}{l}^{pc}\stackrel{}{l}^{pc}.$$
(13)
is identical to Eq. (12).

The fabric tensor in Eq. (12) is symmetric by definition and thus consists of up to three independent scalar quantities in two dimensions. The first of them, the trace (or volumetric part) $`F_V=\mathrm{tr}(𝐅^p)=(F_{\mathrm{max}}+F_{\mathrm{min}})`$, is the number of contacts of particle $`p`$, with the major and the minor eigenvalues $`F_{\mathrm{max}}`$ and $`F_{\mathrm{min}}`$, respectively. One gets from Eqs. (12) or (13) the number of contacts of particle $`p`$

$$\mathrm{tr}(𝐅^p)=\underset{c=1}{\overset{𝒞^p}{}}\mathrm{tr}(\stackrel{}{n}^c\stackrel{}{n}^c)=𝒞^p,$$
(14)
since the scalar product of $`\stackrel{}{n}^c`$ with itself is unity by definition. The second scalar, the deviator $`F_D=F_{\mathrm{max}}F_{\mathrm{min}}`$, accounts for the anisotropy of the contact network to first order, and the third, the angle $`\varphi _F`$, gives the orientation of the “major eigenvector”, i.e. the eigenvector corresponding to $`F_{\mathrm{max}}`$, with respect to the radial direction. In other words, the contact probability distribution is proportional to the function $`[F_V+F_D\mathrm{cos}(2(\varphi \varphi _F))]`$ tsoungui98b ; dubujet98 ; schollmann99 , when averaged over many particles, an approximation which is not always reasonable mehrabadi88 .

#### 4.1.2 Averaged fabric tensor

Assuming that all particles in $`V`$ contribute to the fabric with a weight of their volume $`V^p`$ one has

$$𝐅=𝐅^p=\frac{1}{V}\underset{pV}{}w_V^pV^p\underset{c=1}{\overset{𝒞^p}{}}\stackrel{}{n}^c\stackrel{}{n}^c.$$
(15)
The values of $`F_V`$, $`F_D/F_V`$, and $`\varphi _F`$, as obtained from our simulations, are plotted in Figs. 6(a-c). The trace of the fabric tensor and thus the mean number of contacts increases with increasing distance from the inner ring, and is reduced in the vicinity of the walls. With increasing mean density, the trace of $`𝐅`$ is systematically increasing, while the deviatoric fraction seems to decrease; this means that a denser system is slightly more isotropic concerning the fabric. The major eigendirection is tilted counterclockwise by somewhat more than $`\pi /4`$ from the radial outward direction, except for the innermost layer and for the strongly fluctuating outer region.

In analogy to the trace of the fabric for a single particle, the trace of the averaged fabric is

$$\mathrm{tr}(𝐅)=\frac{1}{V}\underset{pV}{}w_V^pV^p𝒞^p=𝒞^p,$$
(16)
which, in the case of a regular, periodic contact network of almost identical particles with $`a_pa`$, reduces to the sum over all of particles in $`V`$ with the prefactor $`\nu `$ defined in Eq. (8). Now, one has a relation between the coordination number, i.e. the mean number of contacts per particle $`𝒞=𝒞(r)=𝒞^p/\nu `$, the volume fraction $`\nu `$, and the averaged fabric $`𝐅`$, as a combination of Eqs. (8) and (16):

$$\mathrm{tr}(𝐅)\nu 𝒞.$$
(17)
As a test for the averaging procedure, we plot in Fig. 7 $`\mathrm{tr}(𝐅)`$ against $`\nu 𝒞`$ and obtain all data points from all simulations collapsing close to the identity curve. For a theoretical derivation of the small (about 1 per-cent) deviation due to the polydispersity, see Ref. tsoungui99 .

### 4.2 Micro-mechanical stress tensor

The micro-mechanical stress tensor is derived in a way similar to Ref. kruyt96 . For an arbitrary volume $`V`$ with surface $`V`$, the mean stress is defined as

$$𝝈=\frac{1}{V}_V𝑑V^{}𝝈^{},$$
(18)
where $`𝝈^{}=𝝈^{}(\stackrel{}{x})`$ is the position dependent local stress inside $`V`$ and $`\stackrel{}{x}`$ is the coordinate integrated over. Note that $`𝝈^{}`$ can be a strongly fluctuating function of $`\stackrel{}{x}`$.

If the stress in the pore-space vanishes, the integral can be split into a sum over the stresses $`𝝈^p`$, pre-averaged for particles $`p`$, with the respective center of mass vectors $`\stackrel{}{x}^p`$. The mean stress is thus

$$𝝈=𝝈=\frac{1}{V}\underset{pV}{}w_V^pV^p𝝈^p=\frac{1}{V}\underset{pV}{}w_V^p_{V^p}𝑑V^{}𝝈^{},$$
(19)
which corresponds to a smearing out of $`𝝈^{}`$ over the particles, with $`𝝈^p`$, the average stress in particle $`p`$, to be derived in the following.

One can rewrite the transposed stress tensor so that

$$𝝈^\mathrm{T}=\mathrm{grad}\stackrel{}{\mathrm{x}}𝝈^\mathrm{T}=\mathrm{div}(\stackrel{}{\mathrm{x}}𝝈)\stackrel{}{\mathrm{x}}\mathrm{div}𝝈,$$
(20)
by introducing the unit tensor $`𝐈=\mathrm{grad}\stackrel{}{\mathrm{x}}`$, and by applying the series rule in the first term on the right hand side. Note that the tensors $`𝝈`$ on the right hand side are transposed with respect to the left hand side $`𝝈^\mathrm{T}`$. The transposed stress is used for the following operations for the sake of simplicity, however, using the stress directly should lead to the same result. For a more detailed treatment, see Ref. latzel99a . In static equilibrium and in the absence of body forces, the term $`\mathrm{div}𝝈`$ on the right hand side vanishes and $`𝝈`$ is symmetric.

#### 4.2.1 Mean stress for one particle

Using the definition of the transposed stress tensor in Eq. (20) the integral over one particle in Eq. (19) becomes

$`(𝝈^p)^\mathrm{T}`$ $`=`$ $`{\displaystyle \frac{1}{V^p}}{\displaystyle _{V^p}}𝑑V^{}\mathrm{div}(\stackrel{}{\mathrm{x}}𝝈^{})`$ (21)
$`=`$ $`{\displaystyle \frac{1}{V^p}}{\displaystyle _{V^p}}𝑑S(\stackrel{}{x}𝝈^{})\stackrel{}{n}`$
by application of Gauss’ theorem. In Eq. (21) $`dS`$ is the surface element of $`V^p`$ on $`V^p`$ and $`\stackrel{}{n}`$ is the outwards normal unit vector. Using the definition of the stress vector $`\stackrel{}{t}=𝝈^{}\stackrel{}{n}`$ one arrives at

$$(𝝈^p)^\mathrm{T}=\frac{1}{V^p}_{V^p}𝑑S\stackrel{}{x}\stackrel{}{t}.$$
(22)
A force $`\stackrel{}{f}^c`$ acting at a contact $`c`$ with area $`\delta s^c`$ leads to a stress vector $`\stackrel{}{t}^c=\stackrel{}{f}^c/\delta s^c`$ – even in the limit of small contact area $`\delta s^ca_p`$. Here, we apply the simplifying assumption that the force $`\stackrel{}{f}^c`$ is constant on the surface $`\delta s^c`$, i.e. we do not resolve any details at the contact. Explicitly writing the surface element $`dS`$ as a sum,

$$dS=\underset{c=1}{\overset{𝒞^p}{}}\delta (|\stackrel{}{x}\stackrel{}{x}^c|)\delta s^c,$$
(23)
leads, after integration over the delta functions, to the transposed stress

$$(𝝈^p)^\mathrm{T}=\frac{1}{V^p}\underset{c=1}{\overset{𝒞^p}{}}\stackrel{}{x}^c\stackrel{}{f}^c.$$
(24)
Transposing Eq. (24) leads to an exchange of $`\stackrel{}{x}^c`$ and $`\stackrel{}{f}^c`$, and thus to the expression for the mean stress inside particle $`p`$:

$$𝝈^p=\frac{1}{V^p}\underset{c=1}{\overset{𝒞^p}{}}\stackrel{}{f}^c\stackrel{}{x}^c,$$
(25)
which is, at first glance, dependent on the frame of reference associated with the vector $`\stackrel{}{x}^c`$.

Using vector addition, one has $`\stackrel{}{x}^c=\stackrel{}{x}^p+\stackrel{}{l}^{pc}`$, with the position vector $`\stackrel{}{x}^p`$ of particle $`p`$ and the branch vector $`\stackrel{}{l}^{pc}`$, pointing from the center of mass of particle $`p`$ to contact $`c`$. Inserting this relation in Eq. (25) leads to

$`𝝈^p`$ $`=`$ $`{\displaystyle \frac{1}{V^p}}\left({\displaystyle \underset{c=1}{\overset{𝒞^p}{}}}\stackrel{}{f}^c\right)\stackrel{}{x}^p+{\displaystyle \frac{1}{V^p}}{\displaystyle \underset{c=1}{\overset{𝒞^p}{}}}\stackrel{}{f}^c\stackrel{}{l}^{pc}`$ (26)
$`=`$ $`{\displaystyle \frac{1}{V^p}}{\displaystyle \underset{c=1}{\overset{𝒞^p}{}}}\stackrel{}{f}^c\stackrel{}{l}^{pc},`$
since the first sum vanishes in static equilibrium, where

$$\underset{c=1}{\overset{𝒞^p}{}}\stackrel{}{f}^c=\stackrel{}{0}.$$
(27)
#### 4.2.2 Averaged stress tensor

Inserting Eq. (26) in Eq. (19) gives a double sum over all particles with center inside the averaging volume $`V`$ and all their contacts

$$𝝈=𝝈^p=\frac{1}{V}\underset{pV}{}w_V^p\underset{c=1}{\overset{𝒞^p}{}}\stackrel{}{f}^c\stackrel{}{l}^{pc}.$$
(28)
In other words, stress is pre-averaged over all particles and then averaged over $`V`$.

Inserting Eq. (25) in Eq. (19) leads to the mathematically identical sum

$$𝝈=\frac{1}{V}\underset{pV}{}w_V^p\underset{c=1}{\overset{𝒞^p}{}}\stackrel{}{f}^c\stackrel{}{x}^c.$$
(29)
This expression can be transformed into a sum over all contacts $`c(p,q)V`$, where at least one of the participating particles $`p`$ and $`q`$ lies inside the volume $`V`$. The force $`\stackrel{}{f}^c=\stackrel{}{f}^{pq}=\stackrel{}{f}^{qp}`$, acting at contact $`c`$, from particle $`q`$ on $`p`$ is equal, but opposite in direction to the force acting at the same contact from particle $`p`$ on $`q`$, so that

$`𝝈`$ $`=`$ $`{\displaystyle \frac{1}{V}}{\displaystyle \underset{cV}{}}(w_V^p\stackrel{}{f}^{pq}+w_V^q\stackrel{}{f}^{qp})\stackrel{}{x}^c`$ (30)
$`=`$ $`{\displaystyle \frac{1}{V}}{\displaystyle \underset{cV}{}}(w_V^pw_V^q)\stackrel{}{f}^c\stackrel{}{x}^c,`$
since $`w_V^pw_V^q=0`$ when both particles are completely inside $`V`$. Note that every contact is visited once in Eq. (30) but twice in Eqs. (28) and (29). The expression $`cV`$ means for the particle-center averaging method that contacts contribute only if one particle is inside $`V`$ while the other is outside. In the framework of the slicing method, contacts contribute if at least one particle is cut by $`V`$ so that $`w_V^pw_V^q0`$.

Given an arbitrary averaging volume $`V`$, and all the information about the forces at all contacts from the discrete element simulations, it is obvious that – from the technical point of view – the sum in Eq. (28) can be done relatively easily, whereas the sum in Eq. (30) requires the identification of the contacts at the “surface” of $`V`$, before the possibly shorter summation can be performed. Note that both expressions Eq. (28) and Eq. (30) have been used by different authors, see e.g. kruyt96 ; chang88 ; cundall82 .

The values of the volumetric stress $`\sigma _V=\mathrm{tr}(𝝈)`$, the deviatoric fraction $`\sigma _D/\sigma _V=\mathrm{dev}(𝝈)/\mathrm{tr}(𝝈)`$, as obtained from our simulations, are plotted in Figs. 8(a-b). The volumetric stress is constant, besides fluctuations, whereas the deviatoric fraction is rather large in the shear-zone and decays with increasing distance from the inner ring, like the fabric. Also, the deviatoric fraction of the stress appears slightly reduced in magnitude for increasing global density. Note that the stress varies over almost two orders of magnitude, while the mean density is changed only weakly from $`\overline{\nu }=0.8084`$ to $`0.8194`$, see table 2. In Fig. 8c $`\varphi _\sigma `$ describes the angle which the principal axis of the stress tensor is tilted from the outward direction. The stress tensor is, on average, rotated counter clockwise by somewhat less than $`\pi /4`$ from the outward direction. It is only in the outermost part that strong fluctuations exist around the mean. We now provide evidence that the fabric and stress are not co-linear.

### 4.3 Mean strain tensor

To achieve the material properties of a granular ensemble one is interested in the stress-strain relationship of the material. One of the simplest techniques used, is the application of “Voigt’s hypothesis” which assumes that the strain is uniform and that every particle displacement conforms to the mean displacement field liao97 . Thus, the expected displacement at contact $`c`$ of particle $`p`$, relative to the force free situation, and due to the mean displacement gradient $`\mathit{ϵ}`$, is $`\mathit{ϵ}\stackrel{}{l}^{pc}`$, with $`\stackrel{}{l}^{pc}=a_p\stackrel{}{n}^c`$, and the mean contact deformation $`\delta `$. Note that the linear, symmetric strain $`𝜺=\frac{1}{2}(\mathit{ϵ}+\mathit{ϵ}^\mathrm{T})`$ is not identical to the displacement gradient, in general. In our study, we follow the approach of Liao et al. liao97b , who used $`\stackrel{}{l}^{pq}`$ instead of $`\stackrel{}{l}^{pc}`$, and assume that the actual displacement field does not coincide with, but fluctuates about the mean displacement field. The difference between the actual displacement $`\stackrel{}{\mathrm{\Delta }}^{pc}`$ and the expected displacement is

$$\stackrel{}{\chi }^{pc}=\mathit{ϵ}\stackrel{}{l}^{pc}\stackrel{}{\mathrm{\Delta }}^{pc}.$$
(31)
The actual displacement is directly related to the simulations via $`\stackrel{}{\mathrm{\Delta }}^{pc}=\delta ^c\stackrel{}{n}^c`$.

If one assumes that the mean displacement field best approximates the actual displacement, one can apply a “least square fit” to the total fluctuations

$$S=\frac{1}{V}\underset{pV}{}w_V^p\underset{c=1}{\overset{𝒞^p}{}}(\stackrel{}{\chi }^{pc})^2,$$
(32)
by minimizing $`S`$ with respect to the mean displacement gradient so that

$`{\displaystyle \frac{S}{\mathit{ϵ}}}`$ $`\stackrel{!}{=}`$ $`\mathrm{𝟎}`$
$`=`$ $`{\displaystyle \frac{2}{V}}{\displaystyle \underset{pV}{}}w_V^p{\displaystyle \underset{c=1}{\overset{𝒞^p}{}}}(\mathit{ϵ}\stackrel{}{l}^{pc}\stackrel{}{\mathrm{\Delta }}^{pc}){\displaystyle \frac{}{\mathit{ϵ}}}(\mathit{ϵ}\stackrel{}{l}^{pc}\stackrel{}{\mathrm{\Delta }}^{pc}).`$
These four equations for the four components of $`\mathit{ϵ}`$ in 2D can be transformed into a relation for the mean displacement tensor as a function of the contact displacements and the branch vectors. By assuming that $`\stackrel{}{\mathrm{\Delta }}^{pc}/\mathit{ϵ}=\stackrel{}{0}`$, the expression “$`\frac{}{\mathit{ϵ}}(\mathit{ϵ}\stackrel{}{l}^{pc})`$” becomes a dyadic product “$`\stackrel{}{l}^{pc}`$”, as can be seen by writing down the equation in index notation. Extracting $`\mathit{ϵ}`$ from the sum (what leaves $`\stackrel{}{l}^{pc}\stackrel{}{l}^{pc}`$, the core of the fabric, in the first term) and multiplying the equation with $`𝐀=𝐅^1`$, the inverse of the fabric, see Eq. (15), we find that

$$\mathit{ϵ}=\frac{\pi }{V}\left(\underset{pV}{}w_V^p\underset{c=1}{\overset{𝒞^p}{}}\stackrel{}{\mathrm{\Delta }}^{pc}\stackrel{}{l}^{pc}\right)𝐀.$$
(34)
The values of $`ϵ_V`$, $`ϵ_D/ϵ_V`$, and $`\varphi _ϵ`$, as obtained from our simulations, are plotted in Figs. 9(a-c). The elastic, volumetric deformation gradient of the granulate is localized in the shear zone and the effect is stronger in the less dense systems. Due to dilation it is easier to compress the dilute material closer to the inner ring compared to the outer part. Like the fabric and stress, the strain also becomes more isotropic with increasing mean density.

## 5 Results

At first, we compute the mean-field expectation values for $`𝝈`$ and $`\mathit{ϵ}`$, in order to get a rough estimate for the orders of magnitude of the following results. Replacing, in Eq. (28), $`f^c`$ by its mean $`\overline{f}=k_n\overline{\delta }`$, $`a_p`$ by $`a=\overline{a}`$, and $`\stackrel{}{l}^{pc}`$ by $`a\stackrel{}{n}^c`$, one gets

$$\overline{𝝈}=(k_n\overline{\delta }/\pi a)𝐅.$$
(35)
Performing some similar replacements in Eq. (34), leads to

$$\overline{\mathit{ϵ}}=(\overline{\delta }/a)𝐈,$$
(36)
or

$$\overline{\mathit{ϵ}}=(\pi /k_n)𝝈𝐀.$$
(37)
The material stiffness, $`\overline{E}`$, can be defined as the ratio of the volumetric parts of stress and strain, so that one obtains from Eq. (35) and (36)

$$\overline{E}=(k_n/2\pi )\mathrm{tr}(𝐅).$$
(38)
In Fig. 10 the rescaled stiffness of the granulate is plotted against the trace of the fabric for all simulations. Note that all data collapse almost on a line, but the mean-field value underestimates the simulation data by a few per-cent. Simulation data for different $`k_n`$ and even data from simulations with neither bottom- nor tangential friction collapse with the data for fixed $`k_n`$ and different volume fractions, shown here.

The deviatoric fraction of $`𝝈`$ is plotted in Fig. 11 against the deviatoric fraction of $`𝐅`$. Here, the data are strongly underestimated by the mean-field result in Eq. (35). A similar plot for the deviatoric fraction of the deformation gradient against $`F_D/F_V`$ shows a gathering of the data close to the identity curve, in disagreement with Eq. (36). Thus, we conclude that the deviatoric parts of stress, deformation gradient and fabric are interconnected in a more complicated way than suggested by the simple mean field estimates kruyt96 ; liao97 ; luding99b ; luding99c .

In Fig. 12 the ratio of the deviatoric parts of stress and strain is plotted against the trace of the fabric. Like the material stiffness, both quantities are proportional, only $`G`$ shows much stronger fluctuations and has a proportionality factor of about $`1/3`$. We did not use the traditional definition of the shear modulus kruyt96 , since our tensors are not co-linear as shown below.

In Fig. 13 the orientations of the tensors are plotted for some simulations and in the inner part of the shear cell. In the outer part, the deviatoric fraction is usually around 10 per-cent, i.e. so small that the orientations become noisier. Simulations A and B are skipped here, because the data of the orientations are rather noisy due to the low density which leads to intermittent behavior with strong fluctuations. We observe that all orientation angles $`\varphi `$ show the same qualitative behavior, however, the fabric is tilted more than the stress which in turn is tilted more than the deformation gradient, where the orientations are measured in counter clockwise direction from the radial outward direction.

As a final cross-check, inserting the measured values for $`𝝈`$ and $`𝐀`$ into Eq. (37) leads to the measured values of $`\mathit{ϵ}`$ within a few percent deviation. However, the orientation of the deformation gradient is not well reproduced – it seems very sensitive to small fluctuations.

## 6 Summary and Conclusion

Discrete element simulations of a 2D Couette shear cell were presented and used as the basis for a micro-macro averaging procedure. In the shear cell a shear band is localized close to the inner, rotating cylindrical wall. The boundary conditions were chosen to allow for averaging over large volumes (rings with width $`\mathrm{\Delta }r`$) and over a steady state and thus over long times. The configurations changed rather rapidly in the shear band, whereas the system is frozen in the outer part, a fact which requires either extremely long simulations or a sampling over different initial conditions in order to allow an ensemble averaging in the outer part.

The simplest averaging strategy involves only the particle-centers as carriers of the quantities to be averaged over, whereas a more advanced method assumes the quantities to be homogeneously smeared out over the whole particle which is cut in slices by the averaging volumes. Both methods agree if the averaging volumes are of the particle size (or multiples), but for other sizes differing results are obtained. The slicing method shows discretization effects in the range of averaging volume widths from one to one fifth of a particle diameter, while the particle-center method shows fluctuations due to the choice of the averaging volume in a much wider range.

The material density, i.e. the volume fraction, the coordination number, the fabric tensor, the stress tensor and the elastic, reversible deformation gradient were obtained by the averaging procedures. The fabric is linearly proportional to the product of volume fraction and coordination number. In the shear band, dilation together with a reduction of the number of contacts is observed. The mean stress is constant in radial direction while the deformation gradient decays with the distance from the inner wall. The ratio of the volumetric parts of stress and strain gives the effective stiffness of the granulate, which is small in the shear band and larger outside, due to dilation.

In the shear band, large deviatoric components of all tensorial quantities are found, however, decreasing with increasing distance from the inner wall. The isotropy of the tensors grows only slightly with increasing density and all tensors are tilted counterclockwise from the radial direction by an angle of the order of $`\pi /4`$. The system organizes itself such that more contacts are created to act against the shear. An essential observation is that the macroscopic tensors are not co-linear, i.e. their orientations are different. The orientation of the fabric is tilted most, that of the deformation gradient is tilted least and thus, the material cannot be described by a simple elastic model involving only two Lamé constants (or bulk modulus and Poisson’s ratio) as the only parameters. Alternatives are to cut the system into pieces with different material properties and thus introducing discontinuities herrmann98  or to use the rank four stiffness tensor for anisotropic materials luding99c . The deviatoric parts of stress and deformation gradient are seemingly interconnected via the fabric tensor.

To conclude, we proposed a consistent averaging formalism to obtain a mean quantity $`Q`$ in average over arbitrary volumes $`V`$. Within this framework, we used for $`Q^c`$ the fabric $`𝐅^c=\stackrel{}{n}^c\stackrel{}{n}^c`$, the stress $`𝝈^c=(1/V^p)\stackrel{}{f}^c\stackrel{}{l}^{pc}`$, and the deformation gradient $`\mathit{ϵ}^c=(\delta ^c/a_p)\stackrel{}{n}^c\stackrel{}{n}^c𝐀`$. Future work will involve the extension of the present analysis to measure also the fluctuations of the above quantities and, in addition, other interesting quantities like, e.g. , the plastic strain and non-symmetric parts of the stress tensor due to the effective moments acting on single particles. More systematic parameter studies are currently in progress.