2.1. $\unicode{x1D64C}$ -tensor theory

The orientation of a particle $i$ is defined by the dyad $({\boldsymbol{k}}{_i}\otimes {\boldsymbol{k}}{_i})$ , where ${\boldsymbol{k}}{_i}$ is a unit vector aligned with the grain symmetry axis. For an assembly of $N$ particles, orientation is defined by the traceless, symmetric tensor

(2.1) \begin{equation} {\unicode{x1D64C}} = \frac {1}{N}\sum _{i=1}^N \left ( {\boldsymbol{k}}_i \otimes {\boldsymbol{k}}_i\right )-\dfrac {1}{3}\unicode{x1D644}, \end{equation}

where $\unicode{x1D644}$ is the identity matrix. This tensor was first introduced in the context of liquid crystals (Mottram & Newton Reference Mottram and Newton2014; Liu & Wang Reference Liu and Wang2016), and has all eigenvalues in $[-1/3,2/3]$ .

The orientational tensor is assumed to be governed by a balance law of the form (Xu Reference Xu2022)

(2.2) \begin{equation} \dot {{\unicode{x1D64C}}} = \unicode{x1D64E} +\varGamma {\unicode{x1D643}}, \end{equation}

where $\dot {{\unicode{x1D64C}}}$ is the material time derivative of the orientational tensor, $\unicode{x1D64E}$ accounts for the co-rotation and stretching effect of the flow and $\unicode{x1D643}$ , proportional to $\rho _p\nu T$ , is responsible for the relaxation to the minimum of the free energy density. The coefficient $\varGamma$ , which is proportional to $1/ (\rho _p\nu T^{1/2} d_v )$ , represents rotational diffusion and is actually the inverse of a rotational viscosity (Mottram & Newton Reference Mottram and Newton2014; Liu & Wang Reference Liu and Wang2016).

The expression for $\unicode{x1D64E}$ in the liquid crystal literature (Xu Reference Xu2022), which was also formulated in Vescovi et al. (Reference Vescovi, Nadler and Berzi2024b ) for steady-state homogeneous shear flows of cylinders, reads

(2.3) \begin{equation} \unicode{x1D64E} = \left ({\boldsymbol{\varOmega }}+ \phi {\unicode{x1D63F}}\right )\left ({\unicode{x1D64C}}+\dfrac {1}{3}\unicode{x1D644}\right ) -\left ({\unicode{x1D64C}}+\dfrac {1}{3}\unicode{x1D644}\right )\left ({\boldsymbol{\varOmega }}-\phi {\unicode{x1D63F}}\right ) - 2\phi \left ( {\unicode{x1D64C}} : {\unicode{x1D63F}} \right ) \left ({\unicode{x1D64C}}+\dfrac {1}{3}\unicode{x1D644}\right )\!, \end{equation}

where ${\boldsymbol{\varOmega }}=(\boldsymbol{\nabla }{\boldsymbol{v}} - \boldsymbol{\nabla }{\boldsymbol{v}}^{\textrm{T}})/2$ is the spin tensor, ${\unicode{x1D63F}}=(\boldsymbol{\nabla }{\boldsymbol{v}} + \boldsymbol{\nabla }{\boldsymbol{v}}^{\textrm{T}})/2$ is the deformation rate of the velocity field $\boldsymbol{v}$ , and the flow alignment parameter $\phi$ was determined by fitting to discrete simulations (Vescovi et al. Reference Vescovi, Nadler and Berzi2024b ).

The usual assumption is that the free energy density, $\mathcal{F}$ , of the particles is (Mottram & Newton Reference Mottram and Newton2014)

(2.4) \begin{equation} \mathcal{F} = \dfrac {a}{2}{\textrm{tr}} ({\unicode{x1D64C}\,}^2)-\dfrac {b}{3}{\textrm{tr}} ({\unicode{x1D64C}\,}^3)+\dfrac {c}{4}{\textrm{tr}}^2\big ({\unicode{x1D64C}}^2\big )+\dfrac {L}{2}\big \lvert \boldsymbol{\nabla }{\unicode{x1D64C}}\big \rvert ^2. \end{equation}

The first three terms, which contain the coefficients $a$ , $b$ and $c$ , represent the Landau–de Gennes form for the bulk free energy density (Mottram & Newton Reference Mottram and Newton2014, a fourth-order degree polynomial). They ensure that $\mathcal{F}$ has two minima: one is , the isotropic state; and one is , the nematic state. The last term on the right-hand side of (2.4) instead represents the (entropic) elastic component of the free energy density and is introduced to penalise the distortion of the alignment. Here, $L$ is the Frank coefficient in the one-constant approximation of the Oseen–Frank (entropic) elastic energy (Heymans & Schilling Reference Heymans and Schilling2017).

The variational principle implies that ${\unicode{x1D643}}=-\partial \mathcal{F}/\partial {\unicode{x1D64C}}+\boldsymbol{\nabla }\boldsymbol{\cdot } (\partial \mathcal{F}/\partial \boldsymbol{\nabla \!}{\unicode{x1D64C}} )$ is the term responsible for driving $\unicode{x1D64C}$ to equilibrium. From (2.4),

(2.5) \begin{equation} \dfrac {{\unicode{x1D643}}}{\rho _p\nu T} = - \tilde {a} {\unicode{x1D64C}}+\tilde {b}\left [{\unicode{x1D64C}}^2-\dfrac {1}{3}{\textrm{tr}} ({\unicode{x1D64C}\,}^2)\unicode{x1D644}\right ]- \tilde {c}{\unicode{x1D64C}}\, {\textrm{tr}} ({\unicode{x1D64C}\,}^2 )+\dfrac {1}{ T}\boldsymbol{\nabla }\boldsymbol{\cdot }\big ( \tilde Ld_v^2 T \boldsymbol{\nabla }{\unicode{x1D64C}}\big ). \end{equation}

The parameters $\tilde a$ , $\tilde b$ and $\tilde c$ are the dimensionless form of the corresponding parameters $a$ , $b$ and $c$ in the bulk free energy. The parameter $\tilde L=L/ ( \rho _p\nu T d_v^2 )$ is the dimensionless form of the Frank elastic coefficient and depends on the shape factor. Unlike the common assumption of the Beris–Edwards model (Xu Reference Xu2022) that $L$ is spatially homogeneous, here we account for the possible spatial variations of at least the temperature in the expression for $L$ , and write the elastic term as the divergence of a flux of $\unicode{x1D64C}$ .

In the case of liquid crystals composed of hard particles, $\tilde b$ and $\tilde c$ are positive functions of the shape factor and independent of the solid volume fraction. On the other hand, $\tilde a=\tilde \alpha (\nu ^*-\nu )$ , where $\tilde \alpha \gt 0$ depends on the shape factor and $\nu ^*$ is the specific volume fraction at which $\tilde a$ changes sign (a necessary condition for the isotropic–nematic transition to take place (Mottram & Newton Reference Mottram and Newton2014). Given that the difference between $\nu ^*$ and the volume fraction at which the isotropic–nematic transition actually takes place is small (Mottram & Newton Reference Mottram and Newton2014), we can confuse $\nu ^*$ with the latter and fit the data obtained in Monte Carlo simulations of hard ellipsoids (Odriozola Reference Odriozola2012) with the following function of the shape factor (see Appendix A):

(2.6) \begin{equation} \nu ^* \approx -0.77r_g^2+0.04r_g+0.71. \end{equation}

2.2. DEM simulations and kinetic theory

We have performed DEM simulations, using MercuryDPM (www.mercurydpm.org, Weinhart et al. Reference Weinhart2019), of a densely packed assembly of $N$ identical frictionless ellipsoids of mass density $\rho _p$ and restitution coefficient $e_n=0.95$ (the negative of the ratio of post- to pre-collisional relative velocity normal to the plane of contact). Particles interact through linear spring-dashpots. More details about the contact force model are provided in Appendix B. The ellipsoids are contained in a box cell with four solid flat frictionless walls, in the absence of gravity. Two parallel walls perpendicular to the $y$ -direction vibrate with constant amplitude $\lambda =0.1d_v$ and constant angular frequency $\omega$ , while the other two parallel walls perpendicular to the $x$ -direction do not move. Periodic boundary conditions are used in the $z$ -direction. The dimensions of the box along the three directions are $l_x=l_y=10 d_v$ and $l_z=5 d_v$ .

We initially placed the ellipsoids at random in the box and let the system evolve. Unless otherwise stated, the average solid volume fraction is $\bar \nu = N\pi d_v^3/ (6l_xl_yl_z )=0.4$ . With this, we anticipate that the local solid volume fraction is in the range 0.28–0.58 everywhere in the box cell, well below the critical value at which we expect the presence of rate-independent components of stresses (Berzi, Buettner & Curtis Reference Berzi, Buettner and Curtis2022) and below the value at which ellipsoids such as those employed in the present simulations experience a nematic–smectic transition (Odriozola Reference Odriozola2012). In this condition, the only velocity scale of the system is $\lambda \omega$ . Hence, without loss of generality, we have set $\rho _p$ , $d_v$ and $\lambda \omega$ equal to one in the simulations.

If the system reaches a steady state with zero mean velocity, and the balance of orientation (2.2) takes the simple form

(2.7) \begin{equation} -\boldsymbol{\nabla }\boldsymbol{\cdot }\big (-\tilde {L}d_v^2T\boldsymbol{\nabla }{\unicode{x1D64C}}\big ) = T \left \{\tilde {a} {\unicode{x1D64C}}-\tilde {b}\left [{\unicode{x1D64C}}^2-\dfrac {1}{3}{\textrm{tr}} ({\unicode{x1D64C}\,}^2)\unicode{x1D644}\right ]+ \tilde {c}{\unicode{x1D64C}}\, {\textrm{tr}} ({\unicode{x1D64C}\,}^2 )\right \}\!. \end{equation}

The left-hand side of (2.7) is the negative of the divergence of the orientational flux.

Similarly, the granular temperature is governed by the balance of kinetic energy associated with the fluctuations in translational velocity (Garzó & Dufty Reference Garzó and Dufty1999), which, in the absence of shearing and in the steady state, reduces to a balance between the diffusion of energy and its dissipation via inelastic collisions:

(2.8) \begin{equation} -\!\boldsymbol{\nabla }\boldsymbol{\cdot } \big (- {\unicode{x1D646}}\, \boldsymbol{\nabla }T\big ) = \dfrac {6\big (1-e_n\big )p}{\pi ^{1/2}d_v}f\big (r_g\big )T^{1/2}, \end{equation}

where $p$ is the pressure. The left-hand side of (2.8) is the negative of the divergence of the fluctuation energy flux, in which $\unicode{x1D646}$ is the diffusivity tensor. The right-hand side of (2.8) is the rate of collisional dissipation and includes a function of the shape factor, $f (r_g )$ , that modifies the classical expression for spherical grains to account for the additional dissipation induced by particle rotation. An explicit expression,

(2.9) \begin{equation} f\big (r_g\big )=\left (0.65\dfrac {1+\lvert r_g\rvert }{1-\lvert r_g\rvert }+1.39\right ) \end{equation}

has been proposed in the case of elongated cylinders by Buettner et al. (Reference Buettner, Guo and Curtis2017). With respect to the expression proposed by Buettner et al. (Reference Buettner, Guo and Curtis2017), here we use the absolute value of the shape factor to deal with both prolate and oblate grains. Notice that we expect (2.9) to overestimate dissipation in the case of ellipsoids: cylinders have a larger excluded volume with respect to ellipsoids, leading to higher frequency of collisions and larger dissipation; also, the presence of sharp ends and flat sides should lead to more abrupt collisions in the case of cylinders, and a stronger coupling between rotational and translational motion. Also notice that modifications of the collisional dissipation rate to deal with frictional rather than frictionless particles are available (Buettner, Guo & Curtis Reference Buettner, Guo and Curtis2020).

In the case of spheres, the diffusivity is isotropic and $\unicode{x1D646}$ can be written as

(2.10) \begin{equation} {\unicode{x1D646}} = \dfrac {2M_\infty d_v p}{\pi ^{1/2}\left (1+e_n\right )T^{1/2}}\unicode{x1D644}. \end{equation}

Here, with respect to the expression derived by Garzó & Dufty (Reference Garzó and Dufty1999), we have made use of the fact that, in the dense limit, the dependence of the diffusivity on the solid volume fraction can be equivalently expressed as $p/ [2 (1+e_n )T ]$ , (Berzi Reference Berzi2024), and

(2.11) \begin{equation} M_\infty = 4.5\left [\dfrac {1+e_n}{2}+\dfrac {9\pi }{8}\dfrac {\left (1+e_n\right )^3\left (2e_n-1\right )}{16\left (1+e_n\right )-7\left (1-e_n^2\right )}\right ]\!. \end{equation}

Note that we have introduced a prefactor equal to 4.5 with respect to the usual expression for $M_\infty$ , otherwise the diffusivity would be underestimated even in the case of spheres. This underestimation was already noticed in Vescovi et al. (Reference Vescovi, de Wijn, Cross and Berzi2024a ) in the case of steady, homogeneous shearing of frictionless spheres between bumpy walls.

In the case of uniaxial grains with some preferential orientation, the anisotropy implies that the fluctuation energy can diffuse more efficiently along some directions. Similarly to what we have done for the shear viscosity in the case of shearing flows of uniaxial grains (Vescovi et al. Reference Vescovi, Nadler and Berzi2024b ), we introduce a linear dependence of the diffusivity tensor on the orientational tensor,

(2.12) \begin{equation} {\unicode{x1D646}} = \dfrac {2M_\infty d_v p}{\pi ^{1/2}\left (1+e_n\right )T^{1/2}}\big (\unicode{x1D644}+\varphi {\unicode{x1D64C}}\big ). \end{equation}

The coefficient $\varphi$ must depend on the shape of the particle and the surface friction and accounts for the anisotropy in the thermal diffusivity between the direction parallel and perpendicular to the alignment. Previous experimental and theoretical work in the context of liquid crystals showed that the diffusivity is enhanced in the direction parallel (perpendicular) to the alignment for prolate (oblate) particles, i.e. $\varphi$ is positive (negative) (Simões et al. Reference Simões, Palangana, Steudel, Kimura and Gómez2008). It also showed that $\varphi$ should be of order unity.

Assuming that the pressure $p$ in steady state is homogeneous, (2.7) and (2.8) simplify to

(2.13) \begin{equation} \delta ^2 \boldsymbol{\nabla }\boldsymbol{\cdot } (T\boldsymbol{\nabla }{\unicode{x1D64C}})= T\left \{ (\nu ^*-\nu ) {\unicode{x1D64C}}-\dfrac {\tilde {b}}{\tilde {\alpha }}\left [{\unicode{x1D64C}\,}^2-\dfrac {1}{3}{\textrm{tr}} ({\unicode{x1D64C}\,}^2)\unicode{x1D644}\right ]+ \dfrac {\tilde {c}}{\tilde {\alpha }}{\unicode{x1D64C}}\, {\textrm{tr}} ({\unicode{x1D64C}\,}^2)\right \} \end{equation}

and

(2.14) \begin{equation} \epsilon ^2 \boldsymbol{\nabla }\boldsymbol{\cdot } \big [\big (\unicode{x1D644}+\varphi {\unicode{x1D64C}}\big )\boldsymbol{\nabla }T^{1/2}\big ] = T^{1/2}, \end{equation}

respectively, where $\delta = d_v (\tilde L/\tilde \alpha )^{1/2}$ and

(2.15) \begin{equation} \epsilon = d_v \left [\dfrac {2M_\infty }{3\big (1-e_n^2\big )f (r_g)}\right ]^{1/2} \end{equation}

are two characteristic lengths for the diffusion of orientation and velocity fluctuations.

In the context of liquid crystals, the Beris–Edwards model (Beris & Edwards Reference Beris and Edwards1994; Xu Reference Xu2022) assumes incompressibility and therefore treats the solid volume fraction $\nu$ in (2.13) as a constant. Here, we will show that allowing for variations in the solid volume fraction permits the reproduction of striking features of the numerical simulations, such as the isotropic–nematic transition taking place inside the domain, at some distance from the boundaries. To account for compressibility, we must introduce an equation of state.

Explicit expressions for the equation of state in the isotropic phase of hard convex bodies have been proposed in the literature (Vega, Lago & Garzon Reference Vega, Lago and Garzon1994; Vega Reference Vega1997), without including the role of the coefficient of restitution. Recently, it has been shown that the particle shape has only a small influence on the equation of state in DEM simulations of cylinders under simple shearing (Berzi et al. Reference Berzi, Thai-Quang, Guo and Curtis2016), and that the Carnahan & Starling (Reference Carnahan and Starling1969) expression,

(2.16) \begin{equation} p = \rho _p T \dfrac {\nu +\nu ^2+\nu ^3-\nu ^4}{\left (1-\nu \right )^3}, \end{equation}

is sufficiently accurate up to moderate volume fractions, at least if the coefficient of restitution is not too far from one. To facilitate the numerical integration of the differential equations, here we expand this expression around the average volume fraction $\bar \nu$ and express $\nu$ in (2.13) as

(2.17) \begin{equation} \nu =\bar \nu + \left [\dfrac {p}{\rho _p T } - \dfrac {\bar \nu +\bar \nu ^2+\bar \nu ^3-\bar \nu ^4}{\left (1-\bar \nu \right )^3}\right ]\dfrac {\left (1-\bar \nu \right )^4}{\bar \nu ^4-4\bar \nu ^3+4\bar \nu ^2+4\bar \nu +1}. \end{equation}

It is worth mentioning that the equation of state (2.16) is not valid in the nematic phase (Frenkel & Mulder Reference Frenkel and Mulder1985), but we claim that it is still sufficiently accurate if the solid volume fraction is not too far from the value at the isotropic–nematic transition.

2.3. Boundary conditions

For the granular temperature, the balance of fluctuation energy phrased at the boundary prescribes that the component of the fluctuation energy flux normal to the wall and directed inward must be equal to $I-D$ where $I\approx (\lambda \omega )^2 p/ T_w^{1/2}$ is the energy input from a vibrating wall (Soto Reference Soto2004), and $D= (2/\pi )^{1/2} (1-e_w) p T_w^{1/2}$ is the energy dissipation from the impacts of particles on a wall (Richman Reference Richman1988; Jenkins Reference Jenkins2007), where $e_w$ is the coefficient of restitution between the wall and the particles. The subscript $w$ indicates quantities evaluated at the walls. Here we consider $e_w = e_n = 0.95$ .

With the flux of fluctuation energy given in (2.8), the boundary condition for the square root of the granular temperature at the walls is, then,

(2.18) \begin{equation} {\boldsymbol{n}}\boldsymbol{\cdot } \big [(\unicode{x1D644}+\varphi {\unicode{x1D64C}}) \boldsymbol{\nabla }T^{1/2}\big ] \big \rvert _w = \left [1 -\dfrac {\beta }{\zeta }\dfrac {(\lambda \omega )^2}{ T_w} \right ] \dfrac {T_w^{1/2}}{\beta }, \end{equation}

where $\boldsymbol{n}$ is the unit inward normal to the boundary. Here,

(2.19) \begin{equation} \zeta = d_v \dfrac {4M_\infty } {\left (2\pi \right )^{1/2}\left (1+e_n\right )} \end{equation}

and

(2.20) \begin{equation} \beta = d_v \dfrac {4M_\infty } {2^{1/2}\left (1-e_w\right )\left (1+e_n\right )} \end{equation}

are two characteristic lengths associated with the energy input from vibrating walls and the energy dissipated in wall collisions, respectively. The factor $2^{1/2}$ in the denominator of (2.19) allows us to better fit the results of discrete simulations in the case of spheres.

For the boundary condition on the orientation, we adopt the usual assumption that the normal inward flux of $\unicode{x1D64C}$ at the boundary is proportional to the jump in $\unicode{x1D64C}$ there. This is known as weak anchoring in the context of liquid crystals (Mottram & Newton Reference Mottram and Newton2014). Then, with the flux of orientation given in (2.7),

(2.21) \begin{equation} {\boldsymbol{n}}\boldsymbol{\cdot }\big ( \tilde L d_v^2T\boldsymbol{\nabla }{\unicode{x1D64C}}\big )\big \rvert _w = \dfrac {W}{\rho _p\nu _w } [\![\! {\unicode{x1D64C}} ]\!]_w , \end{equation}

where $[\![\! {\unicode{x1D64C}} ]\!]_w = {\unicode{x1D64C}}_w - {\unicode{x1D64C}}_b$ is the jump in the orientational tensor at the boundary, ${\unicode{x1D64C}}_b$ is the preferred orientation at the boundary, and $W$ is the anchoring strength (Mottram & Newton Reference Mottram and Newton2014). Here, we assume that

(2.22) \begin{equation} W = \rho _p\nu T_w\dfrac {T_w}{\left (\lambda \omega \right )^2}\tilde W d_v, \end{equation}

that is, the anchoring strength is an Onsager-like term due to loss of entropy near the wall (proportional to $\rho _p\nu _w T_w d_v$ , Poniewierski & Hoiyst Reference Poniewierski and Hoiyst1988) dynamically weakened in case of wall vibration (the correction is the ratio of $\rho _p\nu T_w$ over the average kinetic energy of particles moving with the wall, $\rho _p\nu (\lambda \omega )^2$ ). The dimensionless coefficient $\tilde W$ is a function of the particle shape. Then, (2.21) becomes

(2.23) \begin{equation} {\boldsymbol{n}}\boldsymbol{\cdot } \big (T\boldsymbol{\nabla }{\unicode{x1D64C}}\big )\big \rvert _w = \dfrac {T_w^2}{\left (\lambda \omega \right )^2}\dfrac {[\![\! {\unicode{x1D64C}} ]\!]_w}{\gamma } \end{equation}

where $\gamma = d_v (\tilde L/\tilde W )$ is another characteristic length which depends on the particle shape. It should be noted that at the non-moving walls, the Neumann boundary condition (2.23) reduces to a Dirichlet boundary condition, i.e. zero orientational jump at the boundary.

The preferred orientation ${\unicode{x1D64C}}_b$ on flat walls is different for oblate and prolate particles. Oblate particles ( $r_g \lt 0$ ) prefer to orient normal to flat walls (homeotropic anchoring, Poniewierski & Hoiyst Reference Poniewierski and Hoiyst1988). Hence, ${\unicode{x1D64C}}_b ={\boldsymbol{n}}\otimes {\boldsymbol{n}}-\unicode{x1D644}/3$ and the orientational jump is

(2.24) \begin{equation} [\![\! {\unicode{x1D64C}} ]\!]_{w,r_g\lt 0} = {\unicode{x1D64C}}_w - {\boldsymbol{n}} \otimes {\boldsymbol{n}}+\dfrac {1}{3}\unicode{x1D644} . \end{equation}

Prolate particles ( $r_g \gt 0$ ) prefer to orient tangent to flat walls, with no preferential direction there. This is known as degenerate planar anchoring (Mottram & Newton Reference Mottram and Newton2014). Here, we propose a modification of previous expressions (Fournier & Galatola Reference Fournier and Galatola2005) and set ${\unicode{x1D64C}}_b={\unicode{x1D64C}\,}^\perp -{\textrm{tr}} ({\unicode{x1D64C}\,}^\perp ){\unicode{x1D64B}}/2$ , where ${\unicode{x1D64C}\,}^\perp ={\unicode{x1D64B}} ({\unicode{x1D64C}}_w+\unicode{x1D644}/3 ){\unicode{x1D64B}}-\unicode{x1D644}/3$ is the projection of the orientational tensor onto the plane tangent to the wall and ${\unicode{x1D64B}}=\unicode{x1D644}-{\boldsymbol{n}} \otimes {\boldsymbol{n}}$ is the projection operator. This expression has the advantage with respect to that proposed by Fournier & Galatola (Reference Fournier and Galatola2005) that it allows for both uniaxial and biaxial ordering at the boundary, without constraining the order parameters, and that ${\unicode{x1D64C}}_b$ is symmetric, traceless and with all eigenvalues in $[-1/3,2/3]$ . The orientational jump, then, reads

(2.25) \begin{equation} [\![\! {\unicode{x1D64C}} ]\!]_{w,r_g\gt 0} = {\unicode{x1D64C}}_w-{\unicode{x1D64B}}{\unicode{x1D64C}}_w{\unicode{x1D64B}}+\dfrac {1}{2}{\textrm{tr}}\big ({\unicode{x1D64B}}{\unicode{x1D64C}}_w{\unicode{x1D64B}}\,\big ){\unicode{x1D64B}}-\dfrac {1}{2}{\unicode{x1D64B}}+\dfrac {1}{3}\unicode{x1D644}. \end{equation}

The two second-order differential equations (2.13) and (2.14), complemented with the boundary conditions (2.18) and (2.23), with (2.24) or (2.25), are fully determined once the two dimensionless ratios, $\tilde b/\tilde \alpha$ and $\tilde c/\tilde \alpha$ , the critical volume fraction $\nu ^*$ , the orientation-enhanced diffusivity parameter $\varphi$ , the pressure $p$ , and the five characteristic lengths, $\zeta$ , $\beta$ , $\gamma$ , $\delta$ and $\epsilon$ are given.

Both $\tilde b/\tilde \alpha$ and $\tilde c/\tilde \alpha$ depend on the particle shape and are evaluated by fitting the results of DEM simulations in a different flow configuration (simple shearing of true cylinders, as detailed in Appendix A).

As mentioned, the dependence of $\nu ^*$ on the shape factor has been determined in Monte Carlo simulations of hard ellipsoids (Odriozola Reference Odriozola2012).

We have verified that the orientation-enhanced diffusivity parameter has no discernible influence on the results. This is because one needs both spatial variations of the granular temperature and alignment to observe anisotropy in the diffusion of fluctuation kinetic energy. We will show that the gradients of $T$ are significant along the direction perpendicular to the vibrating walls, but the vibration greatly weakens the alignment, while the non-moving walls promote alignment, but act as adiabatic boundaries. Then, here, for simplicity, we set $\varphi =0$ , with the understanding that it could play a more significant role in alternative configurations.

We adjust the value of $p$ so that the spatial average of the solid volume fraction given by (2.17) in terms of the temperature field matches the value $\bar \nu$ imposed in the simulations.

Finally, the five characteristic lengths scale with the equivalent particle diameter $d_v$ , and three of them ( $\zeta$ , $\beta$ and $\epsilon$ ) are known functions of $e_n$ , $e_w$ and $r_g$ from previous studies. Only $\delta$ and $\gamma$ need to be fitted with the current DEM simulations. The list of material parameters, together with their units and their physical interpretation, is summarised in table 1.

Table 1. List of material parameters in the governing equations.

Even if this work focuses on systems composed of identical particles, we expect that poly-dispersity would only quantitatively change the results, by affecting the values of the parameters in the balance of the $\unicode{x1D64C}$ -tensor and the solid volume fraction at the isotropic–nematic transition, but not the qualitative behaviour of the system.

We use a routine implemented in Matlab to solve the system. Unless stated otherwise, the results of the DEM simulations shown in the plots are temporal means once the system reaches a steady state.