2.1. Geometry and operating conditions

The computational domain of the CFD-DEM simulations is represented in figure 1(a). The internal radius of the cylindrical column is $R_c = 0.077$ m. Above height $H_C = 12 R_c$, the radius is enlarged progressively over a height $H_E = 3.7 R_c$ until $R_E = 1.65 R_c$. The total height of the column is $H_T = 22.2 R_c$. A pressurised gas (nitrogen) is injected with uniform velocity $V_f$ on the bottom boundary (inlet). A pressure $12$ bar is prescribed on the top boundary (outlet). This elevated pressure is not necessary for the development of the present wall shear stress model, which should also apply to other pressure conditions and in particular to dry granular flows, where the interstitial gas plays a negligible role. The lateral boundaries are considered no-slip walls for the gas. The variations of temperature and density within the column are assumed negligible. The gas temperature is $T_g = 298$ K. The gas density is $\rho _g = 13.595$ kg m$^{-3}$. The gas dynamic viscosity is $\eta _g = 1.7982 \times 10^{-5}$ Pa s.

Figure 1. (a) Computational domain of the CFD-DEM simulations, and (b) bottom view of the CFD-DEM mesh.

The bed is filled with smooth, mono-disperse, ideally spherical particles of diameter $d_p = R_c/88 = 875 \times 10^{-6}$ m and density $\rho _p = 740$ kg m$^{-3}$. The minimum fluidisation velocity of the particles is 0.095 m s$^{-1}$ according to the correlation of Wen & Yu (Reference Wen and Yu1966), and 0.10 m s$^{-1}$ according to the correlation of Thonglimp, Hiquily & Laguerie (Reference Thonglimp, Hiquily and Laguerie1984). For the particles, the lateral boundaries and the gas inlet act as walls, and are treated identically. Thus no particle can leave the bed during the simulation. Two cases from Fede et al. (Reference Fede, Moula, Ingram, Dumas and Simonin2009, Reference Fede, Simonin and Ingram2016) are considered, as reported in table 1: case F3, with total solid mass $M_p = 2.5$ kg and fluidisation velocity $V_f = 0.32\ {\rm m}\ {\rm s}^{-1}$; and case F4, with total solid mass $M_p = 3.8$ kg and fluidisation velocity $V_f = 0.16\ {\rm m}\ {\rm s}^{-1}$. A frictionless variant of case F3 is also considered, in which the sliding friction between particles and walls is neglected. The two cases are in the bubbling fluidisation regime. The bed is denser and associated with smaller bubbles in case F4, given its larger total solid mass and lower fluidisation velocity. The average bed height is $h_b \approx 4.8R$ in case F3, and $h_b \approx 6.6R$ in case F4.

Table 1. Operating conditions and numerical parameters of the CFD-DEM simulations.

2.2. Governing equations and numerical method

Following the CFD-DEM approach, the particulate phase is modelled as a discrete set of particles, whereas the gaseous phase is modelled as a continuous medium that can penetrate the particles. The motion of the particles is modelled using Newton's second law of motion, with a soft-sphere treatment of the contact forces (Cundall & Strack Reference Cundall and Strack1979). The motion of the gas is modelled using a large-eddy simulation formalism, under the low Mach number approximation and without density or temperature variations. A detailed description of the governing equations used for the gaseous and solid phases is found in Dufresne et al. (Reference Dufresne, Moureau, Lartigue and Simonin2020) and Nigmetova et al. (Reference Nigmetova, Masi, Simonin, Dufresne and Moureau2022). The most relevant modelling hypotheses are recalled below for clarity. First, one may note that particle rotation is not considered, for consistency with the two-fluid modelling assumptions discussed in § 4.1. Although the effect of particle rotation is neglected in the present study, it would be valuable to carry out an investigation of its effect in the future, as it might affect the dynamics of the particles near the side walls. Second, the normal contact force exerted on a particle by another particle or by a wall is computed using damped linear spring model, while the tangential contact force is computed using a Coulomb sliding friction model. The collision model is characterised by spring stiffness $k_n$, the spring damping coefficient $\gamma _n$, and the coefficient of sliding friction $\mu$. In the soft-sphere approach, the spring stiffness is generally not determined based on the Young's modulus of the particle material, as this would lead to an impractically small collision time (van der Hoef et al. Reference van der Hoef, Ye, van Sint Annaland, Andrews, Sundaresan and Kuipers2006)

(2.1)\begin{equation} T_c = {\rm \pi}\sqrt{m_{eff}/k_n}. \end{equation}

In the present CFD-DEM simulations, the spring stiffness is set to $k_{n,p} = 200\ {\rm N}\ {\rm m}^{-1}$ for particle–particle contacts. Since particle–particle contacts are limiting for the simulation time step, a higher value $k_{n,w} = 300\ {\rm N}\ {\rm m}^{-1}$ is used for for particle–wall contacts, without additional computational cost. The spring damping coefficient $\gamma _n$ is given by

(2.2)\begin{equation} \gamma_n = \frac{- \ln e_c}{\sqrt{{\rm \pi}^2 + (\ln e_c)^2}}\,\sqrt{k_n/m_{eff}}, \end{equation}

where $e_c$ is the normal restitution coefficient associated with the contacts. For particle–wall contacts, the restitution coefficient is denoted $e_w$ and set to $e_w = 1$, which corresponds to purely elastic contacts. For particle–particle contacts, the restitution coefficient is denoted $e_c$ and set to a lower value $e_c = 0.9$ to take into account the energy dissipated during particle collisions. As reported in table 1, CFD-DEM simulations are performed with wall sliding friction coefficient $\mu _w = 0$ (frictionless wall contacts) and $\mu _w = 0.3$ (non-zero wall sliding friction). The coefficient of sliding friction is set to zero for particle–particle contacts in all cases.

The simulations are performed using a finite-volume method. The gas is advanced using the time-staggered projection method of Chorin (Reference Chorin1968) and the explicit TRK4 scheme of Kraushaar (Reference Kraushaar2011), which combines linearly a fourth-order Runge–Kutta scheme and a two-step Taylor–Galerkin scheme (Colin & Rudgyard Reference Colin and Rudgyard2000). A fourth-order centred scheme is used for spatial discretisation. The simulation is performed on an O-grid mesh that contains 350 000 cells in total, as represented in figure 1(b). The cell size in the vertical direction ($z$) is regular and equal to $\Delta z = 5 d_p$. The typical cell size in the radial direction ($r$) is also $\Delta r = 5 d_p$. The grid sensitivity analysis performed in Nigmetova et al. (Reference Nigmetova, Masi, Simonin, Dufresne and Moureau2022) has shown that this cell size is relevant. In the outer region, the number of cells in the angular direction ($\theta$) is $N_\theta = 80$. This implies that the circular geometry of the column is approximated by an octacontagon (A.4).

The particles are advanced using an explicit second-order Runge–Kutta scheme for time advancement. A linked-cell data structure is used to search for potential collision partners (Lubachevsky Reference Lubachevsky1991; Komiwes et al. Reference Komiwes, Mege, Meimon and Herrmann2006). A Voronoi-region decomposition of the wall features is used to solve the wall contacts. This is performed using the massively parallel code YALES2 (Moureau, Domingo & Vervisch Reference Moureau, Domingo and Vervisch2011). A detailed description of the CFD-DEM methodology is given in Dufresne et al. (Reference Dufresne, Moureau, Lartigue and Simonin2020). The gas is advanced in time using a time step $\Delta t_f$. The particles are advanced in time using a time step $\Delta t_p < \Delta t_f$, during which the gas advancement is frozen. The computation of the fluid and particle time steps is detailed in Nigmetova et al. (Reference Nigmetova, Masi, Simonin, Dufresne and Moureau2022).

3.1. Continuum modelling

To represent a dispersed phase as a continuous medium, one may consider that the motion of the particles is analogous to the motion of the particles in a gas. Namely, continuum equations and constitutive relationships relevant for an Eulerian modelling of the particles may be derived following the approach used in formal kinetic theory (Chapman & Cowling Reference Chapman and Cowling1990), that is, from the transport equation of the particle probability density in phase space (Buyevich & Shchelchkova Reference Buyevich and Shchelchkova1979; Morioka & Nakajima Reference Morioka and Nakajima1987; Reeks Reference Reeks1991; Simonin Reference Simonin2000; Zaichik, Oesterlé & Alipchenkov Reference Zaichik, Oesterlé and Alipchenkov2004), which is a Boltzmann-type equation that accounts for the interaction with the gaseous phase, external forces and particle–particle collisions. In this work, the one-particle phase-space density is denoted $f_p(\boldsymbol {v}, \boldsymbol {x} ; t)$ and defined such that $f_p(\boldsymbol {v}, \boldsymbol {x} ; t) \,{{\rm d}\boldsymbol {v}} \,{{\rm d} \kern0.7pt \boldsymbol {x}}$ is the expected number of particles whose centre is within a range ${{\rm d} \kern0.7pt \boldsymbol {x}}$ of location $\boldsymbol {x}$, and within a range ${{\rm d}\boldsymbol {v}}$ of velocity $\boldsymbol {v}$. Since the particles have a large inertia, we can neglect the effect of the fluid turbulence on the particle motion. For any particle-valued variable $\varphi$, the expected value of $\varphi$ at time $t$ may be expressed as

(3.1)\begin{equation} \langle\varphi\rangle(\boldsymbol{x}, t) = \frac{1}{n_p(\boldsymbol{x}, t)} \int_{-\infty}^{+\infty} \varphi(\boldsymbol{v}, \boldsymbol{x} ; t)\,f_{p}(\boldsymbol{v}, \boldsymbol{x} ; t) \,{{\rm d}\boldsymbol{v}}, \end{equation}

where $n_p(\boldsymbol {x}, t) = \int _{-\infty }^{+\infty } f_p(\boldsymbol {v}, \boldsymbol {x} ; t) \,{{\rm d}\boldsymbol {v}}$ is the number of particle centres per unit volume. The expectation $\langle \cdot \rangle$ is conditioned implicitly on a given realisation of the gaseous phase, and corresponds to an ensemble average over a distribution of particulate-phase realisations (Février, Simonin & Squires Reference Février, Simonin and Squires2005). Conservation equations associated with the particle velocity moments may be derived. Within a particulate medium, up to but excluding wall boundaries, the transport equations of the mean particle density $n_p$, particle velocity $\langle \boldsymbol {u}_p\rangle$ and particle velocity covariance $\boldsymbol{\mathsf{R}}_{\boldsymbol {p}} = \langle \boldsymbol {u}'' \otimes \boldsymbol {u}'' \rangle$ may be expressed as

(3.2)\begin{align} \partial_t (n_p m_p) + \boldsymbol{\nabla} \boldsymbol{\cdot} (n_p m_p \langle\boldsymbol{u}_{p}\rangle) & = 0, \end{align}

(3.3) \begin{align} \partial_t (n_p m_p \langle\boldsymbol{u}_{p}\rangle) + \boldsymbol{\nabla} \boldsymbol{\cdot} (n_p m_p \langle\boldsymbol{u}_{p}\rangle \otimes \langle\boldsymbol{u}_{p}\rangle) & =- \boldsymbol{\nabla}\boldsymbol{\cdot} (n_p m_p \boldsymbol{\mathsf{R}}_{\boldsymbol{p}}) + n_p \,\boldsymbol{f}_{\!p} \nonumber\\ & \quad + n_p \left\langle\,\boldsymbol{f}_{\!g \rightarrow p}\right\rangle + n_p m_p \boldsymbol{g},\end{align}

(3.4) \begin{align}\partial_t (n_p m_p \boldsymbol{\mathsf{R}}_{\boldsymbol{p}}) + \boldsymbol{\nabla}\boldsymbol{\cdot} ( n_p m_p \boldsymbol{\mathsf{R}}_{\boldsymbol{p}} \otimes \langle\boldsymbol{u}_{p}\rangle) & = {}- \boldsymbol{\nabla} \boldsymbol{\cdot} (n_p m_p \boldsymbol{\mathsf{T}}_{\boldsymbol{p}}) \nonumber\\ & \quad - n_p m_p \boldsymbol{\mathsf{R}}_{\boldsymbol{p}} \boldsymbol{\cdot} \boldsymbol{\nabla} \langle\boldsymbol{u}_{p}\rangle + n_p m_p\,\boldsymbol{\nabla} \langle\boldsymbol{u}_{p}\rangle \boldsymbol{\cdot} \boldsymbol{\mathsf{R}}_{\boldsymbol{p}}\nonumber\\ & \quad + n_p \left\langle\,\boldsymbol{f}_{\!p} \otimes \boldsymbol{u}_{p}''\right\rangle + n_p \left\langle\boldsymbol{u}_{p}'' \otimes \boldsymbol{f}_{\!p}\right\rangle\nonumber\\ & \quad + n_p \left\langle\,\boldsymbol{f}_{\!g \rightarrow p} \otimes \boldsymbol{u}_{p}''\right\rangle + n_p \left\langle\boldsymbol{u}_{p}''\otimes \boldsymbol{f}_{\!g \rightarrow p} \right\rangle, \end{align}

where $\boldsymbol {u}_{p}'' = \boldsymbol {u}_{p} - \langle \boldsymbol {u}_{p}\rangle$ is the particle velocity fluctuation, $\boldsymbol{\mathsf{T}}_{\boldsymbol {p}} = \langle \boldsymbol {u}_{p}'' \otimes \boldsymbol {u}_{p}'' \otimes \boldsymbol {u}_{p}''\rangle$ is the third-order particle velocity covariance, $n_p \boldsymbol {f}_{\!g \rightarrow p}$ is the total force per unit volume exerted by the gas, and $n_p \boldsymbol {f}_{\!p}$ is the total force per unit volume due to contacts with other particles.

For convenience, the particle–particle contact force $\boldsymbol {f}_{\!p}$ is expressed formally as the divergence of a particle–particle contact-stress tensor ${\boldsymbol {\varSigma }^{cnt}_{\boldsymbol {p}}}$. This lets us rewrite (3.3) as

(3.5) \begin{align} \partial_t (n_p m_p \langle\boldsymbol{u}_{p}\rangle) + \boldsymbol{\nabla} \boldsymbol{\cdot} (n_p m_p \langle\boldsymbol{u}_{p}\rangle \otimes \langle\boldsymbol{u}_{p}\rangle) & =- \boldsymbol{\nabla}\boldsymbol{\cdot} ((\boldsymbol{\varSigma}^{kin}_{\boldsymbol{p}})^{\rm T} + (\boldsymbol{\varSigma}^{cnt}_{\boldsymbol{p}})^{\rm T})\nonumber\\ & \quad + n_p \left\langle\,\boldsymbol{f}_{\!g \rightarrow p}\right\rangle + n_p m_p \boldsymbol{g}, \end{align}

where $\boldsymbol {\varSigma }^{kin}_{\boldsymbol {p}} = n_p m_p \boldsymbol{\mathsf{R}}_{\boldsymbol {p}}$ is by definition the kinetic stress of the particles. We thereby identify two classes of stresses within the particulate medium, as represented schematically in figure 2(a). The kinetic stress represents the transport of momentum by the translational motion of the particles across a given surface, due to the particle agitation. The contact stress represents the transfer of momentum due to contacts with other particles. In a dense fluidised bed, the stress associated with particle contacts can be expressed formally as the sum of contributions from two modes of contact: short-duration collisions and enduring contacts. Short-duration collisions are the dominant mode of contact in moderately dense beds. In a moderately dense bed, they can be treated as instantaneous collisions and are generally modelled following the KTGF. Enduring contacts are the dominant mode of contact in densely packed beds. In densely packed beds, they can be treated as semi-permanent frictional contacts and are generally modelled using empirical constitutive relationships. In the following, the stress that results from short-duration collisions is referred to as the collisional stress and denoted $\boldsymbol {\varSigma }^{col}_{\boldsymbol {p}}$, and the stress that results from enduring contacts is referred to as the frictional stress and denoted $\boldsymbol {\varSigma }^{fr}_{\boldsymbol {p}}$, such that

(3.6)\begin{equation} \boldsymbol{\varSigma}^{cnt}_{\boldsymbol{p}} = \boldsymbol{\varSigma}^{col}_{\boldsymbol{p}} + \boldsymbol{\varSigma}^{fr}_{\boldsymbol{p}}. \end{equation}

In practice, only the contact stress $\boldsymbol {\varSigma }^{cnt}_{\boldsymbol {p}}$ can be measured in the CFD-DEM simulation, and not the separate contributions $\boldsymbol {\varSigma }^{col}_{\boldsymbol {p}}$ and $\boldsymbol {\varSigma }^{fr}_{\boldsymbol {p}}$. Indeed, $\boldsymbol {\varSigma }^{cnt}_{\boldsymbol {p}}$ corresponds to the contact force per unit area between particles lying on the opposite side of a surface. The computation of the stresses in the CFD-DEM simulation is detailed in (A.2) and (A.3). The collisional and frictional contributions $\boldsymbol {\varSigma }^{col}_{\boldsymbol {p}}$ and $\boldsymbol {\varSigma }^{fr}_{\boldsymbol {p}}$ stem from an abstract decomposition in the TFM, and we are not currently able to use an objective criterion to identify the contacts that should be viewed as collisional or frictional in a CFD-DEM simulation.

Figure 2. Schematic representation of the different classes of (a) particle momentum fluxes across a surface $A_{bulk}$ in direction $n$ at a distance $x_{n,bulk} = d_p/2 + \Delta n$ from the wall surface, and (b) momentum flux exchanges between the particles and the wall. The contact flux is the sum of the collisional flux and the frictional flux.

An analogous classification of the stresses may be established on the wall boundaries, as represented schematically in figure 2(b). The particles in contact with a wall may undergo short-duration collisions and enduring contacts with the wall. The stress vector exerted by an element of wall surface of normal $\boldsymbol {n}$ on the particles, denoted $\boldsymbol {\sigma }^{tot}_{w}$, can thus be expressed formally as the sum of contributions from the two possible modes of contact:

(3.7)\begin{equation} \boldsymbol{\sigma}^{tot}_{w} = \boldsymbol{\sigma}^{col}_{w} + \boldsymbol{\sigma}^{fr}_{w}, \end{equation}

with $\boldsymbol {\sigma }^{col}_{w}$ the collisional particle–wall stress and $\boldsymbol {\sigma }^{fr}_{w}$ the frictional particle–wall stress. There is no kinetic contribution to the particle–wall stress. Indeed, a non-zero kinetic particle–wall stress would imply that some particles leave the domain, i.e. either pass through the wall or deposit onto the wall.

3.2. Continuity of the stress at the wall

Let us consider particulate stress in the vicinity of the location $\boldsymbol {x}_w$ on the wall surface. The unit vector $\boldsymbol {n}$ is normal to the wall at this location, and the unit vector $\boldsymbol {t}$ is parallel to the wall surface and aligned with the tangential velocity of the particles above the wall, $\boldsymbol {t} = \langle u_{p,t}\rangle /\| \langle u_{p,t}\rangle \|$, with $\langle u_{p,t}\rangle = \langle \boldsymbol {u}_{p}\rangle - \langle \boldsymbol {u}_{p,n}\rangle \boldsymbol {n}$ and $\langle \boldsymbol {u}_{p,n}\rangle = \langle \boldsymbol {u}_{p}\rangle \boldsymbol {\cdot } \boldsymbol {n}$. The transverse unit vector $\boldsymbol {s}$ is taken such that $\{\boldsymbol {t}, \boldsymbol {n}, \boldsymbol {s}\}$ is a local orthonormal basis. Assuming hard-sphere collisions, no particle centres are found below $\check {r}_w = 1/2$ at the location $\boldsymbol {x}_{wb} = \boldsymbol {x}_w + (d_p/2)\boldsymbol {n}$. By construction, the stress vector (i.e. the traction) associated with a wall-parallel plane within the particulate medium, $\boldsymbol {\sigma }^{tot}_{p,n} = (\boldsymbol {\varSigma }^{tot}_{\boldsymbol {p}}(\boldsymbol {x}, t))^{\rm T} \boldsymbol {\cdot } \boldsymbol {n}$, tends to the particle–wall stress when approaching the wall boundary. Hence we have

(3.8)\begin{equation} \boldsymbol{\sigma}^{tot}_{w}(\boldsymbol{x}_{wb}, t) = \lim_{\varepsilon \rightarrow 0^+} \boldsymbol{\sigma}^{tot}_{p,n}(\boldsymbol{x}_{wb} + \varepsilon \boldsymbol{n}, t). \end{equation}

For simplicity, we use the notation $\varphi (\boldsymbol {x}_{wb})$ to denote for any field $\varphi$ the limit $\displaystyle \lim _{\varepsilon \rightarrow 0^+} \varphi (\boldsymbol {x}_{wb} + \varepsilon \boldsymbol {n})$. This leads to no confusion since particle fields are not defined at the location $\boldsymbol {x}_{wb}$. We may thereby express as

(3.9)\begin{equation} \boldsymbol{\sigma}^{tot}_{w}(\boldsymbol{x}_{wb}, t) = \boldsymbol{\sigma}^{tot}_{p,n}(\boldsymbol{x}_{wb}, t) = \boldsymbol{\sigma}^{kin}_{p,n}(\boldsymbol{x}_{wb}, t) + \boldsymbol{\sigma}^{cnt}_{p,n}(\boldsymbol{x}_{wb}, t) \end{equation}

the correspondence between the classes of stress within the particulate medium and at the wall. Note that the vector $\boldsymbol {\sigma }^{cnt}_{p,n}$ corresponds to the stress associated with particle–particle contacts, which includes in general both a collisional contribution and a frictional contribution, but includes a frictional contribution only at the location $x_{wb}$ due to the closeness of the wall. Equation (3.9) can be simplified in the limit case of a moderately dense bed or a densely packed bed. Indeed, the particles in contact with a wall are unlikely to simultaneously be in contact with a particle in moderately dense beds, such that (3.9) simplifies to $\boldsymbol {\sigma }^{col}_{w}(\boldsymbol {x}_{wb}, t) = \boldsymbol {\sigma }^{kin}_{p,n}(\boldsymbol {x}_{wb}, t)$, which expresses the collisional stress at the wall based on the kinetic particle stress within the particulate medium. Conversely, the particle velocity covariance can be neglected in densely packed beds, such that (3.9) simplifies to $\boldsymbol {\sigma }^{fr}_{w}(\boldsymbol {x}_{wb}, t) = \boldsymbol {\sigma }^{fr}_{p,n}(\boldsymbol {x}_{wb}, t)$, which expresses the frictional stress at the wall based on the frictional stress within the medium.

It should be noted that the continuity of the stress at the wall applies to only the total stress, as the relative contributions of the kinetic, collisional or frictional stress vary dramatically near the contact point. The near-wall variations of the particle stresses are measured in the CFD-DEM simulations on a series of octacontagonal surfaces at a varying distance $d_p \check {r}_{w,bulk}$ of the lateral wall (A.4). Numerically, we find that the particle–wall normal stress $\sigma ^{tot}_{w,n}$ is close to the total normal particle stress close to the wall (figure 3), which corroborates the continuity of the stress at the wall. However, the kinetic particle stress increases rapidly as the distance to the wall tends to $d_p/2$, while the contact particle stress drops to a small value. At the limit, the kinetic contribution is predominant because most particle–wall contacts are collisional short-duration contacts. Indeed, the contact particle stress at the wall is by construction related to the number of particles infinitesimally close to the wall, which is zero if there are no frictional contacts, as this implies a finite wall particle density. The normal contact particle stress $\boldsymbol {\varSigma }^{cnt}_{p,nn}(\boldsymbol {x}_{wb}, t) = \boldsymbol {n} \boldsymbol {\cdot } \boldsymbol {\sigma }^{cnt}_{p,n}(\boldsymbol {x}_{wb}, t)$ at the height $\check {r}_w = 1/2$ is thus related to the presence of frictional contacts. In the CFD-DEM simulations F3 and F4, that contribution does not exceed 15 % of the wall-normal particle–wall stress at any height of the bed, indicating that enduring contacts are only a small proportion of the contacts. This is confirmed by the direct measurements of the particle–wall contact duration in the CFD-DEM simulation. The contribution of contacts of varying durations to the mean particle–wall stress according to the CFD-DEM simulations is provided in figure 4. Given that the duration of the contacts $T_c$ depends on the spring constant, which is non-physical in CFD-DEM, the figure should be regarded with some care and used for a qualitative assessment rather than a precise reading of the values. It may be seen that the stress increases with elevation from the bottom of the bed, because the particle concentration reaches a maximum in the upper part of the bed. In case F3, more than 90 % of the mean particle–wall stress is due to contacts that are less than three CFD-DEM collision times $T_c$ in duration. Enduring contacts are more important in case F4 due to the larger solid volume fraction and lower particle velocity of that case. Nevertheless, contacts less than seven CFD-DEM collision times $T_c$ account for more than 80 % of the mean particle–wall stress throughout the bed height (figure 4). This suggests that, from a modelling point of view, the collisional contribution to the particle–wall stress is dominant in the two CFD-DEM simulations.

Figure 3. Radial profile very close to the wall in the cases (a) F3 and (b) F4 of the total normal particle stress $\boldsymbol {\varSigma }^{cnt}_{p,nn}$, normal kinetic particle stress $\boldsymbol {\varSigma }^{kin}_{p,nn}$ and normal contact particle stress $\boldsymbol {\varSigma }^{cnt}_{p,nn}$. The filled square symbol at $\check {r}_{w} = 1/2$ indicates the corresponding particle–wall normal stress $\sigma ^{tot}_{w,n}$. (a) Upper part of the bed, $z/R = 3.45$. (b) Bottom of the bed, $z/R = 0.05$.

Figure 4. Vertical profiles in case F4 of the mean particle–wall normal stress $\sigma ^{tot}_{w,r}$ associated with particle–wall contacts that are shorter than a duration. The durations are given in units of the CFD-DEM collision time $T_c$ (see (2.1)), which corresponds to the duration of a short-duration collision between two particles in the CFD-DEM simulation.

The increase in the kinetic normal particle stress $\boldsymbol {\varSigma }^{kin}_{p,nn} = n_p m_p R_{p,nn}$ near the particle–wall contact point can also be explained by a large increase in the particle-centre density near a distance $d_p/2$ from the wall (see figure 6). Indeed, the wall surface has a large influence on the particle statistics, which is related to the quasi-two-dimensional arrangement of the particles in the vicinity of the wall. Section 3.3 describes this effect based on the numerical results of the CFD-DEM simulations, and discusses the prospect of taking into account these wall-specific effects in a TFM simulation and the implications for the modelling of the particle–wall shear stress.

3.3. Near-wall arrangement of the particles

In the vicinity of the container walls, the particle distribution and dynamics are constrained by the wall surface, and therefore display behaviour different to that in the internal region. This subsection analyses these effects for a dense fluidised bed in the bubbling regime, using the CFD-DEM simulations presented in § 2.

Figures 5(a,c) provide a visualisation of the particle arrangement near the wall in the simulation F3. The particles close to cylindrical walls tend to form particle assemblies with a two-dimensional hexagonal arrangement. These assemblies may be observed for each time step on a significant portion of the wall surface, but are local and transient. Indeed, following the passage of a bubble (figure 5c), the flow is locally more dilute and the arrangement of the particles is broken for a brief duration. The influence of the hexagonal arrangement of the particles lingers for a few layers away from the wall. However, as we move farther away from the wall, the influence of the wall diminishes, and the two-dimensional arrangement of the particles is much less orderly, because the particles are no longer constrained by the wall surface. Thus the arrangement of the particles near the wall differs from the arrangement of the particles in the internal region.

Figure 5. Close-up views in case F3 of (a,c) the first layer of particles near the wall ($0< z_r<1$), namely the particles whose centre is located within $z_r = (R_c-r_p)/d_p = 1$ particle diameter of the wall, and (b,d) a thin layer of particles away from the wall ($10< z_r<11$), namely the particles whose centre is $z_r = 10$–$11$ particle diameters away from the wall. (a,b) Instant A: Typical arrangement of the particles ($\alpha _p = 0.59$). (c,d) Instant B: Arrangement of the particles following the passage of a bubble ($\alpha _p = 0.26$).

The structural properties of the near-wall region are consistent with the structural properties of the near-wall region in a fixed randomly packed bed of spheres, which have been described both experimentally and numerically by various authors (Rocke Reference Rocke1971; Tingate Reference Tingate1973; Mueller Reference Mueller1997; Abreu et al. Reference Abreu, Macias-Salinas, Tavares and Castier1999; Zhang et al. Reference Zhang, Thompson, Reed and Beenken2006; Wensrich Reference Wensrich2012; Burtseva et al. Reference Burtseva, Valdez, Werner and Petranovskii2015; von Seckendorff & Hinrichsen Reference von Seckendorff and Hinrichsen2021). The effect of the wall can be characterised by the radial distribution function, as detailed in the supplementary material available at https://doi.org/10.1017/jfm.2024.36. The influence of the wall on the particle distribution can also be characterised by the near-wall radial evolution of the solid volume fraction. Figure 6(a) examines, in case F3, the radial variation of the particle-centre density near the wall for a series of nested octacontagonal shells of thickness $d_p/100$, and figure 6(b) examines the resulting solid volume fraction within each octacontagonal shell (see (A.4)). Due to the wall, the density of particle centres is not uniform. The particles are more likely to be located near specific layers, which after the first period are equally spaced, and less likely to be located in between. Accordingly, there are large oscillations in the solid volume fraction near the wall. Both the oscillations of particle-centre density and solid volume fraction become negligible at a distance $\check {r}_w = 6$ particle diameters away from the wall. Note that since the particle distribution is two-dimensional near the wall, the solid volume fraction on the first oscillation, which exceeds 0.74 in the CFD-DEM, should be compared to that of a circle packing. The period of the oscillations is approximately $0.87 d_p$. This value is consistent with the oscillation period measured experimentally in a fixed randomly packed bed (von Seckendorff & Hinrichsen Reference von Seckendorff and Hinrichsen2021). Figure 6(b) provides a comparison of the CFD-DEM solid volume fraction with the fixed-bed correlation of Mueller (Reference Mueller1992), which can be expressed in the present case as

(3.10)\begin{equation} \alpha_p (\check{r}_w) = \alpha_{inf}\left(1 - {\rm J}_0((7.45 - 11.25/\lambda) \check{r}_w) \exp(-(0.315 - 0.725/\lambda) \check{r}_w)\right)\!, \end{equation}

where ${\rm J}_0$ is the zeroth-order Bessel function of the first kind, $\check {r}_w = (R_c-r)/d_p$ is the normalised distance to the wall, and $\lambda = 2R/d_p = 176$ is the cylinder to particle diameter ratio. Note that we do not use the correlation of Mueller (Reference Mueller1992) to set the bulk solid volume fraction $\alpha _{inf}$, but in place use the CFD-DEM value. The correlation predicts the CFD-DEM oscillations quite well very near to the wall, but underestimates the decrease in the oscillation amplitude away from the wall, compared to the CFD-DEM simulation.

Figure 6. Radial profile at the height $z/R_c=3.45$ in terms of the normalised distance to the wall $\check {r}_w = (R_c - r)/d_p$, in case F3, of (a) the number of particles per unit volume $n_p$ normalised by the mean bed solid fraction, for a series of nested octacontagonal shells of thickness $d_p/100$, and (b) the resulting solid volume fraction $\alpha _p$ within each octacontagonal shell. The solid volume fraction predicted by the correlation of Mueller (Reference Mueller1992) is given for reference.

The dynamics of the particles is also influenced by the wall. For instance, the mean wall-normal particle velocity is oscillating following the variations of particle density near the wall. The mean radial motion of the particles tends to be small within the layers where a large density of particles is found, and larger in the regions in between, where few particles are found (figure 7a). The mean vertical velocity of the particles exhibits similar oscillations. Namely, the fall of the particles in contact with the wall and in the layers where a large density of particles is found tends to be slow compared to the depleted regions. This can be explained by the fact that an orderly three-dimensional arrangement of the particles hinders motion and forbids the presence of particles in between layers. Similarly, the wall-normal particle velocity variance $R_{nn}$ (figure 7b), the vertical particle velocity variance $R_{zz}$ (figure 7d) and the vertical-normal particle velocity covariance $R_{nz}$ (figure 7e) exhibit oscillations and are greater in the regions where few particles are found, since they are associated with a locally less orderly and more dilute flow. The large oscillations of the particle-centre density and of the particle velocity variance are of interested for the modelling of the particle–wall shear stress since they are involved directly in the expression for the kinetic particle stress.

Figure 7. Radial profile at the height $z/R_c=3.45$ in terms of the normalised distance to the wall $\check {r}_w = (R_c - r)/d_p$, in case F3, for a series of nested octacontagonal shells of thickness $d_p/500$ near the wall surface, of (a) the mean wall-normal velocity $\langle u_{p,n} \rangle$, (b) the total wall-normal velocity variance $R_{nn}$, (c) the mean vertical velocity $\langle u_{p,z} \rangle = \langle u_{p,t} \rangle$, (d) the total vertical velocity variance $R_{zz}$, and (e) the total vertical-normal velocity covariance $R_{nz}$ within each octacontagonal shell.

In view of these results, it is clear that the vicinity of the container walls is associated with a specific particle dynamics that is not found in the internal region. The prospect of reproducing this complex near-wall dynamics in a TFM simulation, where the particulate phase is represented as a continuous medium, represents a considerable challenge. First, this would be costly in terms of computational resources, as resolving the near-wall spatial variations in a TFM simulation would require a fine mesh, with a cell size much smaller than a particle diameter. Second, this would be demanding in terms of modelling, as specific source terms would need to be included in the momentum and kinetic energy equations to account for the modification of the particle phase-space density in the near-wall region. In that context, we consider that the Eulerian field of a TFM simulation should filter out these small-scale variations. The continuum equations cannot be resolved too close to the wall in a two-fluid simulation, because certain effects are not accounted for in current models. For instance, the ‘wall shelter effect’ identified by Sakiz & Simonin (Reference Sakiz and Simonin1998), and discussed further in Fede & Simonin (Reference Fede and Simonin2018), leads to the appearance of a mean effective force towards the wall, and causes a significant increase in the number of particle centres per unit volume $n_p$ for wall distances below $(3/2) d_p$. Therefore, in practice, the boundary-condition model for the particulate phase should seek to express the particle–wall shear stress in terms of flow variables that are relatively far from the wall, avoiding values near the particle–wall contact point. Let us consider a location $\boldsymbol {x}_{bulk}$ that is sufficiently far from the wall. A model for the particle–wall shear stress requires, in particular, the modelling of the kinetic particle stress at the location $\boldsymbol {x}_{wb}$. This entails providing a model for the particle-centre density $n_p(\boldsymbol {x}_{wb}, t)$ and the particle velocity covariance $\boldsymbol{\mathsf{R}}_{\boldsymbol {p}}$. For instance, this could be achieved by providing a model for the wall-amplification coefficient $\chi$, which we assume is a function of the solid volume fraction $\alpha _p(\boldsymbol {x}_{bulk})$, to relate $R_{p,nn}(\boldsymbol {x}_{wb}, t)$ to $R_{p,nn}(\boldsymbol {x}_{bulk}, t)$, i.e.

(3.11)\begin{equation} \chi(\alpha_p(\boldsymbol{x}_{bulk})) = \frac{n_p(\boldsymbol{x}_{wb}, t)\,m_p \langle v_p' v_p'\rangle(\boldsymbol{x}_{wb}, t)}{n_p(\boldsymbol{x}_{bulk}, t)\,m_p \langle v_p' v_p'\rangle(\boldsymbol{x}_{bulk}, t)}. \end{equation}

In addition, the modelling of the frictional term should be addressed. This would require dedicated models as no specific frictional model is available to account for the modification of the arrangement of the particles in contact with the wall due to the two-dimensional structure imposed by the wall. The present paper circumvents these modelling requirements by, instead, seeking a relationship between the particle–wall stress and the particle stress within the particulate medium.

3.4. Development of a TFM boundary condition

This subsection uses the momentum balance (3.5) and the results of the CFD-DEM simulation data to develop a TFM model for the particle–wall shear stress. The particle–wall shear stress is the component of the stress parallel to the wall surface. For simplicity, we assume that the particle–wall shear stress vector is aligned with the tangential velocity of the particles at a distance $d_p/2$ from the wall, that is, that the particle–wall stress can be decomposed formally as

(3.12)\begin{equation} \boldsymbol{\sigma}^{tot}_{w}(\boldsymbol{x}_{wb}, t) = \sigma^{tot}_{w,n}(\boldsymbol{x}_{wb}, t)\,\boldsymbol{n} + \sigma^{tot}_{w,t}(\boldsymbol{x}_{wb}, t)\,\boldsymbol{t}, \end{equation}

with $\sigma ^{tot}_{w,t}$ the particle–wall shear stress. This assumption is verified exactly in the purely collisional regime. Indeed, the wall-parallel kinetic particle stress close to the wall, $n_p(\boldsymbol {x}_{wb}, t)\,m_p (\boldsymbol{\mathsf{R}}_{\boldsymbol {p}}(\boldsymbol {x}_{wb}, t) \boldsymbol {\cdot } \boldsymbol {n} - R_{p,nn}(\boldsymbol {x}_{wb}, t)\,\boldsymbol {n})$, is aligned theoretically with the tangential velocity $\langle \boldsymbol {u}_{p,t}\rangle$. If there are enduring contacts, then the alignment of the tangential force exerted by a single particle onto the wall and its tangential velocity is verified for each particle time step in the CFD-DEM simulation. However, assumption (3.12) is verified only if the correlation between the normal force exerted by the wall and the orientation of the instantaneous tangential particle velocity is not strong. Figure 8 confirms that the alignment between the particle–wall shear stress and the cell-averaged tangential particle velocity is generally high in the CFD-DEM simulation.

Figure 8. Relationship between the orientation of the particle–wall shear stress and the cell-averaged tangential particle velocity for various instants and locations along the lateral wall surface, where the orientation is defined as the angle between a vector and the unit vector $\boldsymbol {e}_\theta$. The red line is the identity. Due to the symmetry of the angle space, the clusters of points at the bottom left, bottom right, top left and top right of the figure are all close to each other and to the origin.

A useful relationship for modelling the particle–wall shear stress can be derived from the balance of momentum (3.5) in the vicinity of a wall, in the frame of a quasi-parallel equilibrium boundary layer assumption. By integrating (3.5) over a control volume extending in the wall-normal direction from the near-wall location $\varphi (\boldsymbol {x}_{wb} + \varepsilon \boldsymbol {n})$ to the location $\boldsymbol {x}_{bulk} = \boldsymbol {x}_{wb} + (\Delta n)\boldsymbol {n}$, and infinitesimally small in the two tangential directions $\boldsymbol {t}$ and $\boldsymbol {s}$, we obtain in the limit of an infinitesimally small $\varepsilon$:

(3.13) \begin{align} \boldsymbol{\sigma}^{tot}_{w}(\boldsymbol{x}_{wb}, t) &= \boldsymbol{\sigma}^{kin}_{p,n}(\boldsymbol{x}_{bulk}, t) + \boldsymbol{\sigma}^{cnt}_{p,n}(\boldsymbol{x}_{bulk}, t) + n_p(\boldsymbol{x}_{bulk}, t)\,m_p (\langle\boldsymbol{u}_{p}\rangle \otimes \langle\boldsymbol{u}_{p}\rangle)(\boldsymbol{x}_{bulk}, t) \boldsymbol{\cdot} \boldsymbol{n}\nonumber\\ &\quad + \int_{d_p/2}^{d_p/2+\Delta n} \left( \partial_{t} (n_p m_p \langle\boldsymbol{u}_{p}\rangle) - n_p \left\langle\,\boldsymbol{f}_{\!g \rightarrow p}\right\rangle - n_p m_p \boldsymbol{g} \right) {{\rm d}\kern0.7ptx_n}\nonumber\\ &\quad + \int_{d_p/2}^{d_p/2+\Delta n} \frac{\partial{}}{\partial{x_t}} \left(n_p m_p \left(\langle\boldsymbol{u}_{p}\rangle \otimes \langle\boldsymbol{u}_{p}\rangle\right) + \boldsymbol{\varSigma}^{kin}_{\boldsymbol{p}} + \boldsymbol{\varSigma}^{cnt}_{\boldsymbol{p}}\right) \boldsymbol{\cdot} \boldsymbol{t} \,{{\rm d}\kern0.7pt x_n}\nonumber\\ &\quad + \int_{d_p/2}^{d_p/2+\Delta n} \frac{\partial{}}{\partial{x_s}} \left(n_p m_p \left(\langle\boldsymbol{u}_{p}\rangle \otimes \langle\boldsymbol{u}_{p}\rangle\right) + \boldsymbol{\varSigma}^{kin}_{\boldsymbol{p}} + \left(\boldsymbol{\varSigma}^{cnt}_{\boldsymbol{p}}\right)^{\rm T}\right) \boldsymbol{\cdot} \boldsymbol{s} \,{{\rm d}\kern0.7pt x_n}, \end{align}

where $x_n = (\boldsymbol {x} - \boldsymbol {x}_w) \boldsymbol {\cdot } \boldsymbol {n}$ denotes the wall-normal coordinate. The momentum balance can be simplified in the wall-normal direction. Indeed, assuming that the point $\boldsymbol {x}_{bulk} = \boldsymbol {x}_{wb} + (\Delta n)\boldsymbol {n}$ is not too far from the wall, the mean wall-normal particle velocity, the mean wall-normal gas velocity and the mean gas wall-normal pressure gradient can be assumed to be close to zero. This suggests that the contribution of the interaction with the gas ($\boldsymbol {f}_{\!g \rightarrow p}$) in (3.13) is negligible in the wall-normal direction. Assuming in addition a vertical wall, such that the contribution of gravity is zero in the wall-normal direction, allows us to neglect all the integral terms in (3.13). Hence the wall-normal particle stress balance reads

(3.14)\begin{equation} \sigma^{tot}_{w,n}(\boldsymbol{x}_{wb}, t) = \boldsymbol{\varSigma}^{tot}_{p,nn}(\boldsymbol{x}_{wb}, t) \approx \boldsymbol{\varSigma}^{kin}_{p,nn}(\boldsymbol{x}_{bulk}, t) + \boldsymbol{\varSigma}^{cnt}_{p,nn}(\boldsymbol{x}_{bulk}, t). \end{equation}

The wall-normal profiles of the particle stress balance, provided in figure 9, demonstrate that the total normal particle stress $\boldsymbol {\varSigma }^{tot}_{p,nn} = \boldsymbol {\varSigma }^{kin}_{p,nn} + \boldsymbol {\varSigma }^{cnt}_{p,nn}$ remains approximately constant and equal to the wall-normal particle–wall stress close to the wall. This validates the assumptions leading to (3.14) for the lateral wall of the CFD-DEM simulations. For every height in the CFD-DEM simulation, the contribution of the contact forces to the total particle stress dominates the kinetic contribution, except in the vicinity of $\check {r}_w = 1/2$. Near-wall oscillations of the normal kinetic particle stress and the normal contact particle stress may be observed (figure 9), following the oscillations of the particle-centre density and of the particle velocity variance described in § 3.3. However, these oscillations are relatively small in amplitude, and cancel out in the profile of the total normal particle stress.

Figure 9. Radial profile in case F3 of the total normal particle stress $\boldsymbol {\varSigma }^{cnt}_{p,nn}$, normal kinetic particle stress $\boldsymbol {\varSigma }^{kin}_{p,nn}$ and normal contact particle stress $\boldsymbol {\varSigma }^{cnt}_{p,nn}$. The filled square symbol at $\check {r}_{w} = 1/2$ indicates the corresponding particle–wall normal stress $\sigma ^{tot}_{w,n}$. (a) Upper part of the bed, $z/R = 3.45$. (b) Bottom of the bed, $z/R = 0.05$.

In the tangential direction, these approximations cannot be made, since the mean tangential velocity can be large at the wall. Thus the drag force exerted by the gas and the effect of the gas pressure gradient are not necessarily negligible. Furthermore, the weight of the solid in the control volume should be taken into account in the case of a vertical wall. However, the particle–wall shear stress can be related to the particle–wall normal stress using the coefficient of sliding friction of the particle–wall contacts:

(3.15)\begin{equation} \sigma^{tot}_{w,t}(\boldsymbol{x}_{wb}, t) = \mu_w\, \sigma^{tot}_{w,n}(\boldsymbol{x}_{wb}, t) = \mu_w\, \boldsymbol{\varSigma}^{tot}_{p,nn}(\boldsymbol{x}_{bulk}, t). \end{equation}

This expression will form the basis of the two-fluid modelling. For Coulombian contacts, the relationship (3.15) holds exactly for a single particle for each particle time step of the CFD-DEM simulation. In the general case, the relationship is valid if the orientations of the different particle–wall contacts are perfectly correlated with one another. This is examined in figure 10 based on the CFD-DEM data. A linear best fit of the CFD-DEM data gives $\sigma ^{tot}_{w,t}(\boldsymbol {x}_{wb}, t) = 0.268\, \sigma ^{tot}_{w,n}(\boldsymbol {x}_{wb}, t)$, with a slope 11 % smaller than the friction coefficient $\mu _w$. However, the distribution of the slope is skewed negatively, which suggests that a proportionality coefficient close to $\mu _w$ is preferable, as proposed in (3.15).

Figure 10. Relationship between the particle–wall shear stress and particle–wall normal stress for various instants and locations along the lateral wall surface. The red line is the 0.3 slope corresponding to the coefficient of sliding friction of the CFD-DEM simulation. The green line is the linear best fit, of slope 0.268.

In order to develop a model for the wall-normal particle–particle contact stress in (3.15), it is useful to first assume that the diagonal part of the particle–particle contact-stress tensor is distributed isotropically, that is, to make the approximation

(3.16)\begin{equation} \boldsymbol{\varSigma}^{tot}_{p,nn}(\boldsymbol{x}_{bulk}, t) = P_{p}^{tot} (\boldsymbol{x}_{bulk}, t), \end{equation}

with $P_{p}^{tot} = (1/3) \operatorname {tr}(\boldsymbol {\varSigma }^{tot}_{\boldsymbol {p}})$ the total granular pressure. This is a non-trivial approximation and could not be argued to be valid generally in all gas–solid fluidised beds. For the configurations F3 and F4, the diagonal particle–particle contact stress is, according to the CFD-DEM data, frequently close to isotropy, which supports the assumption of isotropy (figure 11). Furthermore, the two-fluid modelling paradigm used in the present work, which is based on a transport equation for the kinetic energy of particle agitation (§ 4.1), does not consider the anisotropy of the particle agitation and of the granular pressure. Using (3.16) and (3.15), a TFM for the particle–wall shear stress can be obtained by formulating a suitable model for the granular pressure. We use the same model as within the flow, since (3.16) expresses the particle–wall stress in terms of the granular pressure relatively far from the wall. Specifically, the granular pressure is decomposed formally into a kinetic pressure $P_{p}^{kin}$, a collisional pressure $P_{p}^{col}$, and a frictional pressure $P_{p}^{fr}$:

(3.17)\begin{equation} \boldsymbol{\sigma}^{tot}_{w}(\boldsymbol{x}_{wb}, t) = \left( P_{p}^{kin} (\boldsymbol{x}_{bulk}, t) + P_{p}^{col} (\boldsymbol{x}_{bulk}, t) + P_{p}^{fr} (\boldsymbol{x}_{bulk}, t) \right) \left( \boldsymbol{n} - \mu_w\, \frac{\boldsymbol{u}_{p,t}}{\| \boldsymbol{u}_{p,t} \|} \right)\!. \end{equation}

Using (3.17), the particle–wall stress is expressed in terms of the total granular pressure at some distance from the wall. The kinetic, collisional and frictional contribution in (3.17) are given functions of the solid volume fraction $\alpha _p$ and the kinetic energy of particle agitation $q_p^2$, and various models have been proposed in the literature to address these dependencies (Gu et al. Reference Gu, Ozel, Kolehmainen and Sundaresan2019; Si, Shi & Yu Reference Si, Shi and Yu2019). Note that the different contributions are related not to the mode of particle–wall contacts, but rather to the contacts within the flow. Thus each contribution should be included in general, regardless of the specific particle–wall interactions. Finally, it may be noted that although it has been developed in the context of a cylindrical column, the derivation of the present section has neglected the effect of curvature given that the diameter of the column is large compared to the diameter of a particle. The proposed model thus should also be applicable in the case of a flat wall.

Figure 11. Probability density associated with various distributions of the diagonal components of the particle–particle contact-stress tensor, characterised by the ratio of vertical–vertical ($zz$) component to pressure $\varSigma ^{cnt}_{p,zz} (\boldsymbol {x}_{bulk}, t)/P_{p} (\boldsymbol {x}_{bulk}, t)$ and the ratio of radial–radial ($rr$) component to pressure $\varSigma ^{cnt}_{p,rr} (\boldsymbol {x}_{bulk}, t)/P_{p} (\boldsymbol {x}_{bulk}, t)$, in the second cell off the wall in case F3. The dashed lines indicate that two of the diagonal components of the particle–particle contact-stress tensor are equal. Their intersection indicates isotropy.

The proposed model is related closely to several boundary-condition models from the literature. The model of Johnson & Jackson (Reference Johnson and Jackson1987), for a smooth flat wall, can be interpreted as a model that includes only the kinetic and frictional contributions to the particle–wall stress:

(3.18)\begin{equation} \boldsymbol{\sigma}^{tot}_{w}(\boldsymbol{x}_{wb}, t) = \left( P_{p}^{kin} (\boldsymbol{x}_{bulk}, t) + P_{p}^{fr} (\boldsymbol{x}_{bulk}, t) \right) \left( \boldsymbol{n} - \mu_w\,\frac{\boldsymbol{u}_{p,t}}{\| \boldsymbol{u}_{p,t} \|} \right)\!. \end{equation}

The collisional contribution has, however, been shown to be important in the present CFD-DEM simulations. The boundary-condition model of Jenkins (Reference Jenkins1992), in the all-sliding case, may be recast as

(3.19)\begin{equation} \boldsymbol{\sigma}^{tot}_{w}(\boldsymbol{x}_{wb}, t) = \left( \chi(\alpha_p(\boldsymbol{x}_{bulk}))\, P_{p}^{kin} (\boldsymbol{x}_{bulk}, t) \right) \left( \boldsymbol{n} - \mu_w\,\frac{\boldsymbol{u}_{p,t}}{\| \boldsymbol{u}_{p,t} \|} \right), \end{equation}

with $\chi (\alpha _p) = (1/2)\,g_0(\alpha _p)\,(1 + e_w)$, where $g_0$ denotes the radial distribution function, given by $g_0 = (1 - \alpha _p/\alpha _{p}^{max})^{-2.5 \alpha _{p}^{max}}$, with $\alpha _{p}^{max} = 0.64$ the close-packing limit (Lun & Savage Reference Lun and Savage1986). The same remark applies to the boundary-condition model of Schneiderbauer et al. (Reference Schneiderbauer, Schellander, Löderer and Pirker2012), which is based on a similar approach. By neglecting the frictional contribution in (3.17), the present boundary-condition model suggests $\chi (\alpha _p) = 1 + 2 \alpha _p\,g_0(\alpha _p)\,(1 + e_c)$. In the limit case of a dilute bed, the wall-amplification coefficient $\chi$ tends to $(1/2)(1+e_w)$ in the model of Jenkins (Reference Jenkins1992), whereas $\chi$ tends to 1 in the present model, which is consistent with a kinetic boundary-condition model (see (4.5) below). The present boundary-condition model tends to a kinetic boundary-condition model in the limit case of a dilute bed, and to a purely frictional model $\sigma ^{tot}_{w,t} = \mu _w P_p^{fr}$ in the limit of a densely packed bed. This thus provides a candidate general boundary condition that takes into account all modes of particle–wall contacts, whether collisional or frictional.

4.1. Two-fluid modelling

The Eulerian two-fluid method used in the present paper relies on a hybrid modelling approach (Morioka & Nakajima Reference Morioka and Nakajima1987), in which the continuous-phase equations are derived by phase ensemble averaging and the dispersed-phase equations are derived by using the KTGF supplemented with models for fluid and turbulent effects for the dispersed phase. The system of equations is based on the approach of Simonin (Reference Simonin2000) and includes a transport equation for the kinetic energy of particle agitation $q_p^2$. A detailed description of the governing equations is found in Ansart et al. (Reference Ansart, García-Triñanes, Boissière, Benoit, Seville and Simonin2017). The same modelling assumptions as in the CFD-DEM simulations are used to compute the mean gas–particle relaxation time scale (Gobin et al. Reference Gobin, Neau, Simonin, Llinas, Reiling and Sélo2003). The particle–particle normal restitution coefficient is set to $e_c = 0.9$, to be consistent with the CFD-DEM simulations.

The TFM wall boundary condition (3.17) developed in § 3, and several boundary conditions from the literature, have been assessed. Generally, the wall boundary condition for the particles should provide a closure for the particle momentum and agitation kinetic energy fluxes at the wall, or more precisely, for the particles whose centres are at distance $d_p/2$ from the wall. To ensure that there is no solid mass flow across the boundary, the wall-normal particle velocity is set to zero at the walls. Thus only the tangential particle momentum flux requires modelling. In other words, providing a wall boundary condition for the particles entails providing a model for (1) the particle–wall shear stress $\sigma ^{tot}_{w,t} = \eta _p\, \partial _n u_{p,t}|_{wall}$ and (2) the agitation kinetic energy flux $Q_w = \alpha _p \rho _p K_p^{{kin\text -col}}\,\partial _n q_p^2|_{wall}$. The fluxes $\sigma ^{tot}_{w,t}$ and $Q_w$ depend on how the discrete particles interact with the wall and in particular may be functions of the elastic normal and tangential restitution coefficients, the wall friction coefficient (Jenkins Reference Jenkins1992; Schneiderbauer et al. Reference Schneiderbauer, Schellander, Löderer and Pirker2012) and the wall roughness (Hui et al. Reference Hui, Haff, Ungar and Jackson1984; Konan et al. Reference Konan, Simonin and Squires2006). Presently, ideally smooth particles and walls are considered, since there is no roughness in the reference CFD-DEM simulations that we aim to reproduce in a TFM setting.

All the boundary-condition models investigated assume that the particle agitation kinetic energy flux is negligible, $Q_w = 0$. We will therefore focus on the tangential momentum flux. The following models are considered.

(4.1)\begin{gather} \text{No-slip:} \quad U_{p,t}|_{wall} = 0. \end{gather}

(4.2)\begin{gather}\text{Kinetic:} \quad \sigma^{tot}_{w,t} = \mu_w P_{p}^{kin}. \end{gather}

(4.3)\begin{gather}\text{Kinetic-frictional:} \quad \sigma^{tot}_{w,t} = \mu_w \left(P_{p}^{kin} + P_{p}^{fr}\right). \end{gather}

(4.4)\begin{gather}\text{Total-pressure (present):} \quad \sigma^{tot}_{w,t} = \mu_w \left(P_{p}^{kin} + P_{p}^{col} + P_{p}^{fr}\right)\!. \end{gather}

To close the boundary-condition models, the following expressions may be used for the kinetic, collisional and frictional pressures (Ansart et al. Reference Ansart, García-Triñanes, Boissière, Benoit, Seville and Simonin2017; Bennani et al. Reference Bennani, Neau, Baudry, Laviéville, Fede and Simonin2017):

(4.5)\begin{gather} P_p^{kin} = \alpha_p \rho_p \left(\tfrac{2}{3} q_p^2\right)\!, \end{gather}

(4.6)\begin{gather}P_p^{col} = 2 \alpha_p^2 \rho_p g_0\left(1 + e_c\right) \left(\tfrac{2}{3} q_p^2\right)\!, \end{gather}

(4.7)\begin{gather}P_p^{fr} = Fr\,\frac{\left(\alpha_p - \alpha_{p}^{min}\right)^{n_{\!J}}}{\left(\alpha_{p}^{max} - \alpha_p\right)^{m_{\!J}}} \left[\alpha_p \geq \alpha_{p}^{min}\right]\!, \end{gather}

with $g_0 = (1 - \alpha _p/\alpha _{p}^{max})^{-2.5 \alpha _{p}^{max}}$ the radial distribution function. Following the model of Johnson et al. (Reference Johnson, Nott and Jackson1990), the frictional pressure law is a rational function that vanishes for a loose-packed bed and increases asymptotically in the limit of a close-packed bed. The proportionality constant $Fr=0.05$, the exponents $n_{\!J}=2$ and $m_{\!J}=5$, and the activation threshold $\alpha _{p}^{min}$ are empirical parameters. The values proposed originally in Johnson et al. (Reference Johnson, Nott and Jackson1990) are used for $n_{\!J}$ and $m_{\!J}$. The activation threshold $\alpha _{p}^{min}$ is set to 0.63. This high value is relevant for the particles being used, which are ideally spherical and frictionless.

Following Fede et al. (Reference Fede, Simonin and Ingram2016), the no-slip boundary condition assumes the maximum flux transferred by the particles towards the wall. The kinetic boundary condition is a suitable boundary condition in dilute beds, as investigated in Sakiz & Simonin (Reference Sakiz and Simonin1999). It assumes instantaneous fully sliding particle–wall contacts on a smooth flat surface, with a Coulomb sliding friction coefficient $\mu _w$. The kinetic-frictional boundary condition assumes that the particle–wall stress may be expressed as the sum of a kinetic-collisional contribution and a frictional contribution, modelled based on the frictional pressure within the first cell off the wall. It corresponds to the boundary condition of Johnson & Jackson (Reference Johnson and Jackson1987) in the particular case where the specularity coefficient is neglected. This assumption is well justified considering that ideally the particles are spherical and the walls are smooth in the CFD-DEM simulation. The present total-pressure boundary condition can be derived assuming, following the CFD-DEM analysis of § 3, that the particle–wall stress may be expressed based on the total granular pressure, which is due to kinetic, collisional and frictional contributions. It should be noted that the choice of the distance from the wall at which the granular pressure is taken is not critical in the model formulation, given that the total granular pressure is constant in the vicinity of the wall and that the TFMs do not represent the near-wall oscillations of the particle statistics observed in the CFD-DEM simulations (§ 3.3). The total granular pressure within the first computational cell off the wall is therefore used for the sake of simplicity. Due to the high value of the activation threshold $\alpha _{p}^{min}$, the frictional term in (4.4) is not large and can be neglected. The frictional contribution may, however, be dominant in other flow configurations, in particular if the particles are rough or non-spherical.