1. Introduction
Fluidised granular materials find extensive use in industrial processes such as drying, mixing, adsorption, pyrolysis, catalysis and chemical conversion (Jain, Tschisgale & Fröhlich Reference Jain, Tschisgale and Fröhlich2021; Schwarzmeier et al. Reference Schwarzmeier, Rettinger, Kemmler, Plewinski, Núñez-González, Köstler, Rüde and Vowinckel2023; Wang et al. Reference Wang, Dai, Li, Lin, Patrascu and Gao2024a , Reference Wang, Ma, Li, Dai, Patrascu and Gaob ). Particles are often treated as spherical shapes in theoretical models, while most granular materials, like sand, catalysts, biomass particles and plastic particles, are non-spherical (Wang & Shen Reference Wang and Shen2022; Zhao, Zhao & Luding Reference Zhao, Zhao and Luding2023; Giahi & Bergstrom Reference Giahi and Bergstrom2024; Qi et al. Reference Qi, Chen, Fei, Zhou, Liu, Wu, Zhao and Yang2025). These non-spherical particles exhibit significantly different characteristics compared with their spherical counterparts in terms of fluid dynamics, heat transfer, mass transfer and chemical reactions (Gao et al. Reference Gao, Yu, Li, Panday, Xu, Li, Ashfaq, Hughes and Rogers2020; Su et al. Reference Su, Ma, Yan, He and Xu2023; Eghbalmanesh et al. Reference Eghbalmanesh, Fathiganjehlou, Rieder, Hoogendoorn, Miloshevska, Baltussen, Peters, Buist and Kuipers2024; Jiang, Xu & Zhao Reference Jiang, Xu and Zhao2024; Fan et al. Reference Fan, Liu, Heynderickx, Deivendran, Guo and Chen2025). Therefore, a comprehensive understanding of the behaviour of non-spherical particles in fluidised systems is crucial for the effective design, optimisation and scaling of industrial applications (Gao et al. Reference Gao, Yu, Lu and Rogers2021b ; Ma et al. Reference Ma, Zhou, Liu, Chen, Xia and Zhao2022; Sanjeevi, Dietiker & Padding Reference Sanjeevi, Dietiker and Padding2022; Liu et al. Reference Liu, Ai, Wang, Liu, Shuai and Pan2024). However, there is a notable scarcity of research focusing on non-spherical particles in fluidised beds, especially when compared with the extensive studies conducted on spherical particles.
Particle-resolved direct numerical simulation provides a precise and effective approach for characterising the fluidisation of particles (Duan et al. Reference Duan, Yu, Chen and Zhou2022), with the immersed boundary method (IBM) being one of the most established and widely used techniques (Xiao et al. Reference Xiao, Zhang, Luo, Mao and Fan2022). Initially introduced by Peskin et al. (Reference Peskin2002) to investigate biological flows in the human heart, IBM has since become a standard tool for simulating flows around complex geometric shapes. This method is particularly valuable in studies of particle-laden flows, where the interaction between suspended particles and fluid dynamics is intricately coupled (Hori, Rosti & Takagi Reference Hori, Rosti and Takagi2022). Uhlann et al. (Reference Uhlmann2005) advanced this approach by developing the direct-enforced IBM, which applies no-slip boundary conditions directly onto particle surfaces to accurately model particle flows (Vanna, Picano & Benini Reference Vanna, Picano and Benini2020; Ge & Elfring Reference Ge and Elfring2025).
In direct numerical simulations, non-spherical particles (Ruan et al. Reference Ruan, Gorman, Li and Ni2022) are much less commonly studied compared with spherical particles. Shardt & Derksen (Reference Shardt and Derksen2012) developed a method for directly simulating suspensions of non-spherical rigid particles with high solid volume fractions using IBM, which successfully simulated the shape of red blood cells with solid volume fractions up to 45 %, providing a framework for studying other complex systems and phenomena. Tavassoli, Peters & Kuipers (Reference Tavassoli, Peters and Kuipers2015) utilised IBM to investigate fluid–particle heat transfer in non-spherical particles arranged in fixed random arrays, aiming to assess whether the heat transfer characteristics of spherical particles could be applied to non-spherical particles. Their findings indicated that the diameter of a spherocylinder serves as an effective parameter for characterising fluid–particle heat transfer. Jain et al. (Reference Jain, Tschisgale and Fröhlich2021) conducted direct numerical simulations of sediment transport, revealing that prolate spheroids frequently formed spanwise clusters with the highest bed friction, while triaxial ellipsoids, resembling natural sand grains, showed intermediate yet distinct results.
Direct numerical simulations of particles in fluidised beds have seen a steady increase in recent years. For instance, Reddy et al. (Reference Reddy, Sathe, Joshi, Nandakumar and Evans2013) conducted direct numerical simulations in a solid–liquid fluidised bed with a Reynolds number of 200, indicating that the average settling velocity decreased as the number of spherical particles increased. Feng & Musong (Reference Feng and Musong2014) developed a direct numerical simulation method incorporating IBM to study heat transfer in spherical particle flows and discovered that, as fluid ascends through the bed, the fluid temperature increases, leading to a reduction in the particle–fluid temperature gradient and a subsequent decrease in the particle transfer rate. Luo et al. (Reference Luo, Tan, Wang and Fan2016) applied a direct numerical simulation approach using IBM in conjunction with a soft sphere model to simulate a bubbling fluidised bed and revealed that the current Computational Fluid Dynamics-Discrete Element Method (CFD-DEM) drag correlation underestimates the particle drag force in fluidised beds when compared with the existing CFD-DEM drag model. Tang et al. (Reference Tang, Lau, Deen, Peters and Kuipers2016) performed direct numerical simulations on a two-dimensional fluidised bed. They were among the first to compare Direct Numerical Simulation (DNS) results with actual fluidising experimental data, showing that the time-averaged solid motion and pressure fluctuations were in good agreement with the experimental observations. Lan et al. (Reference Lan, Xu, Lu, Liu and Xu2024) developed a method integrating direct numerical simulation into a coarse-grained CFD-DEM solver to model the direct reduction of two iron ores with hydrogen in a fluidised bed, utilising the IBM to facilitate the use of simple meshes in reactors with complex geometries.
Unlike spherical particle fluidisation, non-spherical particle fluidisation presents a much greater complexity and poses significant challenges in achieving a comprehensive understanding. The non-spherical shapes and orientation-dependent interactions of non-spherical particles introduce additional hydrodynamic complexities that are difficult to capture using conventional modelling approaches (Sanjeevi & Padding Reference Sanjeevi and Padding2017, Reference Sanjeevi and Padding2020). Particle-resolved direct numerical simulation (PR-DNS) offers a powerful and accurate means of investigating these intricate phenomena. However, there remains a notable shortage of studies coupling the IBM with DEM-CFD approaches specifically for non-spherical particles, particularly in dense gas–solid flows. Previous efforts have primarily focused on ellipsoidal particles. For instance, Jain et al. (Reference Jain, Tschisgale and Fröhlich2019, Reference Jain, Tschisgale and Fröhlich2021) and Lan et al. (Reference Lan, Zhao, Xu, Zhao, Zhai and Wang2022) explored CFD-DEM-IBM coupling strategies for ellipsoidal particles, addressing aspects of sediment transport and fluidised bed behaviour. Nevertheless, these studies often targeted specific particle shapes, a framework for PR-DNS of particles of different shapes is still lacking in the literature.
In contrast, the present study aims to develop a generalised PR-DNS framework capable of simulating a wide variety of non-spherical particle shapes, including spheres, superellipsoids and cylinders, by coupling the IBM with a superquadric-based DEM-CFD method implemented in the open-source code MFiX (Multi-phase Flow with Interphase eXchanges). The structure of this work is organised as follows: § 2 presents the mathematical framework for coupling IBM with superDEM. Section 3 examines the flow behaviour around a superellipsoid and a cylinder, including analyses of drag, lift and torque coefficients across different Reynolds numbers. The PR-DNS results for non-spherical particle fluidisation are validated through comparisons with fluidisation experiments. Finally, the drag correction coefficients derived from the PR-DNS results are compared against traditional empirical drag models. This work provides valuable simulation insights to support future advances in the fundamental understanding of non-spherical particle fluidisation and the new drag correlations specifically for non-spherical particles in unresolved coarse-grid simulations.
2. Mathematical model
2.1. Immersed boundary method
The immersed virtual body force acting on a particle can be expressed as follows (Guo, Wu & Thornton Reference Guo, Wu and Thornton2013; Tsuji et al. Reference Tsuji, Higashida, Okuyama and Tanaka2014):
\begin{equation} \boldsymbol{F}_{\kern-1.5pt \textit{fp}}=\sum _{{c}\in R_{{p}}}\boldsymbol{f}_{\!{c}}\, ,\end{equation}
where
$R_{{p}}$
is the region where fluid cells are partially or completely covered by particles. Given that the forcing term
$\boldsymbol{f}_{\!{c}}$
equals zero if the fluid cell is fully covered by particles, the integral can be computed by iterating through all boundary cells, where the cells are partially covered by particles
\begin{equation} \boldsymbol{F}_{\kern-1.5pt \textit{fp}}=\sum _{{c}\in BC}\boldsymbol{f}_{\!{c}} .\end{equation}
The forcing term
$\boldsymbol{f}_{\!{c}}$
acting on the fluid cell is calculated as
\begin{align}& \boldsymbol{f}_{\!{c}}=\rho _{\!{f}}\frac{\big(\boldsymbol{u}_{\!{f}}-\boldsymbol{U}_{\boldsymbol{p}}^{\ast }\big)}{{\unicode[Arial]{x0394}} t}V_{{c}} , \end{align}
\begin{align}&\,\,\,\, \boldsymbol{U}_{\boldsymbol{p}}^{\ast }=\varepsilon _{{p}}\boldsymbol{U}_{\!{p}}+\varepsilon _{\!{f}}\boldsymbol{u}_{\!{f}} , \end{align}
\begin{align}&\,\, \boldsymbol{U}_{\!{p}}=\boldsymbol{u}_{\!{p}}+\boldsymbol{\omega }_{{p}}\times \boldsymbol{L}_{\textit{cp}} , \end{align}
where
$\varepsilon _{{p}}$
denotes the volume fraction of the cell that is occupied by the particle,
$V_{{c}}$
represents fluid cell volume,
$\boldsymbol{u}_{\!{p}} , \boldsymbol{\omega }_{p}$
correspond to the translational and rotational velocities of the particle, respectively, and
$\boldsymbol{L}_{\textit{cp}}$
is the vector from the fluid cell centre to the particle centre. The fluid velocity in cells that are completely covered by particles matches the particle velocity.
The torque induced by the fluid that acts on a particle can be calculated as
\begin{equation} \boldsymbol{\tau }_{\textit{fp}}=\sum _{{c}\in BC}\boldsymbol{f}_{\!{c}}\times \boldsymbol{L}_{\textit{cp}} .\end{equation}
2.2. Particle phase superDEM
The particle shape was represented using the superquadric particle method, which was developed and validated in our earlier research. The equation defining the superquadric particle surface is
\begin{equation} f\left(x,y,z\right)=\left[\left(\frac{x}{a}\right)^{\frac{2}{\varepsilon _{2}}}+\left(\frac{y}{b}\right)^{\frac{2}{\varepsilon _{2}}}\right]^{\frac{\varepsilon _{2}}{\varepsilon _{1}}}+\left(\frac{z}{c}\right)^{\frac{2}{\varepsilon 1}}-1=0,\end{equation}
where
$a$
,
$b$
,
$ c$
,
$\varepsilon _{1}$
,
$ \varepsilon _{2}$
are the parameters that determine the size and shape of the particle. By adjusting these five parameters, it is possible to represent almost 80 % of all possible shapes (Williams & Pentland Reference Williams and Pentland1992; Zhong et al. Reference Zhong, Yu, Liu, Tong and Zhang2016).
The motion of the superquadric particle is described using Newton’s laws within a Lagrangian framework. The following equations describe the particle’s position, translational and rotational velocities and orientation quaternions:
\begin{align}&\qquad\qquad\qquad \frac{{\rm d}\boldsymbol{X}_{{i}}}{{\rm d}t}=\boldsymbol{u}_{{p,i}} ,\end{align}
\begin{align}& m_{{i}}\frac{{\rm d}\boldsymbol{u}_{\textit{pi}}}{{\rm d}t}=m_{{i}}\boldsymbol{g}-\rho _{\!{f}}V_{{p}}\boldsymbol{g}+\boldsymbol{F}_{\kern-1.5pt {c,i}}+\boldsymbol{F}_{\kern-1.5pt{\textit{fp},i}} , \end{align}
\begin{align}&\quad \boldsymbol{I}_{{i}}\frac{{\rm d}\boldsymbol{w}_{{i}}}{{\rm d}t}+\boldsymbol{w}_{{i}}\times \boldsymbol{I}_{{i}}\boldsymbol{w}_{{i}}=\boldsymbol{\tau }_{{i,pp}}+\boldsymbol{\tau }_{{i,\textit{fp}}} , \end{align}
\begin{align}&\qquad\qquad\quad \frac{{\rm d}\boldsymbol{q}_{{i}}}{{\rm d}t}=\frac{1}{2}\boldsymbol{q}_{{i}}\times \boldsymbol{\omega }_{{i}} , \end{align}
where
$\boldsymbol{X}_{{i}}$
is the position vector of the particle i,
$m_{{i}}$
and
$\boldsymbol{I}_{{i}}$
are the particle mass and inertia tensor, respectively,
$\boldsymbol{u}_{{p,i}}$
and
$\boldsymbol{w}_{{i}}$
are the translational and rotational velocities, respectively,
$\boldsymbol{F}_{\kern-1.5pt {c,i}}$
denotes the contact force acting on the particle
$i$
due to contact with neighbouring particles and the wall,
$\boldsymbol{F}_{\kern-1.5pt{\textit{fp},i}}$
denotes the immersed virtual body force acting on the particle
$i$
,
$\boldsymbol{\tau }_{{i,pp}}$
and
$\boldsymbol{\tau }_{{i,\textit{fp}}}$
are the torques acting on particle
$i$
due to particle–particle and particle–fluid interactions in the body-fixed coordinate system, respectively, and
$q_{0}$
,
$q_{1}$
,
$q_{2}$
and
$q_{3}$
represent the quaternions. The conversion between the global coordinate system and the local body-fixed coordinate system is achieved through the use of a rotation matrix
\begin{equation}{Q}=\left[\!\!\begin{array}{c@{\quad}l@{\quad}l} 1-2\left(q_{2}^{2}+q_{3}^{2}\right) & 2\left(q_{1}q_{2}-q_{0}q_{3}\right) & 2\left(q_{0}q_{2}+q_{1}q_{3}\right)\\[5pt] 2\left(q_{0}q_{3}+q_{1}q_{2}\right) & 1-2\left(q_{1}^{2}+q_{3}^{2}\right)& 2\left(-q_{0}q_{1}+q_{2}q_{3}\right)\\[5pt] 2\left(-q_{0}q_{2}+q_{1}q_{3}\right)& 2\left(q_{0}q_{1}+q_{2}q_{3}\right)& 1-2\left(q_{1}^{2}+q_{2}^{2}\right) \end{array}\!\!\right]\!.\end{equation}
A two-phase contact detection algorithm was employed to identify particle overlap, utilising a broad phase with bounding spheres and oriented bounding boxes, followed by a narrow phase with superquadric contact detection.
The contact force can be resolved into its normal and tangential components, expressed as
\begin{align}&\qquad\qquad \boldsymbol{F}_{{n},{\textit{ij}}}=-k_{{n},{\textit{ij}}}\boldsymbol{\delta }_{{n},{\textit{ij}}}-\eta _{{n}}\boldsymbol{U}_{{n},{\textit{ij}}} , \end{align}
\begin{align}& \boldsymbol{F}_{{t},{\textit{ij}}}=\left\{\!\!\begin{array}{l@{\quad}l} -k_{{t},{\textit{ij}}}\boldsymbol{\delta }_{{t},{\textit{ij}}}-\eta _{{t}}\boldsymbol{U}_{{t},{\textit{ij}}}, & \left| \boldsymbol{F}_{{t},{\textit{ij}}}\right| \leq \mu _{{p}}\left| \boldsymbol{F}_{{n},{\textit{ij}}}\right|\! ,\\[5pt] -\mu _{{p}}\left| \boldsymbol{F}_{{n},{\textit{ij}}}\right| \boldsymbol{t}_{ij}, & \left| \boldsymbol{F}_{{t},{\textit{ij}}}\right| \gt \mu _{{p}}\left| \boldsymbol{F}_{{n},{\textit{ij}}}\right|\!, \end{array}\right. \end{align}
where
$\boldsymbol{\delta }_{{t},{\textit{ij}}}$
and
$\eta _{{t}}$
represent the tangential displacement and damping coefficient, respectively,
$\mu _{{p}}$
is the friction coefficient and
$k_{{n},{\textit{ij}}}$
and
$k_{{t},{\textit{ij}}}$
are the normal and tangential spring stiffnesses, respectively. The spring stiffness and damping coefficients were computed using the Hertzian model.
The torque resulting from tangential and normal forces is calculated using the following equation:
\begin{equation} \boldsymbol{\tau }_{i,{pp}}=\boldsymbol{L}_{i}\times \big(\boldsymbol{F}_{{n},{\textit{ij}}}+\boldsymbol{F}_{{t},{\textit{ij}}}\big) ,\end{equation}
where
$\boldsymbol{L}_{i}$
represents the distance vector from the particle’s centre to the contact point.
2.3. Gas phase
The gas phase is described using the Eulerian framework, with the Navier–Stokes equations governing its behaviour. The mass and momentum balances for the gas phase are expressed as
\begin{align}&\qquad\qquad\qquad\qquad\quad \boldsymbol{\nabla }\boldsymbol{\cdot } (\boldsymbol{u}_{\!{f}})=0 , \end{align}
\begin{align}& \frac{\partial (\rho _{\!{f}}\boldsymbol{u}_{\!{f}})}{\partial t}+\boldsymbol{\nabla }\boldsymbol{\cdot } (\rho _{\!{f}}\boldsymbol{u}_{\!{f}}\boldsymbol{u}_{\!{f}})=-\boldsymbol{\nabla }P+\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\tau }_{\!{f}}+\rho _{\!{f}}\boldsymbol{g}-\boldsymbol{I}_{\textit{fp}} , \end{align}
where
$\rho _{\!{f}}$
and
$\boldsymbol{u}_{\!{f}}$
are the gas density, and the gas phase velocity, respectively, and
$P , \boldsymbol{\tau }_{\!{f}} , \boldsymbol{g} , \boldsymbol{I}_{\textit{fp}}$
are the pressure, the gas phase stress tensor, the acceleration of gravity and fluid–particle interaction force per unit volume, respectively.
2.4. SuperDEM-CFD-IBM framework
Figure 1 illustrates the coupled superDEM-CFD framework integrated with the IBM model, developed and implemented in the MFiX platform. The superDEM module computes the particle-level dynamics, including motion, contact detection via hierarchical methods (superquadric, bounding sphere and oriented bounding box (OBB) tests) and force calculations. The IBM facilitates bidirectional coupling, exchanging particle-to-fluid and fluid-to-particle forces, while the Computational Fluid Dynamics (CFD) module resolves the fluid phase by generating the mesh, solving governing equations, and updating flow fields. This implementation enables accurate simulations of particle–fluid interactions in complex multiphase systems.
Flow chart of the coupling superDEM-CFD with the IBM.
2.5. Particle suspension simulation set-up
Figure 2 presents the set-up of a three-dimensional simulation for PR-DNS analysis. Figure 2(a) illustrates the schematic diagram of the simulation domain and dimensions with an inlet flow velocity
$U_{\infty }$
, detailing that the particle centre is 20d from the inlet, 40d from the pressure outlet boundary and 20d from the free slip boundary, where d is the volume-equivalent diameter of the particle. Figure 2(b) displays the computational grid, where the grids are refined around the particle. Figures 2(c) and 2(d) illustrate the volume fraction of fluid cells and fluid velocity in the x-o-y plane. Figure 2(e) depicts the forces acting on the particle, including drag force (
$F_{\!D}$
), lift force (
$F_{\!L}$
), torque (T) and angle of attack (
$\phi$
). The drag coefficient
$C_{\!D}$
, the lift coefficient
$C_{\!L}$
and the torque coefficient
$C_{\!T}$
are calculated as follows:
\begin{align}& C_{\!{D}}=\frac{2F_{\!{D}}}{\rho _{\!{f}}U_{\infty }^{2}\pi d_{{p}}^{2}/4} , \end{align}
\begin{align}& C_{{L}}=\frac{2F_{{L}}}{\rho _{\!{f}}U_{\infty }^{2}\pi d_{{p}}^{2}/4} , \end{align}
\begin{align}& C_{{T}}=\frac{2T}{\rho _{\!{f}}U_{\infty }^{2}\pi d_{{p}}^{3}/8} , \end{align}
Set-up of simulation. (a) Schematic diagram of the simulation domain, (b) CFD grid, (c) volume fraction of fluid cell occupied by particle, (d) contour plot of fluid velocity, (e) force diagram.
where
$\rho _{\!{f}}$
represents the density of the fluid,
$U_{\infty }$
represents the inlet velocity and d
p
is the volume-equivalent sphere diameter of the particle.
The computational domain is set to 6 m × 4 m × 4 m in the x, y and z directions. The CFD cells in these directions are 180 × 120 × 120, corresponding to 30 cells ranging from 0 to 1.9 m in all directions and then 60 cells ranging from 1.9 to 2.1 m in all directions; there are 30 cells ranging from 1.9 to 4.0 m in the y and z directions, and 90 cells ranging from 2.1 to 6.0 m in the x direction. The size of the CFD cells is uniform cubic in the three-dimensional region (x = 1.9–2.1 m, y = 1.9–2.1 m, z = 1.9–2.1 m). The strategy of adjusting the time step automatically within user-specified limits is implemented in the MFiX code to reduce the computational time. Table 1 summarises the particle properties and parameters utilised in the single-particle simulations.
Properties of particles and parameters used in the single-particle simulations.
Properties of particles and parameters used in the fluidised bed simulations.
All PR-DNS simulation cases in this work.
2.6. Fluidisation simulation set-up
To validate the PR-DNS for non-spherical particle fluidisation, experimental data from Vollmari, Jasevičius & Kruggel-Emden (Reference Vollmari, Jasevičius and Kruggel-Emden2016) were utilised for comparison, where the fluidisation behaviour of spherical and non-spherical particles was investigated in a rectangular fluidised bed, with a height of 0.41 m and a cross-sectional area of 0.11 × 0.11 m2. Flow straighteners were installed below the distributor to ensure a fully developed inlet gas flow. Pressure drops of bed material particles were measured, which connected two pressure tapping points, one installed at the bed height of 0.31 m and the other 0.02 m below the porous distributor. The superficial gas velocity of pressure drop experiments ranged from 0 to 2.4 m s−1, covering both the fixed bed and fluidised bed regimes. A high-speed digital camera placed at the front of the fluidised bed was used to record videos. The images extracted from the recorded videos were processed using the threshold method to analyse particle height and orientation distributions. In this work, the experimental results from Vollmari et al. (Reference Vollmari, Jasevičius and Kruggel-Emden2016), including the bed pressure drop at different gas velocities and the height and orientation distributions of cylindrical particles, were used to validate the PR-DNS of fluidisation. Table 2 presents the properties of particles and simulation parameters used in the PR-DNS of fluidisation, which are the same as our previous particle-unresolved superDEM-CFD fluidisation simulations (Gao et al. Reference Gao, Yu, Lu, Li and Rogers2021a ). The dimensions for the sphere (2a, 2b, 2c) are all 7.2 mm, while for the cylinder, they are 6.1, 6.1 and 6.2 mm, respectively. The combined mass of the spheres or cylinders in the fluidised bed is 0.535 kg. The particle volume-equivalent sphere diameter/fluid grid resolution is set to 9.4 for the fluidisation of spheres simulation and 8.0 for the fluidisation of cylinder simulations. In this work, all the PR-DNS results were conducted using the HPC stations in our laboratory equipped with AMD Epyc 7763 processors. For the single-particle simulations, 45 threads were used, and a physical time of 30 000 s (for example, at Reynolds number – Re = 0.1) required approximately 80 h. For fluidisation simulations, 192 threads were used, and a physical time of 10 s required approximately 1000 h.
All PR-DNS simulation cases of different particle shapes and conditions in this study are summarised in table 3.
3. Results and discussion
3.1. Flow passing through a sphere
Figure 3 presents the velocity field contours around a spherical object for two different Reynolds numbers. The Reynolds number is calculated using the volume-equivalent diameter of the particle. In figure 3(a), with Re = 1, the flow is characterised by a symmetric and smooth pattern around the sphere, indicating a laminar flow regime where viscous forces dominate. The velocity gradient is gradual, and there is no significant wake formation behind the sphere. In contrast, figure 3(b) for Re = 100 shows a more complex flow pattern. The flow detaches from the surface of the sphere, forming a distinct wake region with lower velocities behind the sphere. The velocity field exhibits slight asymmetry, and the presence of recirculating zones is evident, indicative of the transition toward weak turbulent flow. The colour bar on the right of each plot quantifies the velocity magnitudes, where red represents higher velocities and blue indicates lower velocities. This figure shows that the flow of gas around a spherical object is a stable laminar flow when Re = 1. When Re = 100, the flow changes from a stable laminar flow to an unstable separated flow.
Velocity field contours of U x around a sphere at (a) Re = 1 and (b) Re = 100.
Figure 4(a) presents a comparative analysis of the drag coefficient (
$C_{\!{D}}$
) as a function of the Reynolds number (Re), featuring three different sets of results: Clift (Clift, Grace & Weber Reference Clift, Grace and Weber2013), Stokes (Happel & Brenner Reference Happel and Brenner2012) and the present PR-DNS. Clift’s drag model shows how the drag coefficient varies across a wide range of Reynolds numbers. Stokes’ drag model is applicable in the low Reynolds number regime and shows a rapid decrease in drag coefficient. The current PR-DNS results span the entire range of Reynolds numbers from 0.1 to 1000. It can be observed that the RP-DNS aligns well with Stokes’ model at low Reynolds numbers (Re <<1), and also matches Clift’s results across the broader range. The accuracy of the present PR-DNS was effectively demonstrated by comparing it with existing models. Figure 4(b) illustrates the present PR-DNS of the normalised drag coefficient as a function of the ratio (particle diameter Dp
/fluid cell resolution) for Reynolds numbers Re = 1 and Re = 100. In order to study mesh convergence, the ratio varies between 4 and 30. The particle Reynolds number is calculated as
${\textit{Re}}=\rho _{\!{f}}d_{{p}}U_{\infty }/\mu _{\!{f}}$
. At a Reynolds number of 1, the impact of mesh independence was minimal; however, at a Reynolds number of 100, the effect of mesh independence became more prominent. As Dp
/cell increases, the normalised drag coefficient converges to 1 for both Re = 1 and Re = 100. Beyond Dp
/cell of 20, the values are almost equal to 1, indicating grid independence. The errors are below 10 % when Dp
/cell > 8. These findings align closely with the results reported by Sanjeevi, Kuipers & Padding (Reference Sanjeevi, Kuipers and Padding2018).
(a) A comparison between the present drag coefficient for the flow past a sphere and Clift and Stokes values. (b) Effect of grid resolution on the normalised PR-DNS drag coefficient.
3.2. Flow passing through a superellipsoid and cylinder
Figure 5 presents the velocity contours around a superellipsoid and a cylinder at Re = 90 and
$\phi$
= 45°. The colour bar on the right side indicates the velocity magnitude (m s−1), with red representing higher velocities and blue representing lower velocities. The velocity contours highlight the flow separation and wake regions behind the particle. For the superellipsoid in figure 5(a), the flow smoothly wraps around the shape with a gradual wake formation, whereas for the cylinder in figure 5(b), the flow exhibits a more pronounced separation and a longer wake region. For the negative value of the velocity field U
x
, since at Re = 90, the flow regime is not Stokes flow, and secondary flow exists behind the particle (Sanjeevi et al. Reference Sanjeevi, Kuipers and Padding2018), and this phenomenon was more evident for the superellipsoid and cylinder than a sphere. This visualisation effectively demonstrates the differences in flow characteristics and wake formation between a superellipsoid and a cylinder.
Velocity contours of U x around a superellipsoid (a) and a cylinder (b) at Re = 90.
Studies on the forces acting on non-spherical particles using PR-DNS (Goossens Reference Goossens2019; Mola, Fawell & Small Reference Mola, Fawell and Small2020), including the hydrodynamic coefficients (
$C_{\!{D}} , C_{{L}} , C_{{T}}$
), remain limited (Shao, Zhong & Yu Reference Shao, Zhong and Yu2016). Figure 6 illustrates the drag coefficient (
$C_{\!{D}}$
), lift coefficient (
$C_{{L}}$
) and torque coefficient (
$C_{{T}}$
) of a superellipsoid at two distinct Reynolds numbers: Re = 0.3 and Re = 90, where the present PR-DNS results are compared with the simulation results reported by Hölzer & Sommerfeld (Reference Hölzer and Sommerfeld2008), Ouchene et al. (Reference Ouchene, Khalij, Arcen and Tanière2016) and Kempe & Fröhlich (Reference Kempe and Fröhlich2012). The drag coefficient is shown as a function of the angle of attack, highlighting the relationship between particle orientation and hydrodynamic resistance under varying flow conditions. The drag coefficient at Re = 0.3 increases with the angle for all simulations in figure 6(a), validating the accuracy of the present PR-DNS. Figure 6(b) shows the drag coefficient at Re = 90, which also increases with the angle of attack for all curves, demonstrating the accuracy of the present PR-DNS compared with existing methods. The results from the present PR-DNS align well with those of previous studies, validating the simulation model and indicating a consistent trend across different flow regimes. In figure 6(c), the lift coefficient decreases as the angle increases to around 45 degrees, reaches a minimum and then increases towards 90 degrees for all curves. Figure 6(d) shows a similar trend, decreasing to a minimum of around 45 degrees and then increasing towards 90 degrees. The present PR-DNS results align closely with previous studies, confirming the model’s reliability. In figure 6(e), the torque coefficient decreases as the angle increases to around 45 degrees, reaches a minimum and then increases towards 90 degrees for results, illustrating a consistent trend and validating the simulation model. Figure 6(f) presents the torque coefficient at Re = 90. All the curves display a similar tendency, with the torque coefficient diminishing to a minimum of approximately 45 degrees and then increasing again towards 90 degrees. Ouchene et al. (Reference Ouchene, Khalij, Tanière and Arcen2015) also presented a study on hydrodynamic coefficients for ellipsoidal particles at different Reynolds numbers (Re = 10 and Re = 290), which revealed that the drag coefficient rises along with the increase of the angle at different Reynolds numbers, and the extreme values of both the lift coefficient and the torque coefficient are at 45 degrees. Sanjeevi et al. (Reference Sanjeevi, Dietiker and Padding2022) carried out simulations on prolate spheroids at Reynolds numbers of 10 and 100 and reported the drag, lift and torque coefficients. Their findings also indicated that the drag coefficient increases with the angle, while the extreme values of the lift and torque coefficients occur symmetrically at 45 degrees. The results from the present PR-DNS are in strong agreement with those from previous studies, validating the simulation model and suggesting a consistent trend across different flow regimes, confirming the model’s reliability.
Values of
$C_{\!{D}}$
of a superellipsoid at Re = 0.3 (a) and Re = 90 (b),
$C_{{L}}$
of superellipsoid at Re = 0.3 (c) and Re = 90 (d),
$C_{{T}}$
of superellipsoid at Re = 0.3 (e) and Re = 90 (f).
Figure 7 shows the drag coefficient (
$C_{\!{D}}$
), lift coefficient (
$C_{{L}}$
) and torque coefficient (
$C_{{T}}$
) of cylinders as functions of the angle of attack at Re = 0.3 and Re = 90 and a comparison with the simulation results reported by Hölzer & Sommerfeld (Reference Hölzer and Sommerfeld2008). In figure 7(a), the drag coefficient increases with the angle for both curves, with the present PR-DNS showing a slightly lower drag coefficient compared with Holzer and Sommerfeld’s results. Similar to figure 7(a), the drag coefficient shown in figure 7(b) increases with the angle for both curves, but the drag coefficient by the present PR-DNS is slightly larger than that of Holzer and Sommerfeld’s results. In figure 11(c), both curves show a decrease in the lift coefficient to a minimum value of around 45 degrees, followed by an increase towards 90 degrees, with the present PR-DNS results being slightly higher than those of Holzer and Sommerfeld; figure 7(e) also shows the same trend about the torque coefficient at Re = 0.3. However, at higher Reynolds numbers Re = 90 in figure 7(d), the minimum lift coefficient point of PR-DNS is at 20 degrees, while the minimum torque coefficient point is at 70 degrees in figure 7(f), which is different from the results of Holzer and Sommerfeld.
Values of
$C_{\!D}$
of the cylinder at Re = 0.3 (a) and Re = 90 (b),
$C_{\!L}$
of the cylinder at Re = 0.3 (c) and Re = 90 (d),
$C_{\!T}$
of the cylinder at Re = 0.3 (e) and Re = 90 (f).
Non-spherical particles inclined with respect to a homogeneous flow experience torque, and such orientation will cause complexity when analysing the lift and torque coefficients (Sanjeevi & Padding Reference Sanjeevi and Padding2017, Reference Sanjeevi and Padding2020). There is an observed asymmetry in the lift and torque coefficients of the cylinder as the Reynolds number increases from 0.3 to 90. As the Reynolds number rises toward 90, the flow transitions into a regime separated by the secondary flow, as shown in figure 5(b). For the cylinder, due to the sharp edges compared with a superellipsoid, the flow pattern and evolution of the possible instabilities will be further affected due to the pressure distribution and vortex-shedding dynamics, resulting in the asymmetry in the lift and torque coefficients with respect to the angle of attack.
Besides, as shown in figures 7(d) and 7(f), there is an observed asymmetry in the lift and torque coefficients of the cylinder as the Reynolds number increases from 0.3 to 90. The asymmetry of the lift and torque coefficients for the spherocylinder (AR = 4) and prolate spheroid (AR > 4) at Re = 100 was also reported by Sanjeevi et al. (Reference Sanjeevi, Kuipers and Padding2018), and as the Reynolds number increases, the asymmetry becomes more significant. The asymmetry in lift and torque coefficients of the cylinder can be further explained by investigating the wake structure and vortex dynamics through DNS simulation in the future.
In line with present PR-DNS findings, Lote, Yang & Mao (Reference Lote, Yang and Mao2024) investigated the drag, lift and torque coefficients for aspect ratio 0.6 cylinders. Their results showed that, as the Reynolds number increases, the extreme values of the lift and torque coefficients shift from 45 degrees to 25 degrees and 45 degrees to 75 degrees, respectively. Wang et al. (Reference Wang, Ma, Jiang, Fang, Yin, Tan, Zhang and Chen2024c ) also reported consistent findings, showing that, for longer cylinders at low Reynolds numbers (Re = 10), the minimum lift coefficient occurs at 45 degrees, and as the Reynolds number increases, the angles corresponding to the minimum lift coefficients shift from 45 degrees to 30 degrees. Future studies employing alternative methods are encouraged to elucidate the causes of these observed differences further.
3.3. Flow passing through a superellipsoid and cylinder with different aspect ratios
To further validate the accuracy of the present PR-DNS, the simulated drag coefficients of superellipsoids with different aspect ratios (AR) were compared with analytical solutions and previous results from the literature. The effect of AR on the drag coefficient of superellipsoids and cylinders at different angles of attack was also investigated. All the single particles have the same volume and density. The settings of CFD cells and time steps are the same as mentioned in § 2.5. Table 4 summarises the parameters of superellipsoids and cylinders with different ARs used in the PR-DNS.
Parameters of superellipsoids and cylinders with different ARs used in simulations.
The theoretical solutions of drag coefficients for superellipsoids varying the AR at different angles of attack were derived by Happel & Brenner (Reference Happel and Brenner2012). The detailed formulas are as follows:
\begin{equation} C_{{D,}\phi ={0}^{\circ}}=\frac{64}{{\textit{Re}}\boldsymbol{\cdot }\textit{AR}^{1/3}}\frac{1}{\left[-\frac{2\textit{AR}}{\textit{AR}^{2}-1}+\frac{2\textit{AR}^{2}-1}{\left(\textit{AR}^{2}-1\right)^{3/2}}\boldsymbol{\cdot }\ln \left(\frac{\textit{AR}+\sqrt{\textit{AR}^{2}-1}}{\textit{AR}-\sqrt{\textit{AR}^{2}-1}}\right)\right]} ,\end{equation}
\begin{equation} C_{{D,}\phi ={90}^{\circ}}=\frac{64}{{\textit{Re}}\boldsymbol{\cdot }\textit{AR}^{1/3}}\frac{1}{\left[\frac{\textit{AR}}{\textit{AR}^{2}-1}+\frac{2\textit{AR}^{2}-3}{\left(\textit{AR}^{2}-1\right)^{3/2}}\boldsymbol{\cdot }\ln \big(\textit{AR}+\sqrt{\textit{AR}^{2}-1}\big)\right]} ,\end{equation}
where AR is the aspect ratio of superellipsoids, Re is the Reynolds number calculated using the volume-equivalent diameter of the particle and
$\phi$
is the angle of attack illustrated in figure 2.
Due to the linearity of solutions in the Stokes flow at Re = 0, the lift coefficient
$C_{{L}}$
of superellipsoids can calculated as follows when the AR is greater than unity using theoretical solutions from Happel & Brenner (Reference Happel and Brenner2012):
\begin{equation} C_{{L,\textit{AR}}}=\left(C_{{D,}\phi ={90}^{\circ}}-C_{{D,}\phi ={0}^{\circ}}\right)\sin \phi \cos \phi .\end{equation}
Figure 8 presents the velocity field around superellipsoids and cylinders of different ARs at Re = 0.1. Generally, the velocity fields of superellipsoids and cylinders with different ARs are similar and symmetric in the flow direction. Since the flow is close to Stokes flow at this low Reynolds number, streamlines are close to the particle’s surface without secondary flow, and the minimum circumscribed diameters of superellipsoids and cylinders are close, resulting in a similar velocity distribution.
The velocity of U x around non-spherical particles at Re = 0.1; (a) superellipsoid with AR = 2, (b) superellipsoid with AR = 4, (c) cylinder with AR = 2 and (d) cylinder with AR = 4.
Figure 9 compares the predicted
$C_{\!{D}}$
for superellipsoids of different ARs at different angles of attack with theoretical solutions and the literature results. The predicted
$C_{\!{D}}$
agrees well with theoretical solutions of different ARs and angles of attack, which verifies the accuracy of the current PR-DNS. The predicted
$C_{\!{D}}$
by the current PR-DNS is almost the same as the simulations by Sanjeevi et al. (Reference Sanjeevi, Dietiker and Padding2022) using the multi-relaxation time lattice Boltzmann method. The differences between the predicted
$C_{\!{D}}$
and theoretical solutions are mainly because of the finite inertia (Re = 0.1). As shown in figure 7(a), when
$\phi$
= 0°, as the AR increases, the lift coefficient
$C_{\!{D}}$
decreases first and then increases, and the minimum reaches when the AR is near 2, which was also revealed by Sanjeevi et al. (Reference Sanjeevi, Dietiker and Padding2022). The case
$\phi$
= 0° means the superellipsoid is more aligned with the stream, resulting in the initial drop of
$C_{\!{D}}$
, and then, due to the increasing surface area as the AR increases,
$C_{\!{D}}$
increases. When
$\phi$
= 90°, the lift coefficient
$C_{\!{D}}$
increases because of the increasing surface area as the AR increases.
Comparison of
$C_{\!{D}}$
for superellipsoid with different ARs at Re = 0.1; (a)
$\phi$
= 0° and (b)
$\phi$
= 90°.
Figure 10 compares the predicted
$C_{\!{D}}$
for superellipsoids and cylinders of different ARs at different angles of attack for Re = 0.1 quantitatively. As shown in figure 10, the drag coefficients of cylinders of different ARs are larger than those of superellipsoids. When
$\phi$
= 0°, the drag coefficients of the cylinders are 6.32 % and 6.25 % larger than superellipsoids when AR equals 2 and 4, respectively. When
$\phi$
= 90°, the drag coefficients of cylinders are 3.83 % and 3.82 % larger than superellipsoid, respectively. The cylinder and spherocylinder have almost the same shape except for the two ends of the axis when AR = 4. Their drag coefficients are expected to be similar, which can be seen from the differences in the drag coefficients of the cylinder and spherocylinder when
$\phi$
= 0° and 90°, which are 1.94 % and 2.94 %, respectively.
Values of
$C_{\!{D}}$
of superellipsoids and cylinders with different ARs at Re = 0.1; (a)
$\phi$
= 0° and (b)
$\phi$
= 90°.
Figure 11(a) compares the predicted
$C_{{L}}$
for superellipsoids of different ARs at
$\phi$
= 45° with the theoretical solution and literature results. The predicted
$C_{{L}}$
agrees well with the theoretical solution and literature simulations of different ARs, which further verifies the accuracy of the current PR-DNS. Figure 11(b) compares the predicted
$C_{{L}}$
for superellipsoids and cylinders of different ARs when
$\phi$
= 45° for Re = 0.1 quantitatively. The lift coefficients of cylinders are 11.22 % and 3.91 % smaller than superellipsoids when AR equals 2 and 4, respectively. The ignorable difference in the lift coefficient of the cylinder and spherocylinder is 2.55 % when AR = 4, mainly due to the similar particle shape.
Values of
$C_{{L}}$
of superellipsoids with different ARs at Re = 0.1 (a) compared with theoretical and other results, and (b) compared with cylinders.
3.4. Fluidisation of spheres and cylinders
Both PR-DNS of spheres and cylinders in a fluidised bed were conducted, where the simulation of spherical particle fluidisation serves as the fundamental case, while cylindrical particles are among the most representative shapes in fluidised bed applications. This study focuses on validating the PR-DNS results against experimental data (Vollmari et al. Reference Vollmari, Jasevičius and Kruggel-Emden2016). Figure 12 shows snapshots of the fluidisation behaviour and velocities of spherical and cylindrical particles at a superficial gas velocity of 1.6 m s−1. The y direction is the same as the gravity, and solid and gas velocities are components of the y direction. A total of 10 s of physical time was simulated, and the simulation results of the last 6 s were post-processed. The comparison between spherical and cylindrical particles highlights the differences in the fluidisation dynamics. It can be seen that cylindrical particles are blown higher than spherical particles at the same superficial gas velocity. This observation aligns with the experiment reported by Dai et al. (Reference Dai, Ma, Sheintuch and Gao2024), who showed that the expansion height of cylindrical particles is greater than that of spherical particles at the same fluidising inlet velocity. Since the spherical and cylindrical particles in the fluidised bed were initially packed under gravity. Before the particles fluidised, the inlet gas bypassed through the pore structures formed by the interlocked particles of the fixed bed (Mema et al. Reference Mema, Wagner, van Ommen and Padding2020; Mema & Padding Reference Mema and Padding2021), resulting in the typical channel flow phenomenon (Vollmari et al. Reference Vollmari, Jasevičius and Kruggel-Emden2016; Mahajan et al. Reference Mahajan, Padding, Nijssen, Buist and Kuipers2018; Gao et al. Reference Gao, Yu, Lu, Li and Rogers2021a ). Compared with spheres, cylinders tend to form more complex interlocking and pore structures due to the non-spherical shape, resulting in more pore structures and larger voids for the same cross-section of the bed, meaning the smaller effective cross-section of particles and the smaller minimum fluidised velocity. As a result, the expansion height of cylinders exceeds that of spheres at the same inlet velocity after complete fluidisation.
Snapshots of fluidisation behaviour of spherical and cylindrical particles at 1.6 m s−1. The gas velocity is the field of U y sliced from the central cross-section of the bed.
The reliability of the present PR-DNS method is verified by comparing the simulation results with experimental pressure drops. The particle-unresolved superDEM-CFD simulation results reported in the previous work (Gao et al. Reference Gao, Yu, Lu, Li and Rogers2021a ) are also included in the comparison, as shown in figure 13. In figure 13(a), for the sphere fluidisation simulation, the present PR-DNS is more consistent with the experimental results than the particle-unresolved superDEM-CFD simulations using the Di Felice–Ganser drag model and the Di Felice–Holzer/Sommerfeld drag model in both the fixed bed and fluidised bed regimes. In figure 13(b), for the cylinder fluidisation simulation, both the PR-DNS and the particle-unresolved simulation slightly underestimated the pressure drop in the fixed bed regime. When the superficial gas velocity exceeds 1.0 m s−1, where the cylinders are fluidised, the present PR-DNS aligns better with the experimental data than the particle-unresolved simulations because the PR-DNS allows for precise resolution of the flow around the particles. The pressure drop of spheres predicted by PR-DNS is in the best agreement with experiments, which shows that the current PR-DNS for spheres can adequately analyse the interaction between spheres and gas. However, for cylinders, the pressure drop predicted by PR-DNS is not in the best agreement with the experiments. The pressure drop gradient varies significantly in the superficial gas velocity range of 0.4–0.8 m s−1, because the channel flow effect of non-spherical particles is more significant than that of spheres, it is still challenging to effectively simulate the channel flow phenomenon of non-spherical particles in fluidised beds. Compared with spheres, the current PR-DNS of cylinders is insufficient to capture the channel flow where the velocity gradient is significantly higher than in other regions. Therefore, it is necessary to employ other PR-DNS methods or use an extremely fine fluid grid.
Comparing predicted pressure drop with experimental data using different drag models; (a) sphere, (b) cylinder.
The particle height distribution reflects the bed expansion height and serves as a key metric for further validating the simulation of particles with different shapes, thus assessing the model’s accuracy. Figure 14 shows the comparison of the spherical and cylindrical particle height distributions in the present PR-DNS with experimental data and unresolved superDEM-CFD simulation results at a gas velocity of 1.6 m s−1. For spherical particles, as shown in figure 14(a), both PR-DNS and superDEM-CFD predict similar overall trends in the height distribution. For the cylinder fluidisation simulation, as shown in figure 14(b), the PR-DNS predicted a slightly lower expansion height than the experimental data, while the unresolved simulation predicted a slightly higher expansion height than the experimental data. The PR-DNS predicted a slightly higher expansion height of spheres than that of the unresolved simulations. The experiment and simulation of Dai et al. (Reference Dai, Ma, Sheintuch and Gao2024), which latter was carried out using superDEM-CFD, disclosed a similar tendency. At the same fluidisation velocity, spherical particles exhibited a slightly higher height distribution than cylindrical particles at lower particle heights. Nevertheless, as the particle height rose, the frequency of spherical particles decreased rapidly, while the frequency of cylindrical particles decreased at a slower pace.
Comparison of the predicted particle height distribution; (a) sphere, (b) cylinder.
Figure 15 displays the predicted distribution of cylinder particle orientations alongside experimental data. The orientations are measured in degrees, with 0° indicating a horizontal or lying particle and 90° a vertical or standing particle. The orientations are categorised into 10-degree bins (9 bins total), consistent with the experimental methodology. The frequency within each bin is determined by counting the number of particles it contains. The frequency distribution reveals that both the current PR-DNS and the unresolved simulations closely follow the experimental data. As the orientation angle increases, the frequency diminishes, indicating a general trend captured by both the present PR-DNS and the superDEM-CFD simulation. Wang et al. (Reference Wang, Yang, Yang and Pang2024d ) probed the particle shape and its interactions with the surrounding fluid by adopting a CFD-DEM approach coupled with the super-quadric model for cylindrical particles, and they noticed a similar trend with PR-DNS in the orientation distribution of cylindrical particles, with the frequency decreasing from 0 degrees to 90 degrees.
Comparing predicted cylinder particle orientation distribution with experimental data.
3.5. Drag correlation coefficient
In unresolved CFD-DEM simulations, the flow around each particle is not resolved, and a drag model derived from experiments or PR-DNS is required. In PR-DNS, the flow field around each particle is resolved, allowing for the computation of drag force for each particle. The drag coefficient from PR-DNS can be related to the drag force as follows:
\begin{equation} \beta _{\textit{DNS}}=\frac{\boldsymbol{F}_{\kern-1.5pt{\textit{fp},i}}}{(\boldsymbol{u}_{\kern-1.5pt{f,cg}}-\boldsymbol{u}_{{P}})V_{{i}}/(\varepsilon _{\kern-1.5pt{f,cg}}(1-\varepsilon _{\kern-1.5pt{f,cg}}))} ,\end{equation}
where
$\beta _{\textit{DNS}}$
is the drag correlation coefficient,
$\boldsymbol{F}_{\kern-1.5pt{\textit{fp},i}}$
is the fluid force obtained from PR-DNS and
$\boldsymbol{u}_{\kern-1.5pt{f,cg}}-\boldsymbol{u}_{{P}}$
is the slip velocity between the fluid and the particles calculated in a coarse grid, where fluid grid resolution/particle diameter is set to 1.9. Also
$V_{{i}}$
denotes the particle volume and
$\varepsilon _{\kern-1.5pt{f,cg}}$
is the cell voidage calculated in a coarse grid. The PR-DNS-based drag coefficient can be calculated using (3.4), and the drag correlation coefficients of existing models can be calculated using the expression listed in table 5 with other necessary data calculated in the coarse grids.
Drag correlation coefficient expressions of BVK, Di Felice–Ganser and Di Felice–Holzer/Sommerfeld models.
Figure 16 compares the drag correlation coefficients derived from the PR-DNS and the drag correlation employing different drag models. The results indicate that Di Felice–Holzer/Sommerfeld drag correlation values align better with the drag correlation values calculated using PR-DNS results than the Di Felice–Ganser and BVK (Bokkers, van Sint Annaland & Kuipers Reference Bokkers, van Sint Annaland and Kuipers2004) drag correlations. This finding supports our previous finding that the results of particle-unresolved simulation using the Di Felice–Holzer/Sommerfeld hybrid drag model agree better with the experimental data than those obtained using the Di Felice–Ganser hybrid drag model (Gao et al. Reference Gao, Yu, Lu, Li and Rogers2021a ). Mema & Padding (Reference Mema and Padding2021) utilised the Di Felice–Holzer/Sommerfeld drag model to simulate a fluidised bed with spherocylindrical particles. Their study demonstrated that the Di Felice–Holzer/Sommerfeld drag model exhibited better consistency with experimental results, particularly in the free-board region. Brömmer et al. (Reference Brömmer, Scharnowski, Illana Mahiques, Wirtz and Scherer2024) applied the Di Felice–Holzer/Sommerfeld and Ergun drag models to simulate airflow through narrow openings into a particle aggregate consisting of spherical and dodecahedral particles in a flat-bottom hopper. Their findings indicated that the Di Felice–Holzer/Sommerfeld model provided more accurate predictions for both spherical and complex-shaped particles.
Comparison of the drag correlation coefficients derived from the PR-DNS and the unresolved simulation employing different drag models; (a) sphere, (b) cylinder.
4. Conclusions
In conclusion, a novel PR-DNS method for non-spherical particles was developed in the open-source code MFiX to simulate the suspension and fluidisation of particles. Through coupling superquadric DEM-CFD with the IBM, the model allows for the analysis of fluid dynamics coefficients (
$C_{\!{D}} , C_{{L}} , C_{{T}}$
) for superellipsoids and cylinders within a range of Reynolds numbers. The pressure drop, height distribution and orientation distribution of cylinder particles simulated in the fluidised bed were compared with experimental data, offering a novel numerical approach to explaining experimental results. Detailed conclusions are as follows:
-
(i) The coupled model shows great predictivity of the drag, lift and torque coefficients of non-spherical particles at different Reynolds numbers, with the angle of attack ranging from 0° to 90°, and superellipsoid and cylinder particles with different ARs up to 4.
-
(ii) The PR-DNS results show that the particle shape has a significant effect on the lift and torque of cylindrical particles as the Reynolds number increases from 0.3 to 90, resulting in the asymmetry of the lift and torque coefficients with the angle of attack, which provides new theoretical guidance for the direct numerical simulation of non-spherical particles.
-
(iii) The PR-DNS results of sphere fluidisation agree better with the experimental data than those of the cylinder fluidisation in the fixed-bed regime, with the superficial gas velocity ranging from 0 to 1.0 m s−1, indicating the need for further study of the non-spherical particle fluidisation experiment and PR-DNS to capture the channel flow phenomenon sufficiently.
-
(iv) By comparing the drag correlation coefficients of PR-DNS and existing drag models, it was observed that the Di Felice–Holzer/Sommerfeld model aligns more closely with the PR-DNS results, consistent with our previous work that this model matches experimental results better in the particle-unresolved simulation.
Acknowledgements
The work was supported by the Basic and Applied Basic Research Foundation of Guangdong Province (2025A1515011017), Guangdong Provincial Science and Technology Innovation Strategy Special Fund (STKJ202209056), National Natural Science Foundation of China (21FAA02728) and MATEC sprout project (KD2200068).
Declaration of interests
The authors report no conflict of interest.
Open access
