1. Introduction
Fluidization is one of the most important solids processing technologies in industry and it is widely utilized in granular drying, coating, mixing, chemical reaction and combustion, due to its high heat and mass transfer efficiency (Rhodes Reference Rhodes2008). Extensive studies on fluidization have been conducted experimentally and computationally. Particles are assumed to be spherical in most of these studies, particularly computational studies (Pan et al. Reference Pan, Joseph, Bai, Glowinski and Sarin2002; Liu & Hrenya Reference Liu and Hrenya2018; Sippola et al. Reference Sippola, Kolehmainen, Ozel, Liu, Saarenrinne and Sundaresan2018); however, particle shape is known to have a significant impact on particle packing (Wouterse, Williams & Philipse Reference Wouterse, Williams and Philipse2007), particle flow (Guo et al. Reference Guo, Wassgren, Ketterhagen, Hancock, James and Curtis2012) and fluid–particle interaction (Hölzer & Sommerfeld Reference Hölzer and Sommerfeld2008). Thus, fluidization behaviour is also affected by particle shape.
Gas fluidization of non-spherical particles has attracted increasing research interest in recent years, due to the ubiquity of non-spherical particles in engineering practice. Significantly different fluidization behaviours between spherical and non-spherical particles are demonstrated from these studies. As the superficial gas velocity increases, the transition from a fixed bed to a bubbling fluidized bed happens smoothly for dry, spherical particles. However, the flow regime transitions are affected by particle shape (Kruggel-Emden & Vollmari Reference Kruggel-Emden and Vollmari2016). Channel flows can occur before bubbling flows for elongated rods and flat disks (Kruggel-Emden & Vollmari Reference Kruggel-Emden and Vollmari2016; Mahajan et al. Reference Mahajan, Padding, Nijssen, Buist and Kuipers2018b), caused by the stronger interlocking of non-spherical particles compared to spheres. The particle shape affects porosity and fluid drag force, consequently impacting pressure drop, minimum fluidization velocity and bed height (Zhong et al. Reference Zhong, Zhang, Jin and Zhang2009; Hilton, Mason & Cleary Reference Hilton, Mason and Cleary2010; Zhou et al. Reference Zhou, Pinson, Zou and Yu2011; Vollmari et al. Reference Vollmari, Oschmann, Wirtz and Kruggel-Emden2015; Gan, Zhou & Yu Reference Gan, Zhou and Yu2016; Vollmari, Jasevičius & Kruggel-Emden Reference Vollmari, Jasevičius and Kruggel-Emden2016; Ma, Xu & Zhao Reference Ma, Xu and Zhao2017; Mahajan et al. Reference Mahajan, Nijssen, Kuipers and Padding2018a). For example, consider three-dimensional (3-D) ellipsoidal particles in which two of the particle axes are equal (a = b); the aspect ratio (AR) of the ellipsoidal particle is defined as AR = c/a, where c is the length of the third axis. For fine ellipsoidal particles (Gan et al. Reference Gan, Zhou and Yu2016), the dependence of the minimum fluidization velocity, 
${U_{mf}}$, on the particle aspect ratio follows a W-shaped curve for 0 < AR < 3.5, with the peak occurring at AR = 1 (i.e. spheres). A V-shaped curve is observed for the dependence of the minimum bubbling velocity, 
${U_{mb}}$, on AR, with the minimum value of 
${U_{mb}}$ at AR = 1 (Gan et al. Reference Gan, Zhou and Yu2016). Larger bed expansions are obtained for more elongated particles or blockier particles than spherical particles, due to larger porosities (Zhong et al. Reference Zhong, Zhang, Jin and Zhang2009; Kruggel-Emden & Vollmari Reference Kruggel-Emden and Vollmari2016).
In gas-fluidized beds, elongated particles are aligned more horizontally in the lower part of the bed and more vertically in the regions close to the wall, and the particle orientational preference to align in the gas flow direction increases as the particle aspect ratio or the gas velocity increases (Cai et al. Reference Cai, Zhao, Li and Yuan2013; Oschmann, Hold & Kruggel-Emden Reference Oschmann, Hold and Kruggel-Emden2014; Cai, Yuan & Zhao Reference Cai, Yuan and Zhao2016; Ma et al. Reference Ma, Xu and Zhao2017; Mema et al. Reference Mema, Buist, Kuipers and Padding2020). Significant differences between spherical and elongated particles also occur in the particle velocity distributions and the particle circulation patterns (Mema et al. Reference Mema, Buist, Kuipers and Padding2020). Spheres have larger rotational velocities in dense zones than in the freeboard, caused by particle–particle interactions. In contrast, elongated particles rotate more prominently in the freeboard region than in dense zones, unhindered by the presence of other particles (Buist et al. Reference Buist, Jayaprakash, Kuipers, Deen and Padding2017). Particle mixing is also affected by the particle shape. The mixing rate decreases as particle elongation or particle blockiness increases (Ren et al. Reference Ren, Zhong, Jin, Shao and Yuan2013; Oschmann et al. Reference Oschmann, Hold and Kruggel-Emden2014). In fluidized mixtures of spherical and cylindrical particles, the elongated cylindrical particles sink to the lower part of the bed in the non-bubbling flow regime, while no segregation occurs in the bubbling flow regime due to bubble-induced mixing (Boer et al. Reference Boer, Buist, Deen, Padding and Kuipers2018). More asymmetrical bubbles, smaller bubble sizes and lower bubble rising velocities are obtained in fluidized beds of ellipsoidal particles than for beds comprised of spheres (Shrestha et al. Reference Shrestha, Kuang, Yu and Zhou2019), due to the interlocking and stronger particle–particle interactions occurring between ellipsoidal particles. Interparticle cohesion due to the van der Waals force increases the minimum fluidization and bubbling velocities (Gan et al. Reference Gan, Zhou and Yu2016) and reduces the bubble sizes and rising velocities (Shrestha et al. Reference Shrestha, Kuang, Yu and Zhou2020) in fluidized beds of ellipsoidal particles.
For a more comprehensive understanding of fluid–particle two-phase flow behaviour in fluidized beds, computational simulations based on the discrete element method (DEM) and computational fluid dynamics (CFD) are often performed, particularly because transient and particle-scale information can be obtained. In the DEM-CFD method, the fluid cell size is usually larger than the largest dimension of a particle (see figure 1). The average quantities in the fluid cell, such as velocity, density, pressure and porosity, are used to calculate the fluid–particle interaction force through a drag force model. The Hölzer–Sommerfeld drag force model (Hölzer & Sommerfeld Reference Hölzer and Sommerfeld2008) is usually employed for rod-like and disk-like particles, and the particle shape and orientation with respect to the fluid flow direction are taken into account in the model. Besides the drag force 
${\boldsymbol{F}_{drag}}$, a lift force 
${\boldsymbol{F}_{lift}}$ and hydrodynamic torque 
${\boldsymbol{T}_{hyd}}$ also exist, as shown in figure 1. However, the lift force and hydrodynamic torque have usually been neglected in previous dense-phase simulations (Zhong et al. Reference Zhong, Zhang, Jin and Zhang2009; Hilton et al. Reference Hilton, Mason and Cleary2010; Zhou et al. Reference Zhou, Pinson, Zou and Yu2011; Vollmari et al. Reference Vollmari, Oschmann, Wirtz and Kruggel-Emden2015, Reference Vollmari, Jasevičius and Kruggel-Emden2016; Gan et al. Reference Gan, Zhou and Yu2016; Ma et al. Reference Ma, Xu and Zhao2017; Mahajan et al. Reference Mahajan, Nijssen, Kuipers and Padding2018a; Shrestha et al. Reference Shrestha, Kuang, Yu and Zhou2019, Reference Shrestha, Kuang, Yu and Zhou2020), because reliable models of lift force and torque are not available and the effect of surrounding particles on the lift force and torque cannot be accurately considered for non-spherical particles in the DEM-CFD scheme. The effects of lift force and hydrodynamic torque on fluidized beds of elongated sphero-cylinders were investigated by Mema et al. (Reference Mema, Mahajan, Fitzgerald and Padding2019). It was found that smaller particle velocities were obtained without the lift force, and the inclusion of the hydrodynamic torque model led to more elongated particles oriented in the direction perpendicular to the gas flow (Mema et al. Reference Mema, Mahajan, Fitzgerald and Padding2019). Detailed information of the theoretical aspects and typical applications of the DEM and coupled DEM-CFD methods is provided in several review articles (Zhu et al. Reference Zhu, Zhou, Yang and Yu2007; Guo & Curtis Reference Guo and Curtis2015; Lu, Third & Müller Reference Lu, Third and Müller2015; Zhong et al. Reference Zhong, Yu, Liu, Tong and Zhang2016).
An illustration of drag force, lift force and hydrodynamic torque exerted on a non-spherical particle.
The previous studies of fluidization mostly focus on rigid particles that show no significant deformation. However, flexible, fibrous granular materials are common, such as biomass materials (e.g. plant stems), wool, hair, textile fabric and carbon fibres. These flexible granular materials normally have large particle aspect ratios and experience pronounced bending deformation, resulting in different flow behaviours from rigid particles (Guo et al. Reference Guo, Wassgren, Hancock, Ketterhagen and Curtis2015). Computational and experimental studies of deformable, cut tobacco leaves in a fluidized bed dryer have been conducted (Geng et al. Reference Geng, Luo, Li, Yuan, Yuan and Wu2014; Wu et al. Reference Wu, Gao, Yuan, Li, Zhu, Wu, Zhang, Lu and Luo2018, Reference Wu, Dai, Li, Gu, Yuan and Luo2019, Reference Wu, Zhang, Xu, Yuan, Zhu, Li, Wang and Luo2020), and the particle velocities and particle spatial distributions were examined. It was found that the particle spatial distribution was non-uniform in the fluidized bed riser and clusters of particles were more easily created at low gas velocities and high solids mass flow rates.
For a better understanding of the effects of particle flexibility on fluidization, dense gas-fluidized beds of flexible fibres are numerically studied in the present work using a coupled DEM-CFD method. A flexible fibre is modelled as a string of spheres connected by elastic bonds. The gas–fibre interaction force is calculated based on the interaction between the gas and each spherical element, rather than on the interaction between the gas and the whole fibre, as has been done in previous work (Zhong et al. Reference Zhong, Zhang, Jin and Zhang2009; Hilton et al. Reference Hilton, Mason and Cleary2010; Zhou et al. Reference Zhou, Pinson, Zou and Yu2011; Vollmari et al. Reference Vollmari, Oschmann, Wirtz and Kruggel-Emden2015, Reference Vollmari, Jasevičius and Kruggel-Emden2016; Gan et al. Reference Gan, Zhou and Yu2016; Ma et al. Reference Ma, Xu and Zhao2017; Mahajan et al. Reference Mahajan, Nijssen, Kuipers and Padding2018a; Shrestha et al. Reference Shrestha, Kuang, Yu and Zhou2019, Reference Shrestha, Kuang, Yu and Zhou2020). The present treatment of gas–fibre interaction allows the fluid cell size to be slightly larger than the diameter of the elemental spheres, but much smaller than the fibre length. As a result, gas–fibre two-phase flows with fibres of very large aspect ratio can be modelled using the present method. Simulations of fluidized beds with various fibre aspect ratios are performed. The predictions for pressure drop, minimum fluidization velocity, fibre orientation, mixing rate and bed height are analysed and compared with previous experimental and simulation results. By varying fibre bending moduli for fibres with an aspect ratio of AR = 6, the effects of fibre flexibility on the gas–fibre flows in fluidized beds is systematically investigated.
2. Coupled DEM-CFD method
Gas-fluidized beds of flexible fibres are simulated here using the coupled DEM-CFD method. The flexible fibres are modelled using a DEM-based bonded-sphere model, which was developed in previous work (Guo et al. Reference Guo, Wassgren, Hancock, Ketterhagen and Curtis2013). The gas flow is described using a traditional CFD method. Two-way coupling is considered for the gas–fibre interaction. The governing equations of the numerical approach used in this study are presented below.
2.1. Flexible fibre model
The flexible fibre model was developed and validated in previous work (Guo et al. Reference Guo, Wassgren, Hancock, Ketterhagen and Curtis2013). For convenience, the model is also summarized here. A fibre is formed by connecting a number of identical spheres using elastic bonds, as shown in figure 2(a). The bond forces/moments keep the constituent spheres connected within a fibre and resist fibre deformation subject to external forces. According to the bonded-particle model by Potyondy & Cundall (Reference Potyondy and Cundall2004), the bond forces and moments can be calculated incrementally as follows:
\begin{gather}\textrm{d}\boldsymbol{F}_n^b = K_n^b\,\textrm{d}\boldsymbol{\delta }_n^b = \frac{{{E_b}A}}{{{l_b}}}\,\textrm{d}\boldsymbol{\delta }_n^b = \frac{{{E_b}A}}{{{l_b}}}\boldsymbol{v}_n^r\,\textrm{d}t,\end{gather}
\begin{gather}\textrm{d}\boldsymbol{F}_t^b = K_t^b\,\textrm{d}\boldsymbol{\delta }_t^b = \frac{{{G_b}A}}{{{l_b}}}\,\textrm{d}\boldsymbol{\delta }_t^b = \frac{{{G_b}A}}{{{l_b}}}\boldsymbol{v}_t^r\,\textrm{d}t,\end{gather}
\begin{gather}\textrm{d}\boldsymbol{M}_n^b = K_{tor}^b\,\textrm{d}\boldsymbol{\theta }_n^b = \frac{{{G_b}{I_p}}}{{{l_b}}}\,\textrm{d}\boldsymbol{\theta }_n^b = \frac{{{G_b}{I_p}}}{{{l_b}}}\boldsymbol{\omega }_n^r\,\textrm{d}t,\end{gather}
\begin{gather}\textrm{d}\boldsymbol{M}_t^b = K_{ben}^b\,\textrm{d}\boldsymbol{\theta }_t^b = \frac{{{E_b}I}}{{{l_b}}}\,\textrm{d}\boldsymbol{\theta }_t^b = \frac{{{E_b}I}}{{{l_b}}}\boldsymbol{\omega }_t^r\,\textrm{d}t,\end{gather}
in which 
$\textrm{d}\boldsymbol{F}_n^b$, 
$\textrm{d}\boldsymbol{F}_t^b$, 
$\textrm{d}\boldsymbol{M}_n^b$ and 
$\textrm{d}\boldsymbol{M}_t^b$ are the incremental bond normal force, shear force, torsional moment and bending moment, respectively, and 
$\textrm{d}\boldsymbol{\delta }_n^b$, 
$\textrm{d}\boldsymbol{\delta }_t^b$, 
$\textrm{d}\boldsymbol{\theta }_n^b$ and 
$\textrm{d}\boldsymbol{\theta }_t^b$ are the incremental normal deformation, shear deformation, torsional angular deformation and bending angular deformation, respectively, of the bond. Thus, 
$K_n^b$, 
$K_t^b$, 
$K_{tor}^b$ and 
$K_{ben}^b$ represent the normal, shear, torsional and bending stiffnesses, respectively, of the bond. The cylindrical bond has a radius 
${r_b}$, a length 
${l_b}$, a cross-sectional area 
$A = {\rm \pi}r_b^2$, an area moment of inertia 
$I = {\rm \pi}r_b^4/4$ and a polar area moment of inertia 
${I_p} = {\rm \pi}r_b^4/2$. In the present model, the bond radius 
${r_b}$ is assumed to be the same as the constituent sphere radius 
${r_s}$. The bond Young's modulus 
${E_b}$ and shear modulus 
${G_b}$ are related through the Poisson ratio 
$\zeta $, i.e. 
${E_b} = 2(1 + \zeta ){G_b}$. The elastic bond has no mass. The incremental bond deformation can be obtained from the product of the corresponding relative translational/rotational velocity between two bonded spheres (
$\boldsymbol{v}_n^r$, 
$\boldsymbol{v}_t^r$, 
$\boldsymbol{\omega }_n^r$ or 
$\boldsymbol{\omega }_t^r$) and the time step 
$\textrm{d}t$.
Schematic diagrams for the hydrodynamic forces acting on a flexible fibre of AR = 6: (a) the hydrodynamic force is calculated for each constituent sphere based on the local mean variables of the gas; and (b) all the hydrodynamic forces are simplified into a resultant force 
${\boldsymbol{F}_{res}}$ and a resultant torque 
${\boldsymbol{T}_{res}}$ exerted on the centre of the mass of the whole fibre.
The fibre dynamics are determined by the collective motion of the constituent spheres comprising it, and the movement of individual constituent spheres is governed by Newton's second law of motion, i.e.
\begin{gather}{m_s}\frac{{{\textrm{d}^2}{\boldsymbol{x}_s}}}{{\textrm{d}{t^2}}} = {\boldsymbol{F}^c} + {\boldsymbol{F}^b} + {m_s}\boldsymbol{g} + \boldsymbol{F}_d^c + \boldsymbol{F}_d^b + {\boldsymbol{F}^{gs}},\end{gather}
\begin{gather}{J_s}\frac{{{\textrm{d}^2}{\boldsymbol{\theta }_s}}}{{\textrm{d}{t^2}}} = {\boldsymbol{M}^c} + {\boldsymbol{M}^b} + \boldsymbol{M}_d^c + \boldsymbol{M}_d^b,\end{gather}
in which 
${\boldsymbol{x}_s}$ and 
${\boldsymbol{\theta }_s}$ are the translational displacement vector and angular displacement vector, respectively, of the constituent sphere of mass 
${m_s}$ and moment of inertia 
${J_s} = {\textstyle{2 \over 5}}{m_s}r_s^2$. The translational motion of the sphere is driven by the contact force 
${\boldsymbol{F}^c}$ exerted by the spheres that come from a different fibre or are not bonded to it in the same fibre, the bond force 
${\boldsymbol{F}^b}$ exerted by the bonds connecting to it (the sum of 
$\boldsymbol{F}_n^b$ and 
$\boldsymbol{F}_t^b$), and the gravitational force 
${m_s}\boldsymbol{g}$. The moment 
${\boldsymbol{M}^c}$ arises from the tangential component of contact forces 
${\boldsymbol{F}^c}$, and 
${\boldsymbol{M}^b}$ represents the bond moment, including the moment induced by the bond shear force 
$\boldsymbol{F}_t^b$, bending moment 
$\boldsymbol{M}_t^b$ and twisting moment 
$\boldsymbol{M}_n^b$. A contact force 
${\boldsymbol{F}^c}$ can be decomposed into a normal force 
$\boldsymbol{F}_n^c$ and a tangential force 
$\boldsymbol{F}_t^c$. The normal force is computed based on the current normal overlap between two spheres using the Hertz model (Hertz Reference Hertz1882). The tangential force at static friction stage depends on the current normal force, tangential displacement and the loading path according to the Mindlin–Deresiewicz theory (Mindlin & Deresiewicz Reference Mindlin and Deresiewicz1953), and it is capped by the product of the friction coefficient 
$\mu $ and current normal contact force. Details of these contact force models and their implementations are provided in Thornton & Yin (Reference Thornton and Yin1991).
The contact damping force 
$\boldsymbol{F}_d^c$ (in (2.5)) is introduced to account for the energy dissipation when two spheres collide. The normal and tangential components of the contact damping force can be written as
\begin{gather}\boldsymbol{F}_{dn}^c ={-} 2\sqrt {{\textstyle{5 \over 6}}} {\kern 1pt} {\beta _c}\sqrt {m_f^\ast K_n^c} \boldsymbol{v}_n^c,\end{gather}
\begin{gather}\boldsymbol{F}_{dt}^c ={-} 2\sqrt {{\textstyle{5 \over 6}}} {\kern 1pt} {\beta _c}\sqrt {m_f^\ast K_t^c} \boldsymbol{v}_t^c,\end{gather}
in which 
$m_f^\ast $ is the effective mass of two contacting fibres of masses 
${m_{f1}}$ and 
${m_{f2}}$, 
$m_f^\ast{=} {m_{f1}}{m_{f2}}/({m_{f1}} + {m_{f2}})$. The contact normal stiffness 
$K_n^c$ and contact tangential stiffness 
$K_t^c$ have the expressions 
$K_n^c = \textrm{d}F_n^c/\textrm{d}{\delta _n}$ and 
$K_t^c = \textrm{d}F_t^c/\textrm{d}{\delta _t}$, in which 
${\delta _n}$ and 
${\delta _t}$ are contact normal displacement and tangential displacement, respectively. The variables 
$\boldsymbol{v}_n^c$ and 
$\boldsymbol{v}_t^c$ represent the normal and tangential components, respectively, of the relative velocity at the contact point. The negative sign indicates that the damping force direction is opposite to the relative velocity direction. The contact damping coefficient 
${\beta _c}$ determines the extent of energy dissipation during collisions. In (2.6), 
$\boldsymbol{M}_d^c$ is the moment arising from the tangential contact damping force 
$\boldsymbol{F}_{dt}^c$.
Analogous to contact damping, bond damping is proposed to take into account the energy dissipation due to elastic wave propagation within a fibre when the fibre deforms rapidly. Thus, bond damping forces 
$\boldsymbol{F}_d^b$ and damping moments 
$\boldsymbol{M}_d^b$ are exerted on the constituent spheres and they are assumed to be linearly proportional to the fibre deformation rates 
$\boldsymbol{v}_n^r$, 
$\boldsymbol{v}_t^r$, 
$\boldsymbol{\omega }_n^r$, and 
$\boldsymbol{\omega }_t^r$:
\begin{gather}\boldsymbol{F}_{dn}^b ={-} {\beta _b}\sqrt {2{m_s}K_n^b} \boldsymbol{v}_n^r,\end{gather}
\begin{gather}\boldsymbol{F}_{dt}^b ={-} {\beta _b}\sqrt {2{m_s}K_t^b} \boldsymbol{v}_t^r,\end{gather}
\begin{gather}\boldsymbol{M}_{dn}^b ={-} {\beta _b}\sqrt {2{J_s}\; K_{tor}^b} {\kern 1pt} \boldsymbol{\omega }_n^r,\end{gather}
\begin{gather}\boldsymbol{M}_{dt}^b ={-} {\beta _b}\sqrt {2{J_s}K_{ben}^b} {\kern 1pt} \boldsymbol{\omega }_t^r,\end{gather}
in which the bond damping coefficient 
${\beta _b}$ determines the energy dissipation rate. According to previous work (Guo et al. Reference Guo, Wassgren, Hancock, Ketterhagen and Curtis2015), the effective coefficient of restitution for a two-fibre collision depends on the contact damping coefficient 
${\beta _c}$, the bond damping coefficient 
${\beta _b}$ and the elastic modulus for fibre bending 
${E_b}$, and their correlations are analysed based on the simulation data in Guo et al. (Reference Guo, Wassgren, Hancock, Ketterhagen and Curtis2015).
In (2.5), 
${\boldsymbol{F}^{gs}}$ represents the interaction force between the gas and the constituent sphere, and it has the expression (Kafui, Thornton & Adams Reference Kafui, Thornton and Adams2002)
\begin{equation}{\boldsymbol{F}^{gs}} ={-} {V_s}{\bf \nabla }p + {V_s}{\bf \nabla\ \cdot }\ {\boldsymbol{\tau}_{\!g}} + \varepsilon {\boldsymbol{F}_{drag}},\end{equation}
in which 
${V_s}$ is the volume of a sphere, p is the local gas pressure, 
${\boldsymbol{\tau }_{\!g}}$ is the local viscous stress tensor of the gas, which has a linear relationship with shear rate (i.e. Newtonian fluid), 
$\varepsilon$ is the local porosity and 
${\boldsymbol{F}_{drag}}$ is the drag force exerted on the constituent sphere.
2.2. Equations for gas flow
The continuity and momentum equations for gas flow can be written as (Anderson & Jackson Reference Anderson and Jackson1967)
\begin{gather}\frac{{\partial (\varepsilon {\rho _g})}}{{\partial t}} + {\bf \nabla\ \cdot }\ (\varepsilon {\rho _g}{\boldsymbol{u}_g}) = 0,\end{gather}
\begin{gather}\frac{{\partial (\varepsilon {\rho _g}{\boldsymbol{u}_g})}}{{\partial t}} + {\bf \nabla\ \cdot }\ (\varepsilon {\rho _g}{\boldsymbol{u}_g}{\boldsymbol{u}_g}) ={-} {\bf \nabla }p + {\bf \nabla\ \cdot }\ {\boldsymbol{\tau }_{\!g}} - {\boldsymbol{F}_V} + \varepsilon {\rho _g}\boldsymbol{g},\end{gather}
where 
${\boldsymbol{u}_g}$ and 
${\rho _g}$ are the gas velocity and density, respectively. The gas–fibre interaction force per unit volume, 
${\boldsymbol{F}_V}$, is obtained by summing the gas forces 
${\boldsymbol{F}^{gs}}$ acting on all constituent spheres in a fluid cell, 
${n_s}$, and dividing by the volume of the fluid cell, 
${V_{cell}}$:
\begin{equation}{\boldsymbol{F}_V} = \frac{{\sum\limits_{i = 1}^{{n_s}} {{{({\boldsymbol{F}^{gs}})}_i}} }}{{{V_{cell}}}}.\end{equation}
According to Di Felice (Reference Di Felice1994), the drag force acting on a component sphere of diameter 
${d_s}$ can be expressed as
\begin{equation}{\boldsymbol{F}_{drag}} = \frac{1}{2}{C_D}{\rho _g}\frac{{{\rm \pi} d_s^2}}{4}{\varepsilon ^2}|{\boldsymbol{u}_g} - {\boldsymbol{v}_s}|({\boldsymbol{u}_g} - {\boldsymbol{v}_s}){\varepsilon ^{ - (\chi + 1)}},\end{equation}in which the drag force coefficient is written as
\begin{equation}{C_D} = {\left( {0.63 + \frac{{4.8}}{{R{e^{0.5}}}}} \right)^2},\end{equation}
and the sphere Reynolds number 
$Re$ is linearly proportional to the superficial slip velocity between the gas and sphere,
\begin{equation}Re = \frac{{{\rho _g}{d_s}\varepsilon |{\boldsymbol{u}_g} - {\boldsymbol{v}_s}|}}{{{\mu _g}}},\end{equation}
in which 
${\mu _\textrm{g}}$ is the gas shear viscosity. The porosity function 
${\varepsilon ^{ - ({\chi + 1} )}}$ is used to account for the effect of surrounding solids, with the parameter 
$\chi $ a function of the sphere Reynolds number,
\begin{equation}\chi = 3.7 - 0.65 \cdot \exp \left[ { - \frac{{{{(1.5 - {{\log }_{10}}Re)}^2}}}{2}} \right].\end{equation}A semi-implicit finite difference scheme employing a staggered grid is used for discretizing the Navier–Stokes equations on a 3-D Cartesian grid. The pressure and porosity scalars are defined at the centre of each fluid cell, while the fluid velocity components are defined at the cell faces. The fluid velocity inside a fluid cell is obtained by interpolation with the velocities at the cell faces.
In this work, as shown in figure 2(a), the interaction between a fibre and the gas is modelled by calculating the force between each constituent sphere and the gas. As illustrated in figure 2(b), the resultant gas force on a fibre, 
${\boldsymbol{F}_{res}}$, can be obtained by summing the gas forces exerted on all the constituent spheres in a fibre, and these gas–sphere interaction forces can automatically generate a resultant hydrodynamic torque, 
${\boldsymbol{T}_{res}}$, to rotate the fibre. Thus, the effect of hydrodynamic torque on the motion of fibres is considered in the present simulations. There is no need to calculate the resultant hydrodynamic force and torque in order to update the fibre dynamics, because the motion of a fibre is modelled by the collective motion of the constituent spheres, which just requires the gas forces acting on the spheres. Thus, the treatment of the gas–fibre interaction in the present work is different from in prior work (Zhong et al. Reference Zhong, Zhang, Jin and Zhang2009; Hilton et al. Reference Hilton, Mason and Cleary2010; Zhou et al. Reference Zhou, Pinson, Zou and Yu2011; Vollmari et al. Reference Vollmari, Oschmann, Wirtz and Kruggel-Emden2015, Reference Vollmari, Jasevičius and Kruggel-Emden2016; Gan et al. Reference Gan, Zhou and Yu2016; Ma et al. Reference Ma, Xu and Zhao2017; Mahajan et al., Reference Mahajan, Nijssen, Kuipers and Padding2018a; Shrestha et al. Reference Shrestha, Kuang, Yu and Zhou2019, Reference Shrestha, Kuang, Yu and Zhou2020), in which the resultant hydrodynamic force and torque were directly computed by assuming that the particle is immersed in a local homogeneous fluid field. In previous work, the fluid cell size is larger than the length of an elongated particle for the calculation of local mean fluid variables, as shown in figure 1. Therefore, the fluid field resolution can be poor if the particle length is large. Using the method proposed in this work, the fluid cell size is required to be larger than the constituent sphere size, but it can be much smaller than the fibre length (see figure 2). As a result, a smaller fluid cell size can be used in the present simulations than in the previous simulations mentioned above.
It should be noted that different drag forces are obtained for horizontally and vertically aligned fibres using the present numerical approach. As illustrated in figure 3, a fibre is composed of eight spherical elements and the fluid cell size is twice the spherical element size. For a fibre perpendicular to the streamwise direction (see figure 3a), the fluid velocity in each fluid cell is equal to the inflow velocity 
${U_0}$. For a fibre aligned in the streamwise direction (see figure 3b), the first fluid cell has velocity 
${U_0}$, and the fluid velocity of the second fluid cell is reduced to 
${U_1}$ (i.e. 
${U_1} \lt {U_0}$) due to the deceleration by the drag in the first fluid cell. In a similar fashion, the fluid velocity is gradually reduced in the third and fourth fluid cells along the fibre axis (i.e. 
${U_3} \lt {U_2} \lt {U_1} \lt {U_0}$). As a result, the resultant drag force on a horizontal fibre is greater than that on a vertical fibre. This treatment predicts the trend of drag force changing with the fibre orientation. The calculation of the drag forces in the present scheme, in which the local averaged fluid quantities in fluid cells are used, is less accurate than the calculation using the expensive direct numerical simulation (DNS) approach, in which fluid cell size is much smaller than the spherical element size and the fibre surface is resolved for the fluid–solid interaction. However, the calculation of the interaction forces and fluid flows using the present method is more accurate than using the scheme in which the fluid cell size is greater than the fibre length.
Illustration of fluid velocities in fluid cells for the horizontally and vertically aligned fibres: (a) horizontal fibre and (b) vertical fibre.
3. Numerical model and verification
A cubic domain is created for the simulations of gas-fluidized beds, as shown in figure 4(a). A gas inlet boundary with a constant upward velocity (i.e. superficial velocity 
$U$) and a constant pressure is set at the bottom, and a gas outlet boundary with zero gradients of gas velocity and pressure is assigned at the top. Periodic boundaries are specified in the horizontal directions (x and z directions) for the modelling of a small-scale fluidized bed without sidewall boundary effects. A number of fibres are randomly generated without contacts throughout the computational domain. Owing to the gravitational force, the fibres are deposited forming a densely packed bed in the lower region of the domain (see figure 4b). Fluidization is then performed by introducing a gas flow from the bottom boundary. When the superficial gas velocity is sufficiently large, the fibre bed is fluidized with significant fibre motion and bubble formation (see figure 4c).
(a) A sketch of the computational domain for the gas-fluidized bed simulations. (b) A deposited fibre bed with AR = 6 before fluidization and (c) a snapshot during the bubbling fluidization process from the simulation case 10; the orange arrows represent the local averaged gas velocity vectors. The fibres are coloured differently for visual convenience.
A number of simulation cases are performed, as presented in table 1. The base material properties of the fibres in the simulations are close to those of biomass materials, such as wheat straw. In the simulation cases, the fibre aspect ratios and bond bending moduli are varied for investigation of the effects of fibre elongation and flexibility on fluidization behaviour. All fibres have the same material density of 
${\rho _f} = 500\;\textrm{kg}\;{\textrm{m}^{ - 3}}$, the same friction coefficient of 
${\mu _f} = 0.485$, and Young's modulus for the calculation of fibre–fibre contact force 
${E_c} = 4.4 \times {10^9}\;\textrm{Pa}$. Assuming the two rigid fibres are in contact, the contact damping coefficient of 
${\beta _c} = 1.63 \times {10^{ - 2}}$ gives the effective coefficient of restitution 
${e_{eff}} = 0.95$, according to the correlation 
${\beta _c} ={-} \ln {e_{eff}}/\sqrt {{{\rm \pi} ^2} + {{\ln }^2}{e_{eff}}}$. The same contact damping coefficient of 
${\beta _c} = 1.63 \times {10^{ - 2}}$ is also used for the flexible fibres. The bond damping coefficient 
${\beta _b}$ is associated with the ratio of the relative velocities between two bonded spheres at two consecutive equilibrium positions 
${e_b}$ as 
${\beta _b} ={-} \ln {e_b}/\sqrt {{{\rm \pi} ^2} + {{\ln }^2}{e_b}} $. The value 
${\beta _b} = 3.35 \times {10^{ - 2}}$ is employed and leads to 
${e_b} = 0.9$ (Guo et al. Reference Guo, Wassgren, Hancock, Ketterhagen and Curtis2017). The granular flow dynamics is less sensitive to the values of 
${\beta _c}$ and 
${\beta _b}$ for the dense systems compared to the dilute systems, because a larger number of solids contacts in a dense flow can rapidly dissipate the kinetic energy. The present fluidization systems are regarded as dense systems and the values of 
${\beta _c}$ and 
${\beta _b}$ are expected to have limited impact on the macroscopic flow behaviours in the present work. The gas has an initial density of 
${\rho _g} = 1.204\;\textrm{kg}\;{\textrm{m}^{ - 3}}$ and pressure of 
$p = 1.01325 \times {10^5}\;\textrm{Pa}$. The temperature of the gas remains constant at a room temperature of 20 °C during the fluidization process. In coupled DEM-CFD simulations, the time step required for the DEM stability is much smaller than the time step for the CFD scheme, as discussed by Kafui et al. (Reference Kafui, Thornton and Adams2002). Thus, the critical time step is determined by the DEM simulations. According to Guo et al. (Reference Guo, Wassgren, Hancock, Ketterhagen and Curtis2013), the time step for the flexible fibre DEM simulations should be less than the time it takes for an axial extensional/compressional wave to travel a single bond length, and the critical time step is determined as
\begin{equation}\Delta {t_{cri}} = 0.8165 \cdot {l_b}\sqrt {{\rho _f}/{E_b}} .\end{equation}
The real time step used in the simulations is conservatively chosen as a small fraction (<0.5) of 
$\Delta {t_{cri}}$ to ensure numerical stability and accuracy.
Simulation cases considered in the present work.
Owing to the periodic boundary conditions, different computational domains are employed to examine the sensitivity to domain size. For the fibres with AR = 1 (spheres) and AR = 3, two different domains with the dimensions 
$\; l \times b \times h$ of 60 mm × 24 mm × 150 mm and 120 mm × 24 mm × 150 mm are simulated. For the fibres with AR = 4, three different domains of 60 mm × 24 mm × 150 mm, 60 mm × 48 mm × 150 mm and 120 mm × 24 mm × 150 mm are checked. Therefore, the domain size variations in the two periodic directions are examined. The pressure drop is measured as the gas pressure at the inlet boundary minus the gas pressure at the outlet boundary. As the superficial gas velocity increases, the pressure drop initially increases to a peak value and then remains nearly constant at a plateau when the packed bed is fully fluidized. The dimensionless pressure drop 
$\Delta {P^\ast }$, which is defined as the current pressure drop divided by the average pressure drop in the fluidized flow regime (the plateau value), is plotted as a function of superficial gas velocity U in figure 5 for the various fibres in different domains. It can be seen from figure 5 that the pressure drop results are insensitive to the variation of the domain sizes considered. Thus, it is believed that domain-size-independent results can be obtained with the domains employed in the present simulations.
Dimensionless pressure drop as a function of superficial gas velocity with different periodic domain sizes.
A comprehensive validation of the flexible fibre model has been made in previous work for a single fibre deformation (Guo et al. Reference Guo, Wassgren, Hancock, Ketterhagen and Curtis2013) and granular shear flows of assemblies of fibres (Guo et al. Reference Guo, Buettner, Lane, Wassgren, Ketterhagen, Hancock and Curtis2019). As a primary verification of the gas–fibre interaction model, a simple process, that is, gas flow through a static, densely packed bed of spheres (i.e. AR = 1 fibres in the present study), is examined. A shallow sphere bed (case 1) and a deep sphere bed (case 4) with double the fibre number of the shallow bed are considered in the simulations. As shown in figure 6, the pressure drops at various superficial gas velocities obtained from the simulations are in good agreement with the predictions of the Ergun equation (Ergun Reference Ergun1952) for both the shallow and deep beds, verifying the gas–sphere interaction model. To evaluate the present model for the interaction between gas and elongated fibres, the present simulation results are compared with previous numerical and experimental results in figures 7 and 8, which is discussed in the following section.
Comparison of present simulation results with the Ergun equation for the gas flowing through (a) a shallow particle bed and (b) a deep particle bed of spheres.
(a) Dimensionless pressure drop versus superficial velocity for various fibre aspect ratios. (b) Minimum fluidization velocity and initial average porosity as a function of fibre aspect ratio. (c) Dimensionless pressure drop versus dimensionless superficial gas velocity for various fibre aspect ratios.
(a) An illustration of fibre orientation angle, (b) time evolution of the fraction of fibres that are oriented in the streamwise direction (y direction) with angle 
$\theta $ in the specified ranges for AR = 3 fibres at 
$U/{U_{mf}} = 1.22$, (c) probability density distribution of orientation angle 
$\theta $ for fibres with AR = 3 at 
$U/{U_{mf}} = 1.22$, and (d) fraction of fibres with orientation angle 
$\theta $ between 0° and 45° as a function of superficial gas velocity.
4. Fibre elongation
The effect of fibre elongation is studied by considering different fibre aspect ratios (AR = 1, 2, 3, 4 and 6) in the simulations (cases 6–10 in table 1). The fibres of different aspect ratios have the same volume and thus the same equivalent diameter. For the present studies on fibre elongation effects, the fibres of bond bending modulus 
${E_b} = 4.4 \times {10^9}$ Pa exhibit very small deformation and behave like rigid rods in the fluidized beds. The results of pressure drop, minimum fluidization velocity, gas–fibre interaction force, bed height and mixing rate are analysed and compared with previous experimental and numerical observations.
4.1. Pressure drop and minimum fluidization velocity
The dimensionless pressure drop 
$\Delta {P^\ast }$ is plotted as a function of superficial gas velocity U in figure 7(a) for various fibre aspect ratios. The shape of the 
$\Delta {P^\ast }\text{--} U$ curves for the elongated fibres (AR > 1) is similar to that of the spherical particles (AR = 1). However, the fibre aspect ratio has a significant impact on the pressure drops before the bed is fluidized. At a given superficial velocity U before the bed is fluidized, the pressure drop decreases as the fibre aspect ratio increases for fibres with AR ≥ 2, causing the 
$\Delta {P^\ast }\text{--} U$ curve to move to the right. The pressure drop curve of the spherical particles (AR = 1) is located between the AR = 3 and AR = 4 fibres. The minimum fluidization velocity 
${U_{mf}}$, defined as the critical superficial gas velocity at which a packed bed becomes a fluidized bed (as illustrated in figure 7a), also exhibits a strong dependence on the fibre aspect ratio. As shown in figure 7(b), the 
${U_{mf}}$ decreases as AR increases from 1 to 2, then increases as AR increases from 2 to 6. Similar dependences of the pressure drop and minimum fluidization velocity on the fibre aspect ratio were observed in previous simulation studies of ellipsoidal particles (Zhou et al. Reference Zhou, Pinson, Zou and Yu2011) and rod-like particles (Ma et al. Reference Ma, Xu and Zhao2017). It is believed that the difference in the gas–fibre interaction forces, which are influenced by the porosity and particle aspect ratio, plays a critical role in the different fluidization behaviours for the fibres with different aspect ratios, which is discussed later.
By normalizing the superficial gas velocity using the minimum fluidization velocity, 
$U/{U_{mf}}$, the data can be collapsed for various fibre aspect ratios, as shown in figure 7(c). Thus, the dimensionless pressure drop is shown to be a function of the gas velocity ratio, 
$U/{U_{mf}}$, for fibres with various aspect ratios. The experimental results of a fluidized bed with rigid cylinders of AR = 3.5 by Vollmari et al. (Reference Vollmari, Jasevičius and Kruggel-Emden2016) and previous simulation results with AR = 4 cylinders by Ma et al. (Reference Ma, Xu and Zhao2017) fall on the same master curve as the present numerical simulation results. In the simulations by Ma et al. (Reference Ma, Xu and Zhao2017), the Hölzer–Sommerfeld model (Hölzer & Sommerfeld Reference Hölzer and Sommerfeld2008) is used for the gas–whole composite particle interaction, while in the present simulations, the Di Felice model (Di Felice Reference Di Felice1994) is used for the gas–constituent sphere interaction. Despite the difference in the treatment of the gas–fibre interaction forces, similar pressure drop results are obtained in the simulations by Ma et al. (Reference Ma, Xu and Zhao2017) and the present simulations. Thus, the good agreement between the present results and previous experimental and numerical results justifies the treatment of the gas–fibre interaction force in this work, although the local averaged fluid quantities in fluid cells (larger than spherical element size) are used to calculate the interaction forces without considering the detailed fluid flow on the surfaces of the fibres. In addition, the spheres and shorter fibres with 
$AR = 2$ show the fluidization behaviours of Geldart's group B particles, in that the bubbles are formed at the minimum fluidization velocities. However, for the longer fibres with 
$AR \ge 3$, vertical gas flow channels or fibre flow jets are formed in the bed before the bubbling flow, causing bubbling velocities larger than the minimum fluidization velocities. The occurrence of the channelling flow before bubbling flow for the large-aspect-ratio fibres was also observed in previous experiments (Kruggel-Emden & Vollmari Reference Kruggel-Emden and Vollmari2016; Mahajan et al. Reference Mahajan, Padding, Nijssen, Buist and Kuipers2018b).
4.2. Fibre orientation
The orientation of a stiff fibre (of bending modulus 
${E_b} = 4.4 \times {10^9}\;\textrm{Pa}$ in the present study) can be described by the angle, 
$\theta $, between the fibre axis (the line joining the centres of the two spherical elements at the ends of a fibre) and the streamwise direction (y direction), as shown in figure 8(a). The time evolution of the fraction of fibres with the orientation angle 
$\theta $ in specified ranges is shown in figure 8(b) for AR = 3 fibres at 
$U = 1.4\;\textrm{m}\;{\textrm{s}^{ - 1}}$. The fraction for each angle range fluctuates around the average value during the fluidization process, indicating that the probability distribution of the orientation angle does not change much when the bed is fluidized.
The average probability density distribution of the orientation angle 
$\theta $ is plotted in figure 8(c). The probability density increases with increasing 
$\theta $, indicating that more fibres tend to be horizontally aligned. It is also observed that the probability density distribution profile obtained from the present simulation is qualitatively consistent with the profile from the previous simulation by Ma et al. (Reference Ma, Xu and Zhao2017). However, a larger portion of horizontally aligned fibres are obtained in the present simulations compared to the simulations by Ma et al. (Reference Ma, Xu and Zhao2017). This discrepancy may be due to the effect of the hydrodynamic torque acting on the particles, which is considered in the present simulations while not taken into account in the simulations by Ma et al. (Reference Ma, Xu and Zhao2017). According to the work by Mema et al. (Reference Mema, Mahajan, Fitzgerald and Padding2019), the hydrodynamic torque induces more elongated particles to orient in the direction perpendicular to the gas flow, i.e. in the horizontal direction.
The average fraction of the fibres that are more aligned in the streamwise (vertical) direction, i.e. with 
$\theta $ between 0° and 45°, is plotted as a function of the superficial gas velocity in figure 8(d). For the relatively shorter fibres with AR = 3, a small portion of vertically aligned fibres exist in the fixed bed (
$U \lt {U_{mf}}$). A rapid increase in the fraction of the vertically aligned fibres is observed as the superficial gas velocity approaches 
${U_{mf}}$, and the fraction changes slightly when the superficial gas velocity is above 
${U_{mf}}$ for the AR = 3 fibres. In contrast, for the longer fibres with AR = 6, the fraction of vertical fibres still increases with the increasing superficial gas velocity even after 
${U_{mf}}$, which is consistent with previous experimental observation (Mema et al. Reference Mema, Buist, Kuipers and Padding2020). The slower increase to the asymptote of the fraction of vertical fibres for AR = 6 is attributed to a larger superficial gas velocity U required to create enough voidage in the bed for the more elongated fibres to re-orient.
4.3. Gas–fibre interaction force
Figure 9 shows the spatial distributions of the gas–fibre interaction force normalized by the fibre gravitational force, 
${F^{gf}}/{G_f}$, at different time instants during the fluidization process for AR = 6 (case 10 in table 1) at 
$U = 1.7\;\textrm{m}\;{\textrm{s}^{ - 1}}$ (
$U/{U_{mf}} = 1.172$). By plotting the time-averaged gas–fibre interaction force per fibre as a function of the vertical position in the bed (not shown here), it is found that larger gas–fibre interaction forces occur in the lower part of the bed than in the upper part. It is observed from figure 9 that larger gas–fibre interaction forces more likely occur below and above a gas bubble, which is due to the high velocities of the gas flowing into and out of the bubble (Shrestha et al. Reference Shrestha, Kuang, Yu and Zhou2019).
Distributions of gas–fibre interaction force normalized by the fibre gravitational force, 
${F^{gf}}/{G_f}$, during the fluidization process for AR = 6 (simulation case 10) at 
$U = 1.7\;\textrm{m}\;{\textrm{s}^{ - 1}}$ (
$U/{U_{mf}} = 1.172$).
Probability density distributions of the normalized gas–fibre interaction forces for various fibre aspect ratios are plotted in figure 10. At a superficial gas velocity of 
$U = 1\;\textrm{m}\;{\textrm{s}^{ - 1}}$ (below the minimum fluidization velocities of all the fibres considered), the mean gas–fibre interaction force decreases as the fibre aspect ratio AR increases from 2 to 6, and the mean interaction force of the spheres (AR = 1) is close to that of the fibres with AR = 4 (see figure 10a). As the fibre aspect ratio increases, the porosity of the packed bed increases (figure 7b). The increase in porosity contributes to a reduction in the interstitial gas velocity. On the other hand, the average projected area of fibres (on the horizontal plane, which is perpendicular to the streamwise direction) increases as AR increases because most of the elongated fibres tend to be horizontally aligned (with the orientational angle 
$\theta $ between 45° and 90°) in a fixed bed (
$U \lt {U_{mf}}$) as indicated in figure 8(d). As the fibre aspect ratio increases from 1 to 2, the effect of the increase in fibre projected area exceeds the effect of the increase in porosity and, therefore, the average gas–fibre interaction force increases. As the fibre aspect ratio increases from 2 to 6, the effect of the increase in porosity becomes more dominant, and the average gas–fibre interaction force decreases as the interstitial velocity decreases. It is believed that the gas–fibre interaction forces as shown in figure 10(a) lead to the effect of the fibre aspect ratio on the pressure drops and the minimum fluidization velocities as shown in figure 7(a,b). The smaller gas–fibre interaction forces at a superficial velocity below 
${U_{mf}}$ are responsible for the smaller pressure drops and the larger values of 
${U_{mf}}$.
Probability density distributions of the normalized gas–fibre interaction forces, 
${F^{gf}}/{G_f}$, at the superficial gas velocities (a) 
$U = 1.0\;\textrm{m}\;{\textrm{s}^{ - 1}}$, (b) 
$U = {U_{mf}}$ and (c) 
$U = 1.7\;\textrm{m}\;{\textrm{s}^{ - 1}}$ for various fibre aspect ratios, AR. (d) Probability density distribution of the gas–fibre interaction force, 
${F^{gf}}$, normalized by the average interaction force at the minimum fluidization velocity, 
$F_{mf}^{gf}$, for 
$U = 1.7\;\textrm{m}\;{\textrm{s}^{ - 1}}$.
Figure 10(b) shows the probability density distributions of the normalized gas–fibre interaction forces for various fibre aspect ratios at the corresponding minimum fluidization velocities. It can be seen that larger gas–fibre interaction forces are obtained for the larger fibre aspect ratios (AR = 4 and 6) at 
${U_{mf}}$, indicating that larger hydrodynamic forces are required to fluidize the more elongated fibres. Previous work (Guo et al. Reference Guo, Wassgren, Hancock, Ketterhagen and Curtis2015) shows that in dense fibre flows the shear stress and coordination number increase as the fibre aspect ratio increases. As a result, larger external forces should be exerted on the larger-aspect-ratio fibres to make them flow. At a superficial gas velocity above the minimum fluidization velocities, as shown in figure 10(c), the mean gas–fibre interaction force increases with an increase in the fibre aspect ratio. Thus, larger hydrodynamic forces are needed to maintain the fluidized state for the fibres with larger aspect ratios. Probability density distributions of the gas–fibre interaction force, 
${F^{gf}}$, normalized by the average interaction force at the minimum fluidization velocity, 
$F_{mf}^{gf}$, for 
$U = 1.7\;\textrm{m}\;{\textrm{s}^{ - 1}}$ are plotted in figure 10(d). It is interesting to see that the curves with various fibre aspect ratios tend to collapse, indicating that the quantities at the state of minimum fluidization provide a benchmark against which the results of fully fluidized beds are evaluated. In addition, a narrower distribution of the gas–fibre interaction forces is observed for the longer fibres with AR = 6. As shown in figure 9, significant variation of the interaction forces occurs surrounding the gas bubble. At the same superficial velocity 
$U = 1.7\;\textrm{m}\;{\textrm{s}^{ - 1}}$, the gas bubble size is smaller for the larger-aspect-ratio fibres, leading to a smaller variation of the gas–fibre interaction forces. Thus, the narrower probability density distribution of the gas–fibre interaction forces is obtained for the larger-aspect-ratio fibres due to the smaller bubble size.
4.4. Fibre bed height
During the fluidization process, a number of fibres rise up driven by the gas flow and then fall down due to the gravitational force. As a result, the total potential energy due to gravity in a fluidized bed of spheres (AR = 1) exhibits a periodic, sinusoidal evolution with time, as shown in figure 11(a). The total kinetic energy of the spheres also shows a periodic variation with time. However, two unequal peaks are observed in each period for the kinetic energy. The smaller peak in the kinetic energy, which occurs as the potential energy increases, is caused by the acceleration of spheres lifted by the gas flow. The larger peak, which occurs in the decreasing portion of the potential energy, is due to the conversion of potential energy to kinetic energy as the particles fall down. In the fluidization of elongated fibres with AR = 6, the oscillatory evolution of the kinetic and potential energies is also observed, as shown in figure 11(b). However, the magnitudes of the energies vary significantly from period to period for elongated fibres. This energy variation may be due to the intermittent interlocking of elongated fibres, which inhibits fibre motion. The fibres can gain larger energies again once the interlocking structure is unlocked. Owing to stronger fibre–fibre interaction, elongated fibres have a much larger period of energy change (T = 0.4 s) than spheres (T = 0.13 s). It should be noted that the strong periodicity in the energy data is likely due to the small computational domain size, which results in approximately one bubble at a time. A larger computational domain size would result in multiple, simultaneous bubbles and result in less periodicity in the energy plots.
Time evolution of the kinetic energy and gravitational potential energy of the fibres with (a) AR = 1 (spheres) and (b) AR = 6 during the fluidization processes at 
$U/{U_{mf}} = 1.17$.
The height of a fluidized bed, H, is described by
\begin{equation}H = \frac{{2\sum\limits_{i = 1}^{{N_f}} {{y_i}} \; }}{{{N_f}}},\end{equation}
in which 
${y_i}$ is the vertical coordinate of the centre of mass of fibre i and 
${N_f}$ is the total number of fibres. In the fluidization process, the bed height fluctuates periodically. Thus, an average bed height during at least five fluctuating periods, H, is used to characterize the height of a fluidized bed at a specified superficial gas velocity. The bed expansion ratio, 
$H/{H_{mf}}$, is plotted as a function of the normalized superficial gas velocity, 
$U/{U_{mf}}$, in figure 12, in which 
${H_{mf}}$ represents the average bed height at the minimum fluidization velocity. For the shorter fibres with AR ≤ 4, the bed shows an expansion before the bed is fluidized (
$U/{U_{mf}} \lt 1$). While for the longer fibres with AR = 6, the bed height remains unchanged before the fluidization, which is due to the large porosity of the packed bed. When the bed is fluidized (
$U/{U_{mf}} \gt 1$), the bed height increases with increasing superficial gas velocity due to the increase in the hydrodynamic forces acting on the fibres. The bed expansion ratio decreases as AR increases from 1 to 2, and the expansion ratio increases as AR increases from 2 to 6. The effect of fibre aspect ratio on the bed expansion is similar to its effect on the minimum fluidization velocity (see figure 7b). Thus, a larger minimum fluidization velocity corresponds to a larger bed expansion ratio at a given value of 
$U/{U_{mf}}$.
Bed expansion ratio, 
$H/{H_{mf}}$, as a function of normalized superficial gas velocity, 
$U/{U_{mf}}$, for various fibre aspect ratios.
4.5. Mixing rate
Solids mixing rate is critical for the operation of fluidized beds, such as in coating and chemical reactions. In the present work, the solids mixing behaviour in the flow (i.e. vertical) direction is investigated. The fibre bed at any time instant can be partitioned by a horizontal plane into an upper region and a lower region that both contain the same number of fibres. All the centres of mass of the fibres in the upper and lower regions are, respectively, above and below the horizontal plane. Thus, the vertical position of this horizontal plane changes with time as the bed fluctuates. The fibres in the upper region at time 
$t + \Delta t$ are compared with the fibres in the upper region at time t, and the number of new fibres that move into the upper region from the lower region during the time period 
$\Delta t$ is recorded as 
$\Delta N_{upper}^{new}$. The mixing rate, 
${R_{mix}}$, is defined as
\begin{equation}{R_{mix}} = \frac{{\Delta N_{upper}^{new}}}{{N_{upper}^{tot} \cdot \Delta t}} \times 100\,\%,\end{equation}
in which 
$N_{upper}^{tot}$ is the total number of fibres in the upper region, which is equal to half of the total number of fibres in the bed. The time evolution of the mixing rate, 
${R_{mix}}$, for the fibres of AR = 6 (case 10 in table 1) at 
$U = 1.9\;\textrm{m}\;{\textrm{s}^{ - 1}}$ (
$U/{U_{mf}} = 1.31$) is shown in figure 13. A pronounced periodic fluctuation of the mixing rate is observed. The largest mixing rates occur when the fibre bed is expanded and the air bubbles are about to burst on the top surface, while the mixing rates are smaller for the more packed beds, as shown by the inserted images in figure 13.
Time evolution of the mixing rate, 
${R_{mix}}$, for the fibres of AR = 6 (case 10 in table 1) at 
$U = 1.9\;\textrm{m}\;{\textrm{s}^{ - 1}}$ (
$U/{U_{mf}} = 1.31$). The inserts show snapshots of the fibre flow pattern at the corresponding time instants.
The mean mixing rate, 
${R_{mix}}$, can be determined for a specified superficial velocity as a time-averaged value of the mixing rates during a sufficiently long period of fluidization. Figure 14(a) shows the mean mixing rate 
${R_{mix}}$ varying with the superficial gas velocity for various fibre aspect ratios. In general, the mean mixing rate exhibits a linear increase with the increase in the superficial gas velocity. Nevertheless, the mean mixing rate can level off or decrease at very high superficial gas velocities. It is observed that the flow regime transitions from bubbling fluidization to slugging fluidization when 
$U/{U_{mf}} \gt 1.5$. In slugging fluidization, big air bubbles, which fill most of the column cross-section, lift the fibres periodically. The plug-like movement of fibres in slugging fluidization prevents convective mixing of the fibres. At a given superficial gas velocity, the mean mixing rate decreases as the fibre aspect ratio increases from 2 to 6, and the mean mixing rate of the spheres (AR = 1) is between those of the fibres with AR = 3 and AR = 4. Thus, the effect of the fibre aspect ratio on the mixing rate is similar to the effect of the fibre aspect ratio on the minimum fluidization velocity. By normalizing the superficial gas velocity using the minimum fluidization velocity, 
$U/{U_{mf}}$, the mixing rate data collapse for various fibre aspect ratios, as shown in figure 14(b).
Mean mixing rate versus (a) superficial gas velocity and (b) normalized superficial gas velocity for various fibre aspect ratios.
5. Fibre flexibility
Fibre flexibility has a significant impact on the packing density (Guo et al. Reference Guo, Buettner, Lane, Wassgren, Ketterhagen, Hancock and Curtis2019) and the stresses in volume-fixed fibre flows (Guo et al. Reference Guo, Wassgren, Hancock, Ketterhagen and Curtis2015). In the present work, the effect of fibre flexibility on fluidized beds is investigated by performing a series of simulations with the fibre aspect ratio fixed at AR = 6 and the fibre bond bending modulus varying between 
${E_b} = 4.4 \times {10^2}\;\textrm{Pa}\; $ and 
${E_b} = 4.4 \times {10^9}\;\textrm{Pa}$ (cases 10–16 in table 1), leading to the bending flexibility term 
${E_b}I$ (
$I = {\rm \pi}r_b^4/4$ is the area moment of inertia) in the range between 
$4.9 \times {10^{ - 10}}\;\textrm{N}\;{\textrm{m}^2}$ and 
$4.9 \times {10^{ - 3}}\;\textrm{N}\;{\textrm{m}^2}$ (i.e. 
${E_b}I \in [4.9 \times {10^{ - 10}},\; \; 4.9 \times {10^{ - 3}}]\;\textrm{N}\;{\textrm{m}^2}$). The commonly encountered biomass materials and rubber/plastic wires have the values of 
${E_b}I \in [1.0 \times {10^{ - 5}},1.0 \times {10^{ - 2}}]\;\textrm{N}\;{\textrm{m}^2}$ (Pusca, Bobancu & Duta Reference Pusca, Bobancu and Duta2010; Guo et al. Reference Guo, Buettner, Lane, Wassgren, Ketterhagen, Hancock and Curtis2019), braided textile ropes have 
${E_b}I \in [1.0 \times {10^{ - 7}},\; \; 1.0 \times {10^{ - 4}}]\;\textrm{N}\;{\textrm{m}^2}$ and wool and hair have 
${E_b}I \in [1.0 \times {10^{ - 10}},\; \; 1.0 \times {10^{ - 6}}]\;\textrm{N}\;{\textrm{m}^2}$ (Żak & Kobielarz Reference Żak and Kobielarz2010).
5.1. Pressure drop and minimum fluidization velocity
Images of the packed beds of softer fibres with 
${E_b} = 4.4 \times {10^2}\;\textrm{Pa}$ and stiffer fibres with 
${E_b} = 4.4 \times {10^9}\;\textrm{Pa}\; $ are shown in figure 15(a,b), respectively. It can be seen that the bed of softer fibres is more densely compacted, with a smaller bed height than the bed of stiffer fibres. The more flexible fibres are easier to bend to fill the voids among the fibres, leading to a denser packing. As a result, the porosity generally decreases as the fibre bending modulus decreases, as illustrated in figure 15(c). The porosity reaches an upper limit when the fibre bending modulus is sufficiently large (
${E_b} \ge 4.4 \times {10^6}\;\textrm{Pa}$).
Snapshots of the initial fibre bed with the fibre bending moduli for fibres with AR = 6: (a) 
${E_b} = \; 4.4 \times {10^2}\;\textrm{Pa}$ and (b) 
${E_b} = \; 4.4 \times {10^9}\;\textrm{Pa}$. (c) The variation of the average porosity of the initial fibre bed with the fibre bending modulus, 
${E_b}$.
The dimensionless pressure drop 
$\Delta {P^\ast }$ is plotted as a function of superficial gas velocity U for various fibre bending moduli 
${E_b}$ in figure 16(a). The same 
$\Delta {P^\ast }\text{--} U$ correlation is obtained for fibre beds with the larger bending moduli (
${E_b} \ge 4.4 \times {10^6}\;\textrm{Pa}$), which have similar porosities. The 
$\Delta {P^\ast }\text{--} U$ curve moves to the left side as the fibre bending modulus 
${E_b}$ decreases from 
$4.4 \times {10^6}\;\textrm{Pa}$. Thus, the minimum fluidization velocity 
${U_{mf}}$ remains constant for 
${E_b} \ge 4.4 \times {10^6}\;\textrm{Pa}$ and it decreases as the bending modulus 
${E_b}$ decreases for 
${E_b} \lt 4.4 \times {10^6}\;\textrm{Pa}$, as shown in figure 16(b). The dimensionless pressure drop 
$\Delta {P^\ast }$ is plotted as a function of normalized superficial gas velocity, 
$U/{U_{mf}}$, in figure 16(c). The data of 
$\Delta {P^\ast }$ versus 
$U/{U_{mf}}$ collapse onto the same master curve for stiffer fibres with 
${E_b} \ge 4.4 \times {10^4}\;\textrm{Pa}$. When 
$U/{U_{mf}}$ is less than one (i.e. the bed is not fluidized), concave 
$\Delta {P^\ast }$ versus 
$U/{U_{mf}}$ curves are obtained for stiffer fibres with 
${E_b} \ge 4.4 \times {10^4}\;\textrm{Pa}$ and convex curves are obtained for softer fibres with 
${E_b} \lt 4.4 \times {10^4}\;\textrm{Pa}$. At a given velocity ratio 
$U/{U_{mf}}$ less than one, larger dimensionless pressure drops 
$\Delta {P^\ast }$ are obtained for the smaller bending moduli when 
${E_b} \lt 4.4 \times {10^4}\;\textrm{Pa}$. It is believed that the effect of the fibre bending modulus observed in figure 16 is caused by the difference in porosity as shown in figure 15(c), which influences the gas–fibre interaction forces as discussed in the following section.
(a) Dimensionless pressure drop versus superficial gas velocity for various fibre bending moduli. (b) Minimum fluidization velocity as a function of fibre bending modulus, 
${E_b}$. (c) Dimensionless pressure drop versus normalized superficial gas velocity for various fibre bending moduli. The fibre aspect ratio is AR = 6.
5.2. Gas–fibre interaction force
Probability density distributions of the normalized gas–fibre interaction forces, 
${F^{gf}}/{G_f}$, for various fibre bending moduli, 
${E_b}$, are plotted in figure 17. At the superficial gas velocity of 
$U = 1.0\;\textrm{m}\;{\textrm{s}^{ - 1}}$ (lower than all the minimum fluidization velocities considered), similar distributions are obtained for the stiffer fibres, those with fibre bending moduli 
${E_b} \ge 4.4 \times {10^6}\;\textrm{Pa}$. The mean gas–fibre interaction force increases as the fibre bending modulus decreases for the more flexible fibres with 
${E_b} \lt 4.4 \times {10^6}\;\textrm{Pa}$ due to the reduction in the bed porosity with decreasing fibre bending modulus (see figure 15c). As a result, larger pressure drops and smaller minimum fluidization velocities are obtained for the more flexible fibres with 
${E_b} \lt 4.4 \times {10^6}\;\textrm{Pa}$ before the bed is fluidized (figure 16a,b).
Probability density distributions of the normalized gas–fibre interaction forces, 
${F^{gf}}/{G_f}$, at the superficial gas velocities (a) 
$U = 1.0\;\textrm{m}\;{\textrm{s}^{ - 1}}$, (b) 
$U = {U_{mf}}$ and (c) 
$U = 1.7\;\textrm{m}\;{\textrm{s}^{ - 1}}$ for various fibre bending moduli, 
${E_b}$. The fibre aspect ratio AR is 6.
At the minimum fluidization velocities (
$U = {U_{mf}}$) or above (
$U = 1.7\;\textrm{m}\;{\textrm{s}^{ - 1}}$), as shown in figure 17(b,c), smaller mean gas–fibre interaction forces are obtained for the fibres with the smaller bending moduli (
${E_b} \lt 4.4 \times {10^6}\;\textrm{Pa}$). Thus, smaller hydrodynamic forces are required to fluidize the more flexible fibres. According to the previous work by Guo et al. (Reference Guo, Wassgren, Hancock, Ketterhagen and Curtis2015), the fibre–fibre interaction is weaker and the shear stress is smaller in flows of more flexible fibres. Thus, the more flexible fibres can flow subject to smaller hydrodynamic forces.
5.3. Fibre deformation
The extent of fibre bending deformation was quantified by a deforming parameter in previous work (Guo et al. Reference Guo, Wassgren, Hancock, Ketterhagen and Curtis2015). For fibre i the deforming parameter 
${S_i}$ is defined as
\begin{equation}{S_i} = \frac{1}{{{\phi _{max}}}}\sqrt {\frac{1}{K}\mathop \sum \limits_{j = 1}^K \phi _{ij}^2} ,\end{equation}
where 
${\phi _{ij}}$ is the angle between the two connected bonds as shown in figure 18, K is the number of angles 
${\phi _{ij}}$ in the fibre i, and 
${\phi _{max}}$ is the maximum angle that 
${\phi _{ij}}$ can achieve, which is equal to 
$2{\rm \pi} /3$ (see figure 18).
Schematic illustration for the definition of fibre bending deformation parameter S (adapted from Guo et al. Reference Guo, Wassgren, Hancock, Ketterhagen and Curtis2015).
To evaluate the overall fibre deformation for a bed of 
${N_f}$ fibres, the mean deforming parameter is employed:
\begin{equation}{S_{mean}} = \frac{1}{{{N_f}}}\mathop \sum \limits_{i = 1}^{{N_f}} {S_i}.\end{equation}
The mean deforming parameter 
${S_{mean}}$ has a value between zero and one. A larger value of 
${S_{mean}}$ reflects a larger degree of fibre bending deformation. As shown in figure 19(a), the mean deforming parameter 
${S_{mean}}$ does not change much with the superficial gas velocity for the more flexible fibres with 
${E_b} \le 4.4 \times {10^6}\;\textrm{Pa}$, and 
${S_{mean}}$ increases only slightly with the superficial gas velocity for the stiffer fibres with 
${E_b} \gt 4.4 \times {10^6}\;\textrm{Pa}$. During the fluidization process, the bending deformation of the fibres increases as the fibre bending modulus decreases, as shown in figure 19(b).
Mean fibre bending deformation parameter, Smean, variation with (a) normalized superficial gas velocity and (b) fibre bending modulus, 
${E_b}$, at the minimum fluidization velocity, 
${U_{mf}}$. The fibre aspect ratio is AR = 6.
5.4. Fibre bed height
The bed expansion ratio, 
$H/{H_{mf}}$, as a function of normalized superficial gas velocity, 
$U/{U_{mf}}$, is plotted in figure 20 for the various fibre bending moduli, 
${E_b}$. Before the beds are fluidized (
$U/{U_{mf}} \lt 1$), bed expansion is observed for the beds of more flexible fibres with 
${E_b} \le 4.4 \times {10^4}\;\textrm{Pa}$, which have smaller porosities. In contrast, the bed height remains constant for the stiffer fibres with 
${E_b} \ge 4.4 \times {10^6}\;\textrm{Pa}$ before the bed is fluidized. The significant bed expansion before the minimum fluidization velocity is a distinct feature for the very flexible fibres. The increase in the average porosity, due to the fibre bed expansion, causes a reduction in the pressure drop compared to a fixed bed without the bed height change. As a result, the shape of the pressure drop curves deviates from the Ergun relationship (Ergun Reference Ergun1952) for the more flexible fibres (
${E_b} = \; 4.4 \times {10^2}\;\textrm{Pa}$, 
$4.4 \times {10^3}\;\textrm{Pa}$ and 
$4.4 \times {10^4}\;\textrm{Pa}$), as shown in figure 16(a). When the packed beds are fluidized (
$U/{U_{mf}} \gt 1$), smaller bed expansion ratios are obtained for the more flexible fibres with 
${E_b} \le 4.4 \times {10^4}\;\textrm{Pa}$, due to the smaller porosities.
Bed expansion ratio as a function of normalized superficial gas velocity for various fibre bending moduli, 
${E_b}$. The fibre aspect ratio is AR = 6.
5.5. Mixing rate
The effect of fibre flexibility on the mixing rate, 
${R_{mix}}$, is shown in figure 21. A roughly linear increase in the mixing rate with increasing superficial gas velocity, U, is observed for each fibre bending modulus 
${E_b}$, as shown in figure 21(a). The 
${R_{mix}}\text{--} U$ curves move to the left as the fibres become more flexible (i.e. as 
${E_b}$ decreases), due to the smaller minimum fluidization velocities (see figure 16b). The reduction of the mixing rate at high superficial gas velocity for the fibres with 
${E_b} = 4.4 \times {10^2}\;\textrm{Pa}$ is due to the transition from the bubbling flow regime to the slugging flow regime, as discussed in § 4.5. By plotting the mean mixing rate against the normalized superficial gas velocity, 
$U/{U_{mf}}$, the mixing rate data collapse, as shown in figure 21(b). Based on figures 14 and 21, the mixing rate is essentially determined by the normalized superficial gas velocity, 
$U/{U_{mf}}$. The fibre aspect ratio and fibre flexibility influence the mixing rate by affecting the minimum fluidization velocity, 
${U_{mf}}$.
Variation of the mean mixing rate with (a) superficial gas velocity and (b) normalized superficial gas velocity for the various fibre bending moduli, 
${E_b}$. The fibre aspect ratio is AR = 6.
6. Conclusion
A DEM-based flexible fibre model is coupled with CFD to simulate fluidized beds of flexible fibres. In the present DEM-CFD method, a fibre is represented by a string of spheres that are connected by elastic bonds. The interaction between a fibre and the gas is modelled by calculating the interaction force between each constituent sphere and the gas. The present treatment of the fibre–gas interaction is different from those in previous work (Zhong et al. Reference Zhong, Zhang, Jin and Zhang2009; Hilton et al. Reference Hilton, Mason and Cleary2010; Zhou et al. Reference Zhou, Pinson, Zou and Yu2011; Vollmari et al. Reference Vollmari, Oschmann, Wirtz and Kruggel-Emden2015, Reference Vollmari, Jasevičius and Kruggel-Emden2016; Gan et al. Reference Gan, Zhou and Yu2016; Ma et al. Reference Ma, Xu and Zhao2017; Mahajan et al., Reference Mahajan, Nijssen, Kuipers and Padding2018a; Shrestha et al. Reference Shrestha, Kuang, Yu and Zhou2019, Reference Shrestha, Kuang, Yu and Zhou2020) in which the resultant interaction force between the elongated particle and the gas is directly calculated and acts on the centre of the mass of the particle. The fluid cell size needs to be larger than the fibre length in the previous work, while the fluid cell size can be much smaller than the fibre length in the present studies. The proposed DEM-CFD method is verified, as the present simulation results of fluidized beds are in good agreement with previous theories, experimental results and simulation results.
In packed beds, most of the elongated fibres tend to be aligned horizontally. When the packed beds are fluidized, ~30 % of the fibres (with 3 ≤ AR ≤ 6) are aligned in the vertical direction due to the gas flow-induced alignment of the elongated fibres in the streamwise direction. For fibres of the same volume, as the fibre aspect ratio increases, the average projected area of a fibre on the horizontal plane (perpendicular to the streamwise direction) increases in the packed beds, which contributes to the increase in the gas–fibre interaction forces. On the other hand, as the fibre aspect ratio increases, the porosity of the packed bed increases, which contributes to the decrease in the gas–fibre interaction forces. When the superficial gas velocity is smaller than the minimum fluidization velocity, the mean gas–fibre interaction force increases as the fibre aspect ratio increases from 1 to 2, for which the increase in the fibre projected area is more dominant. The mean gas–fibre interaction force decreases as the fibre aspect ratio continues to rise from 2 to 6, for which the increase in the bed porosity plays a bigger role. As a result, the minimum fluidization velocity decreases as the fibre aspect ratio increases from 1 to 2 and it increases as the fibre aspect ratio increases from 2 to 6. When the bed is fluidized (
$U \gt {U_{mf}}$), the mean gas–fibre interaction force generally increases with increasing fibre aspect ratio, indicating that larger hydrodynamic forces are required to fluidize the longer fibres, which have stronger fibre–fibre interactions. Larger bed expansion is obtained for the larger-aspect-ratio fibres at a given normalized superficial gas velocity 
$U/{U_{mf}}$ due to larger porosities and larger minimum fluidization velocities. The curve of the mixing rate versus superficial gas velocity moves to the left as the fibre aspect ratio increases from 1 to 2, and the curve moves to the right as the fibre aspect ratio increases from 2 to 6. Nevertheless, the curves of the mixing rate versus the normalized superficial gas velocity 
$U/{U_{mf}}$ collapse for various fibre aspect ratios.
The fibre flexibility has an impact on the fluidization behaviour. Larger fibre deformation occurs for more flexible fibres in fluidized beds. The porosity of a packed fibre bed increases as the fibre bending modulus 
${E_b}$ increases before it reaches a plateau value. Thus, before the bed is fluidized, larger gas–fibre interaction forces are obtained for the more flexible fibres due to the smaller porosities. As a result, the pressure drop–superficial gas velocity curve moves to the left and the minimum fluidization velocity decreases as the fibres become more flexible. In the plot of dimensionless pressure drop versus the normalized superficial gas velocity 
$U/{U_{mf}}$, the curves collapse for the stiffer fibres (
${E_b} \ge 4.4 \times {10^4}\;\textrm{Pa}$), while larger dimensionless pressure drops are observed for the more flexible fibres with 
${E_b} = 4.4 \times {10^2}\;\textrm{Pa}$ and 
$4.4 \times {10^3}\;\textrm{Pa}$ before the bed is fluidized (
$U/{U_{mf}} \lt 1$). At the minimum fluidization velocity or above, smaller gas–fibre interaction forces are required to fluidize the beds of more flexible fibres, due to the weaker fibre–fibre interactions. The fibre bed expansion occurs for the more flexible fibres before fluidization. When the bed is fluidized, the bed expansion ratio is smaller for the more flexible fibres, due to the smaller porosities. The mixing rate–superficial gas velocity curve moves to the left as the fibres become more flexible and the data collapse in the plot of the mixing rate versus the normalized superficial gas velocity 
$U/{U_{mf}}$. As a result, the mixing rate is essentially determined by the velocity ratio, 
$U/{U_{mf}}$. The effects of fibre aspect ratio and fibre flexibility on the mixing rate are correlated to their effects on the minimum fluidization velocity, 
${U_{mf}}$.
In the future, experimental work using the positron emission particle tracking (PEPT) technique (Parker Reference Parker2017) can be conducted to track the trajectory of a single fibre in the fluidization system. Thus, the drag force model can be assessed by comparing the individual fibre motion between the simulations and the experiments. Gas flow passing a single, flexible fibre can be simulated using the DNS approach in order to develop a more accurate gas–fibre interaction model by taking into account the effect of fibre deformation. In addition, DEM-DNS simulations and experimental work can be performed on the gas flow through fixed fibre beds of various porosities, so that the gas–fibre interaction model can be further advanced by more accurately considering the effect of the presence of surrounding fibres.
Acknowledgement
The National Science Foundation of China (Grants 11872333 and 91852205 for Y.G., Grant 11632016 for J.L., and Grant 91752117 for Z.Y.) and the Zhejiang Provincial Natural Science Foundation of China (Grant LR19A020001 for Y.G.) are acknowledged for their financial support.
Declaration of interests
The authors report no conflict of interest.
