2.1. Grain shape generation and description

In concurrent studies on drag coefficients of grains, the ideal geometry with fully convex shapes is of high priority. However, realistic grain surfaces usually demonstrate fractality manifested by irregular tortuosity over a wide range of length scales (Barclay & Buckingham Reference Barclay and Buckingham2009). In this study, we adopt the inverse analysis of SH with prescribed power spectrum to generate such multi-scaled morphology features (Wei et al. Reference Wei, Wang, Pereira and Gan2021); the power spectrum is enriched with fractality and root mean square roughness. A rough grain constructed with star-shaped surficial points can be approximated using the SH function ${Y}_{n}^{m}(\theta ,\varphi )$ , via denoting its radial length in the spherical coordinate system as

(2.1) \begin{equation}r(\theta ,\varphi )=\sum\nolimits_{n=0}^{\infty }\sum\nolimits_{m=-n}^{n}{c}_{n}^{m}\,{Y}_{n}^{m}(\theta ,\varphi ),\end{equation}

where $\theta \in [0,\pi ]$ and $\varphi \in (0,2\pi ]$ are the latitudinal and longitudinal coordinates, respectively, and ${c}_{n}^{m}$ are the SH coefficients of spherical wavenumber $n$ and order $m$ . With the help of Parseval’s theorem, at each SH frequency, the power spectrum is $L_{n}=\|f_{n}\|=\sqrt{\sum_{m=-n}^{n}\|{c}_{n}^{m}\|^{2}}$ , with $\|\cdot\|$ being the L ₂ norm. The average radius length $\overline{r}$ , equal to its $c_{0}$ -determined sphere $R_{0}$ , reads as $\overline{r}\approx R_{0}=c_{0}\,Y_{0}(\theta ,\varphi )$ , where $Y_{0}(\theta ,\varphi )=1/(2\sqrt{\pi })$ . To quantify how a rough grain surface is globally different from $\overline{r}$ , the relative roughness is defined as the ratio $R_{r}=S_{q}/\overline{r}$ , i.e. the root mean square $S_{q}=\sqrt{(\sum_{n=1}^{\infty }\sum_{m=-n}^{n}\|{c}_{n}^{m}\|^{2})/(4\pi )}$ over $\overline{r}$ . With $R_{r}$ and $c_{0}$ in hand, the grain volume can be correspondingly determined as $V_{\!p}=c_{0}^{3}(1+3R_{r}^{2})/(6\sqrt{\pi })$ , which eases the definition of effective grain diameter $d_{\textit{eq}}=\sqrt[3]{6V_{\!p}/\pi }=c_{0}\sqrt[3]{1+3R_{r}^{2}}/\sqrt{\pi }$ ; the grain surface area is approximately $S_{a}=c_{0}^{2}(1+\pi R_{r}^{\pi /2}D_{\!f}^{3.874})/20$ (Wei et al. Reference Wei, Wang, Pereira and Gan2021). Meanwhile, the multi-scale features – denoting how a rough grain surface is globally different from $\overline{r}$ – can be characterised by fractal dimension $D_{\!f}$ ; it is originated from the logarithmic linear relations between $L_{n}$ and $n$ , i.e. $L_{n}\propto n^{\beta }$ , where $\beta =-2H$ is the slope of the regression plot of $\log(L_{n})$ versus $\log(n)$ , and $H$ is the Hurst coefficient used to calculate $D_{\!f}=3-H$ (Russ Reference Russ2013). With increasing $n$ , more and more fine surficial details could be depicted. Considering the extremely low ${\textit{Re}}=10^{-5}$ , $n_{max}=15$ is fixed considering the related cut-off or minimum wavelength $\lambda (n)={\pi }\overline{r}/n={\pi }\overline{r}/15$ according to the famous Jean’s formula.

Figure 1.

Rough grain morphology features and the numerical set-up. (a) Various rough grains in the $R_{r}$ – $D_{\!f}$ space. The colour bar represents the ratio of the grain radial length $r(\theta ,\varphi )$ to its mean radial length $\overline{r}$ . (b) The schematic of a free stream passing a rough grain. Here, $\boldsymbol{p}_{{1}}$ , $\boldsymbol{p}_{{2}}$ , and $\boldsymbol{p}_{{3}}$ are directions of the longer semi-major, median and semi-minor axes of the rough grain, respectively; they are parallel with the global $x$ -, $y$ - and $z$ -axes, respectively, of which the positive directions are indicated by the arrow directions in the left bottom corner. The purple solid line denotes the semi-major grain axis after a rotation around the centre $O$ . Here, $\boldsymbol{p}_{1}{'}$ is the vector of the semi-major axis $\boldsymbol{p}_{{1}}$ after rotation, and $\boldsymbol{p}_{{1},\boldsymbol{xy}}{'}$ is its projection on the $\boldsymbol{x}$ – $y$ plane. Also, $\Theta \in [0,\pi ]$ , defined as the angle between $\boldsymbol{p}_{{3}}$ or $\boldsymbol{z}$ and $\boldsymbol{p}_{1}{'}$ , and $\Phi \in [0,2\pi )$ , defined as the angle between $\boldsymbol{p}_{{1},\boldsymbol{xy}}{'}$ and $\boldsymbol{p}_{{1}}$ or $\boldsymbol{x}$ , are the zenith angle and azimuth angle of $\boldsymbol{p}_{1}{'}$ , respectively; $\vartheta$ is the angle crossed by $\boldsymbol{p}_{1}{'}$ and the opposite direction of the flow, which is indicated by the blue arrow. (c) Mesh sensitivity study on the benchmark case, i.e. Stokes drag of a sphere. Here, $h$ is the mesh size, and $R_{0}$ is the volume-equivalent sphere radius. The selected mesh size for the systematic study is indicated by the star.

As illustrated in figure 1(a), where $c_{0}=2\sqrt{\pi }$ is set to depict a unit sphere with $\overline{r}=1$ , higher $R_{r}$ and $D_{\!f}$ could induce more irregularities in general shapes and local waviness, respectively. For natural granular materials, the value of $D_{\!f}$ for quartz sands is mostly in the range $[2.1,2.3]$ , and the value of $R_{r}$ is approximately 0.1 (Wei et al. Reference Wei, Wang, Pereira and Gan2021). Extremely fine values of $(R_{r},D_{\!f})=(0.01,2.05)$ could be encountered in artificial spherical particulates, such as glass beads; ultra-high $R_{r}\gt 0.15$ and $D_{\!f}\gt 2.4$ could appear in some macromolecules and human brains (Yotter et al. Reference Yotter, Thompson, Nenadic and Gaser2010). Without loss of generality, a series of simulations are conducted for rough grains with $R_{r}\in [0,0.35]$ and $D_{\!f}\in [2.1,2.5]$ . The classical shape parameters, including ratio ( $\mathcal{A}_{\mathscr{r}}$ ), elongation ( $\mathcal{E}_{\mathscr{i}}$ ), flatness ( $\mathcal{F}_{\mathscr{i}}$ ), roundness ( $\mathcal{R}$ ), sphericity ( $\mathcal{S}$ ) and convexity ( $\mathcal{C}_{\mathcal{X}}$ ), of these aspherical grain shapes are demonstrated in tables 1, 2, 3, 4, 5 and 6 of Appendix A. The specific definitions of these shape indices can be found in Zhao & Wang (Reference Zhao and Wang2016). Notably, $R_{r}$ cannot increase infinitely, as higher $R_{r}$ could generate negative $r(\theta ,\varphi )$ , which is contrary to the star-like shape reconstructed by SH; when expanded in the spherical coordinate system, $r(\theta ,\varphi )$ is an injective rather than bijective function. As a supplement, the ratio of minimum to maximum $r(\theta ,\varphi )$ is also provided as a newly defined aspect ratio $\mathcal{A}_{\mathcal{R}}=\min _{} (r(\theta ,\varphi ))/\max _{} (r(\theta ,\varphi ))$ in table 7 of Appendix A. For ellipsoids, $\mathcal{A}_{\mathcal{R}}=\mathcal{A}_{\mathscr{r}}$ , which is not the case for rough grains in figure 1(a).

To determine $C_{\!D}$ for a rough grain, a numerical model is developed in the next subsection. The analysis of rotational dependence for drag coefficients is performed with two rotations relative to the flow direction: azimuthal rotation $\Phi$ and polar rotation $\Theta$ as in figure 1(b). Principal component analysis (PCA) is also utilised to determine major, median and minor axes with respect to $\boldsymbol{p}_{{1}}$ , $\boldsymbol{p}_{{2}}$ and $\boldsymbol{p}_{{3}}$ (Fonseca Reference Fonseca2011), respectively. Notably, minimum and maximum $A_{\!p}$ are respectively represented by projected areas vertical to $\boldsymbol{p}_{{1}}$ and $\boldsymbol{p}_{{3}}$ . In total, 49 periodic rotations from an initially random orientation, of both $\Phi$ and $\Theta$ from $0$ to $2\pi$ with gap ${\pi }/{3}$ , plus 2 featured directions ( $\boldsymbol{p}_{{1}}$ and $\boldsymbol{p}_{{3}}$ ) are considered for each morphology. A sensitivity study of the orientation number has been performed (see figure 9(c) of Appendix B), suggesting that the statistical features of drag can be captured with this orientation number ( $N=51$ ).

2.2. Governing equations and numerical configuration

When relative movement occurs between the grain and its surrounding fluid, it experiences a resultant force. The force may not align with the flow due to the asymmetric morphology, and can be split into streamwise drag force and spanwise lift force. Here, our study exclusively focuses on the drag force. The total drag force $F_{\!D}$ is obtained by a surface integral of the total stress $\boldsymbol{\sigma }$ on the grain surface $\partial \Omega$ projected on streamwise direction $\boldsymbol{e}_{s}$ :

(2.2)

Here, $\boldsymbol{\sigma }=-p\boldsymbol{I}+\boldsymbol{\tau }$ is a sum of the pressure tensor $p\boldsymbol{I}$ and the shear stress $\boldsymbol{\tau }=[\boldsymbol{\nabla }\boldsymbol{u}+(\boldsymbol{\nabla }\boldsymbol{u})^{\rm T}]$ , where $\boldsymbol{u}$ is flow velocity. Therefore, the total drag force can be considered to be composed of two components, i.e. the skin friction and the pressure drag . We define ${C}_{\!D}^{p}$ and ${C}_{\!D}^{f}$ using ${F}_{\!D}^{p}$ and ${F}_{\!D}^{f}$ in the same manner as $C_{\!D}={F_{\!D}}/{(({1}/{2})\rho _{\!f}| \boldsymbol{u}_{\infty }| ^{2}A_{\!p})}$ , with $A_{\!p}$ being the projected area, corresponding to effects of pressure and viscous stresses on drag, respectively. For all grains in this study, ${\textit{Re}}={\rho _{\!f}}U_{\infty }d_{\textit{eq}}/\mu _{\!f}$ is approximately $10^{-5}\ll 1$ , where $\rho _{\!f}$ , $\mu _{\!f}$ and $U_{\infty }$ are respectively the density, viscosity and relative velocity of the fluid. In such a Stokesian limit, $C_{\!D}={24}/{\textit{Re}}$ . We henceforth use the Reynolds-normalised drag coefficients $C_{\overline{D}}={\textit{Re}}\,C_{\!D}$ , ${C}_{\overline{D}}^{p}={\textit{Re}}\,{C}_{\!D}^{p}$ and ${C}_{\overline{D}}^{f}={\textit{Re}}\,{C}_{\!D}^{f}$ ; for a sphere in the Stokes limit, $C_{\overline{D}}=24$ , ${C}_{\overline{D}}^{p}=8$ and ${C}_{\overline{D}}^{f}=16$ . The insensitivity of $C_{\overline{D}}$ , ${C}_{\overline{D}}^{p}$ and ${C}_{\overline{D}}^{f}$ to the selected Reynolds number ${\textit{Re}}\approx 10^{-5}$ is confirmed via simulations for spherical, ellipsoidal and rough grains, as indicated in figure 9(d) of Appendix B.

The twofold integrals in (2.2) can be obtained by solving the continuity and momentum equations. No separation of the boundary layer occurs, and the flow field can be adequately described by the Stokes equations

(2.3) \begin{equation}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}=0 \quad\text{and}\quad \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\tau }=\boldsymbol{\nabla }p.\end{equation}

In the DNS, lengths, velocities, stresses and forces are non-dimensionalised by $D_{\!p}$ , $U_{\infty }$ , ${\mu _{\!f}}U_{\infty }/d_{\textit{eq}}$ and ${\mu _{\!f}}U_{\infty }d_{\textit{eq}}$ , respectively. The flow region, as shown in figure 1(b), is a rectangular domain with $[L_{x},L_{y},L_{z}]=[18d_{\textit{eq}},20d_{\textit{eq}},20d_{\textit{eq}}]$ along, respectively, the streamwise ( $x$ -axis) and spanwise ( $y$ - and $z$ -axes) directions. A sensitivity study of the domain size was conducted on the roughest grain with $(R_{r},D_{\!f})=(0.35,2.5)$ (see figure 9(b) of Appendix B), and it is confirmed that with the selected $[L_{x},L_{y},L_{z}]$ , the domain size effect is negligible. At the inlet, outlet and sidewalls, a translational velocity $U_{{\infty }}$ along the streamwise direction is imposed. The grain remains stationary at the domain centre, and a no-slip boundary condition is imposed on the grain surface. The numerical settings are verified by the very good agreement between theoretical and simulated $C_{\!D}$ of perfectly smoothed spheres in figure 1(c).

The finite element scheme with the boundary-fitted method in COMSOL Multiphysics® is adopted to solve (2.3). The fluid domain is discretised using more than 600 000 tetrahedron elements with mixed-order shape functions, i.e. Taylor–Hood elements (Taylor & Hood Reference Taylor and Hood1973), to guarantee the divergence stability. The stationary solution of the equation set is obtained by solving a saddle point problem. For the detailed discussion of the weak form and corresponding numerical stabilisation techniques, one can refer to Hughes & Mallet (Reference Hölzer and Sommerfeld1986) and Hauke & Hughes (Reference Hauke and Hughes1994). The element size $h$ is grossly identical to that of generic triangular meshes on grain surfaces: $h=1.2\times {10^{-2}}d_{\textit{eq}}$ is calibrated using the Stokes drag and is sufficient for converged results, as confirmed in figure 1(c). The convergence is also tested in the case of the roughest grain, $(R_{r},D_{\!f})=(0.35,2.5)$ , with a relative error to the finest mesh of less than 6 % (see figure 9(a) of Appendix B), and further validated against a theoretical solution of ellipsoidal grains, as shown in the following sections.

3.1. Theoretical formula for $C_{d}$ versus $A_{\!p}$ of spheroidal grains

At a specific low ${\textit{Re}}$ , $C_{d}$ of a spheroid is an injective function of the eccentricity $e$ , i.e. its sphericity, and an inclination angle $\Theta$ (Happel & Brenner Reference Happel and Brenner1983). Since the streamwise projected area $A_{\!p}$ can be also explicitly denoted by the combination of $e$ and $\Theta$ , it is promising to transform $C_{d}(e,\Theta )$ to $C_{d}(A_{\!p})$ . A spheroidal spheroid surface can be represented in the Cartesian coordinate system as

(3.1) \begin{equation}\frac{x^{2}+y^{2}}{\mathscr{a}^{2}}+\frac{z^{2}}{\mathscr{a}^{2}\left(1-\epsilon \right)^{2}}=1,\end{equation}

where $0\lt | \epsilon | \ll 1$ . Positive $\epsilon \gt 0$ depicts an oblate spheroid with eccentricity $\sqrt{1-(1-\epsilon )^{2}}$ , while negative $\epsilon \lt 0$ depicts a prolate spheroid with eccentricity $\sqrt{(1-\epsilon )^{2}-1}$ . Notably, in Happel & Brenner (Reference Happel and Brenner1983) and many other subsequent studies on drag forces on spheroidal shapes, the aspect ratio is always defined as the length ratio of the equatorial semi–axis to the polar semi–axis, ${1}/({1-\epsilon })$ . By contrast, in this study for both ellipsoidal and rough grains, the aspect ratio is defined as the length ratio of minor axis to major axis determined by PCA, or of the minimum radial length to the maximum radial length (see Appendix A). Such a near-spherical geometry facilitates the use of analytical Stokes flow solutions, such as stream function formulations, which are otherwise intractable for strongly anisotropic shapes (Jeffery Reference Jeffery1922). According to Vickers (Reference Wadell1996), the projected area in the direction $(\ell, \mathscr{m}, \mathscr{n})=(\sin \Theta \cos \Phi , \sin \Theta \sin \Phi , \cos \Theta)$ of an ellipsoid $({x^{2}}/{\mathscr{u}^{2}})+({y^{2}}/{\mathscr{v}^{2}})+({z^{2}}/{\mathscr{w}^{2})}=1$ is $A_{\!p}=\pi \sqrt{\ell^{2}\mathscr{v}^{2}\mathscr{w}^{2}+\mathscr{m}^{2}\mathscr{u}^{2}\mathscr{w}^{2}+\mathscr{n}^{2}\mathscr{u}^{2}\mathscr{v}^{2}}$ . Thus the projected area for a spheroid described in (3.1) is

(3.2) \begin{equation}A_{\!p}=\pi \mathscr{a}^{2}\sqrt{{\left(1-\epsilon \right)^{2}}\cos ^{2} \Phi +\sin ^{2} \Phi }.\end{equation}

This equation indicates that for spheroidal shapes, $A_{\!p}$ is dependent on only one rotation angle rather than two rotation angles for triaxial-ellipsoidal or asymmetrical shapes shown in figure 1(b). Using the crosswise area of the volume-equivalent sphere $A_{\textit{eq}}=\pi \mathscr{a}^{2}(1-\epsilon )^{2/3}$ , to calculate $C_{\!D}=F_{\!D}/(({1}/{2})\rho _{\!f}| \boldsymbol{u}_{\infty }| ^{2}A_{\textit{eq}})$ instead of $A_{\!p}$ , the so-called sine-squared law exists (Ouchene et al. Reference Ouchene, Khalij, Arcen and Tanière2016; Sanjeevi & Padding Reference Sanjeevi and Padding2017):

(3.3) \begin{equation}C_{\!D}\left(\Theta \right)=C_{\!D}\left(\Theta =0\right)+\left(C_{\!D}\left(\Theta = \dfrac{\pi }{2}\right)-C_{\!D}\left(\Theta =0\right)\right) \cdot \sin ^{2} \Theta,\end{equation}

i.e.

(3.4) \begin{equation}F_{\!D}\left(\Theta \right)=F_{\!D}\left(\Theta =0\right)+\left(F_{\!D}\left(\Theta = \dfrac{\pi }{2}\right)-F_{\!D}\left(\Theta =0\right)\right) \cdot \sin ^{2} \Theta,\end{equation}

where the argument $\Theta$ denotes $C_{\!D}$ or $F_{\!D}$ of the spheroidal grain with inclined angle $\Theta$ .

For an oblate spheroid ( $0\lt \epsilon \ll 1$ ), the maximum drag force $F_{D,\textit{max}}=F_{\!D}(\Theta ={\pi }/{2})=6\pi \mu _{\!f}\mathscr{a} | \boldsymbol{u}_{\infty }| [1-\epsilon /5+\mathcal{O}(\epsilon ^{2})]$ occurs at streaming flow parallel to the long axis, while the minimum drag force $F_{D,min}=F_{\!D}(\Theta =0)=6\pi \mu _{\!f}\mathscr{a}| \boldsymbol{u}_{\infty }| [1-2\epsilon /5+\mathcal{O}(\epsilon ^{2})]$ occurs at streaming flow perpendicular to the long axis (Payne & Pell Reference Payne and Pell1960; Happel & Brenner Reference Happel and Brenner1983). If the high-order term $\mathcal{O}(\epsilon ^{2})$ is omitted, then the project-area-based drag coefficient is

(3.5) \begin{equation}C_{\overline{D}}\left(\Theta \right)={\textit{Re}}\cdot C_{D,\Theta }={\textit{Re}} \cdot \frac{F_{D,\Theta }}{\dfrac{1}{2}\rho _{\!f}\left| \boldsymbol{u}_{\infty }\right| ^{2}A_{\textit{eq}}}=24\times\dfrac{\left(1-\dfrac{2}{5}\epsilon +\dfrac{1}{5}\epsilon \sin ^{2} \Theta \right)\left(1-\epsilon \right)^{\frac{1}{3}}}{\sqrt{{\left(1-\epsilon \right)^{2}}\cos ^{2} \Theta +\sin ^{2} \Theta }}.\end{equation}

We aim to represent $C_{\overline{D},\Theta }$ as a function of $\overline{A}_{\!p}:=A_{\!p}/A_{\textit{eq}}=\sqrt{{(1-\epsilon )^{2}}\cos ^{2} \Phi +\sin ^{2} \Phi }\,(1-\epsilon )^{-2/3}$ . After representing $\cos ^{2} \Theta$ using $\overline{A}_{\!p}$ ,

(3.6) \begin{equation}\sin ^{2} \Phi =\frac{{\overline{A}}_{\!p}^{2}\left(1-\epsilon \right)^{\frac{4}{3}}-\left(1-\epsilon \right)^{2}}{1-\left(1-\epsilon \right)^{2}},\end{equation}

$C_{\overline{D}}(\Theta )$ reads as a combination function of $\overline{A}_{\!p}$ and $\epsilon$ :

(3.7) \begin{equation}C_{\overline{D}}\left(\overline{A}_{\!p},\epsilon \right)=\frac{24\left(1-\epsilon \right)^{-\frac{1}{3}}}{\overline{A}_{\!p}}\left\{1-0.4\epsilon +\frac{\dfrac{1}{5}\epsilon \left[\overline{A}_{\!p}^{2}\left(1-\epsilon \right)^{{4}/{3}}-\left(1-\epsilon \right)^{2}\right]}{2\epsilon -\epsilon ^{2}}\right\}.\end{equation}

After omitting high-orders terms ( $\mathcal{O}(\epsilon ^{2})$ , $\mathcal{O}(\epsilon ^{3})$ , and so on) via the Taylor series approximations $(1-\epsilon )^{\mathcal{N}}\approx 1-\mathcal{N}\epsilon$ , we arrive at

(3.8) \begin{equation}C_{\overline{D}}\left(\overline{A}_{\!p},\epsilon \right)\approx \frac{24}{\overline{A}_{\!p}}\left(0.9+0.1\overline{A}_{\!p}^{2}\right)=\frac{24}{\overline{A}_{\!p}}\left[1+0.1\left(\overline{A}_{\!p}^{2}-1\right)\right].\end{equation}

Let $\overline{A}_{\!p}=1+\delta$ with $| \delta | \ll 1$ . Then $\ln \overline{A}_{\!p}=\ln (1+\delta )\approx \delta$ by Taylor’s formula, and only consideration of $\mathcal{O}(\delta )$ ; under the same condition, $\overline{A}_{\!p}^{2}=(1+\delta )^{2}\approx 1+2\delta \approx 1+2\ln \overline{A}_{\!p}\rightarrow \overline{A}_{\!p}^{2}-1\approx 2\ln \overline{A}_{\!p}$ . It is also well known that $\text{e}^{x}=\lim _{n\rightarrow {\infty }} (1+({x}/{n}))^{n}$ if $x\approx 0$ and $\text{e}^{x}\approx 1+x$ . Since $\epsilon \approx 0\rightarrow \ln \overline{A}_{\!p}\approx 0$ , we get $\text{e}^{0.2\ln \overline{A}_{\!p}}=\overline{A}_{\!p}^{0.2}\approx 1+0.2\ln \overline{A}_{\!p}\rightarrow 2\ln \overline{A}_{\!p}\approx \overline{A}_{\!p}^{2}-1\approx 10(\overline{A}_{\!p}^{0.2}-1)$ . Taking $\overline{A}_{\!p}^{2}-1\approx 10(\overline{A}_{\!p}^{0.2}-1)$ in (3.8), we have

(3.9) \begin{equation}C_{\overline{D}}\left(\overline{A}_{\!p}\right)\approx 24\overline{A}_{\!p}^{-0.8}.\end{equation}

Notably, when $-1\ll \epsilon \lt 0$ , (3.1) represents the surface of a prolate spheroid, and (3.3)–(3.9) still hold for the scaling relationship $C_{\overline{D}}\approx 24\overline{A}_{\!p}^{-0.8}$ . In fact, both oblate and prolate spheroidal grains are special cases of triaxial-ellipsoidal grains. Compared to one-angle-dependent $C_{\overline{D}}(\Theta )$ of spheroidal grains, the drag coefficient $C_{\overline{D}}(\Theta ,\Phi )$ of triaxial-ellipsoidal grains is two-angle-dependent. Albeit such a big difference, $C_{\overline{D}}\propto \overline{A}_{\!p}^{-0.8}$ is also the case. For the better flow of logic, we put the derivation process in Appendix C.

The derived power law $C_{\overline{D}}\propto \overline{A}_{\!p}^{-0.8}$ is validated in figure 2 by the theoretical data points of both spheroidal and ellipsoidal grains. The more the grain shape departs from a perfect sphere, the more its data deviate from $C_{\overline{D}}\approx 24\overline{A}_{\!p}^{-0.8}$ ; however, even when the short axis length is only half of that of the long axis for spheroidal grains, $C_{\overline{D}}$ can be well linear to $\overline{A}_{\!p}^{-0.8}$ on a log-log scale. For brevity, we use only one parameter, $\varepsilon$ , to define the triaxial-ellipsoidal surface ${x^{2}}/{a^{2}}+{y^{2}}/({a^{2}(1-\varepsilon))}+{z^{2}}/({a^{2}(1-\varepsilon)^{2})}=1$ in figure 2(b). The data appear to be logarithmically linear for lower $\varepsilon$ . Interestingly, the relation between $C_{\overline{D}}$ and $\overline{A}_{\!p}$ becomes gradually not injective with increasing $\varepsilon$ , as illustrated in the left-hand inset of figure 2(b). Albeit the emerging $(C_{\overline{D}},\overline{A}_{\!p})$ data form a scattered band rather than a clear line, the $C_{\overline{D}}$ – $\overline{A}_{\!p}$ relation can be well approximated by a log-linear function. Motivated by the high efficiency of this power law in ellipsoidal grains, we use obtained DNS results to examine whether this power law is suitable for rough grains, for which it is challenging to derive $C_{\overline{D}}$ .

Figure 2.

Analytical results of $C_{\overline{D}}$ versus $\overline{A}_{\!p}$ for (a) spheroidal and (b) ellipsoidal grains. The black stars represent the sphere data. The colour bar is shared by both $\epsilon \in [-0.5,0.5]$ and $\varepsilon \in [0.5,1]$ . The solid line denotes the fitting line via the formula $C_{\overline{D}}\propto \overline{A}_{\!p}^{-0.8}$ . Zoom-ins in all plots are for enlargements around the point $(\overline{A}_{\!p},C_{\overline{D}})=(1,24)$ . In (a), the left bottom inset with $0\lt \epsilon \lt 0.5$ is for oblate spheroids marked by the square, while the right upper inset with $-0.5\lt \epsilon \lt 0$ is for prolate spheroids marked by the circle.

3.2. Universality of $C_{\overline{D}}\propto \overline{A}_{\!p}^{-0.8}$

Mutual verifications between the derived power law $C_{\overline{D}}\propto \overline{A}_{\!p}^{-0.8}$ and DNS are demonstrated in figure 3(a), for which the fitting parameters are listed in Appendix D. The database contains one spherical, six oblate and five prolate spheroidal, and four triaxial-ellipsoidal grains, of which shape indices are summarized in table 8. The aspect ratio $\mathcal{A}_{\mathcal{R}}$ , defined as the ratio of the minimum to maximum radial length (see Appendix A), of oblate spheroids is $[0.1,0.9]$ , while the range of prolate spheroids is $[0.59,0.91]$ . Four triaxial ellipsoids are with the same $\mathcal{A}_{\mathscr{r}}$ , defined as the minimum to the maximum principal dimension in range $[0.51,0.97]$ , for the rough particles $(R_{r},D_{\!f})\in \{(0.2,2.2),(0.2,2.5),(0.35,2.1),(0.35,2.5)\}$ . All ellipsoidal grains show good alignment with the line $C_{\overline{D}}=24\overline{A}_{\!p}^{-0.8}$ .

Figure 3.

Log-log plots of $C_{\overline{D}}$ versus $\overline{A}_{\!p}$ for (a) the spherical, spheroidal and triaxial-ellipsoidal cases and (b–f) each $R_{r}$ . The dashed lines represent the fitting results using $C_{\overline{D}}\propto \overline{A}_{\!p}^{-0.8}$ . All fittings in (b–f) have R-squared values larger than 0.90, which are summarised in table 11 of Appendix D. The colour bar is shared by $\epsilon \in [-0.7,0.9]$ , $\mathcal{A}_{\mathcal{R}}\in [0.5,0.8]$ and $D_{\!f}\in [2.1,2.5]$ . The black star marks the value of a sphere. The $x$ -axis range in all plots is $[0.9\min _{} \overline{A}_{\!p},1.1\max _{} \overline{A}_{\!p}]$ . Zoom-ins in all plots are for enlargements around the point $(\overline{A}_{\!p},C_{\overline{D}})=(1,24)$ .

Figure 4.

(a) Comparisons of $C_{\overline{D}}$ and (b) the ratio of pressure drag to friction drag between DNS results and the deterministic model from Leith (Reference Leith1987). In both plots, dashed lines represent the equality $y=x$ . Solid markers denote cases where the flow is perpendicular to the projected area in the plane spanned by $\boldsymbol{p}_{{1}}$ – $\boldsymbol{p}_{{2}}$ or $\boldsymbol{p}_{{2}}$ – $\boldsymbol{p}_{{3}}$ . The red star is for the spherical case. Zoom-ins in all plots are for enlargements around the data point for the spherical case.

In line with testing whether analytics from minimal geometries apply to more complex shapes, the model of Leith (Reference Leith1987) may be the most accepted. Therein, it is postulated that since for a sphere in the Stokesian regime, two-thirds of the drag is due to viscous effects and one-third to the pressure drag, the same should hold for aspherical objects. With further consideration of $F_{\!D}$ as the integrations of fluid pressure and friction over the object surface, the former and the latter should be associated with the area projected normal to the fluid motion and the area of the total surface, respectively:

(3.10) \begin{equation}F_{\!D}=6\pi \mu \left| \boldsymbol{u}_{\infty }\right| \left(\underset{\text{pressure}}{\underbrace{\frac{1}{3}R_{{A_{\!p}}}}}+\underset{\text{viscosity}}{\underbrace{\frac{2}{3}R_{{S_{a}}}}}\right)\!,\end{equation}

where $R_{{A_{\!p}}}=\sqrt{{A_{\!p}}/{\pi }}$ and $R_{{S_{a}}}=\sqrt{({S_{a}}/{4\pi }})$ are, respectively the radius lengths of the projected and surface-area-equivalent spheres. Importing this equation to measure the drag coefficient, we get

(3.11) \begin{equation}C_{\!D}=\dfrac{24\pi \left(\dfrac{1}{3}R_{{A_{\!p}}}+\dfrac{2}{3}R_{{S_{a}}}\right)R_{\textit{eq}}}{Re\,A_{\!p}}\quad\rightarrow\quad C_{\overline{D}}=\frac{\left(8R_{{A_{\!p}}}+16R_{{S_{a}}}\right)R_{\textit{eq}}}{R_{{A_{\!p}}}^{2}},\end{equation}

and the ratio of pressure to friction drag is $R_{A_{p}}/(2R_{S_a})$ . Interestingly, as demonstrated in figure 4(a), except for the spherical case, results predicted by Leith’s model systematically depart from those from the DNS; this is because the pressure– and friction–drag contributions are nearly independent, as shown in figure 4(b). The ratio of the two departs well from the postulated value ${R_{{A_{\!p}}}}/({2R_{{S_{a}}}})$ . It is worth noting that to further improve the accuracy of Leith’s model in predicting $C_{\overline{D}}$ of realistic shaped grains, Hölzer & Sommerfeld (Reference Hughes and Mallet2008) utilised the lengthwise sphericity instead of the crosswise sphericity in the original model to account for the pressure part in total drag. Compared with the latter, the former is more complicated and involves in quantifying the average cross-sectional area parallel to the flow direction. This replacement can only slightly improve the predictive accuracy in the rough drag, as demonstrated in figure 10(a) of Appendix E for the roughest grain $(R_{r}=0.35,\ D_{\!f}=2.5)$ ; it is less robust when compared with $C_{\overline{D}}\propto \overline{A}_{\!p}^{-0.8}$ . This indicates the unreasonableness of directly splitting the overall drag by the fashion of Stokes’ law, when coping with the interaction between a rough grain and the fluid around it.

Then we further test our power law $C_{\overline{D}}\propto \overline{A}_{\!p}^{-0.8}$ against the asymmetrical rough grains. From figures 3(b–f), it can be seen that when $\overline{A}_{\!p}$ values are within $[0.9,1.8]$ for all rough grains shown in figure 1(a), the power law still holds. However, the scale factor $a$ in $C_{\overline{D}}=a\overline{A}_{\!p}^{-0.8}$ is no longer 24, and depends on the grain morphology (see table 10 in Appendix D). Generally, the rougher a grain is, the smaller the scale factor may become, though it is not strictly monotonic with $R_{r}$ or $D_{\!f}$ , and the smallest $a$ appears at $R_{r}=0.20,\ D_{\!f}=2.5$ . More importantly, for each rough grain considered, most orientations result in $C_{\overline{D}}$ values lower than that of the spherical case, indicating a net drag reduction, as shown in figures 3(b–f). Roughness causes concave surficial features, which are in nature different from purely convex geometries, such as ellipsoids. It apparently leads to two factors contributing to this drag reduction according to $C_{\overline{D}}=a\overline{A}_{\!p}^{-0.8} : \overline{A}_{\!p}$ is mostly smaller than 1, and $a$ is smaller than 24. In the following sections, after discussing which is the dominant role between pressure and friction drag in drag reduction, more evidence from micromechanical analysis will be provided.

Taken together, the simulations calibrate the 0.8 power law over a finite morphological and geometrical range. For the considered ellipsoids, the aspect ratio $\mathcal{A}_{\mathcal{R}}$ (defined as the ratio of minimum to maximum radial length; see Appendix A) spans $[0.1,0.9]$ for oblate spheroids, and $[0.59,0.91]$ for prolate spheroids, and the PCA-determined equivalent-ellipsoid aspect ratio $\mathcal{A}_{\mathscr{r}}$ of triaxial ellipsoids lies in $[0.51,0.97]$ . The rough grains span $R_{r}\in [0.01,0.35]$ and $D_{\!f}\in [2.1,2.5]$ . For all these shapes with various orientations such that $\overline{A}_{\!p}\in [0.9,1.8]$ , the Reynolds-number-normalised drag coefficient can be well described by $C_{\overline{D}}=\alpha \overline{A}_{\!p}^{-0.8}$ with R-squared values much higher than 0.9, as shown in figure 3 and Appendix D. Outside such an interval for $\overline{A}_{\!p}$ , system deviations from the fit are encountered. Therefore, the scaling law should not be used for quantitative extrapolation there, such as the extremely high or low values of $\overline{A}_{\!p}$ for an oblate spheroid with $\mathcal{A}_{\mathcal{R}}=0.1$ , as demonstrated in figure 3(a).

We have also benchmarked our 0.8 power law against representative classical drag correlations for non-spherical particles, summarised in table 12 in Appendix E. The main outcome is that the present projected-area-based closure remains competitive for axisymmetric bodies in the creeping-flow limit, and becomes particularly advantageous for two-angle-dependent grains, such as triaxial ellipsoids and rough grains. The detailed comparison is deferred to § 6.1, with the indicative proof in figure 10 in Appendix E.

Figure 5.

Survival probability function $P_{s}$ of (a–e) $C_{\overline{D}}$ , (f–j) ${C}_{\overline{D}}^{f}$ and (k–o) ${C}_{\overline{D}}^{p}$ at varying rotation angles for each grain morphology in the $R_{r}$ – $D_{\!f}$ space. The columns from left to right correspond to $R_{r}$ values 0.01, 0.05, 0.10, 0.20 and 0.35. In each plot, the solid mark represents data of the flow perpendicular to the projected area on the plane crossed by $\boldsymbol{p}_{{1}}$ and $\boldsymbol{p}_{{2}}$ or $\boldsymbol{p}_{{2}}$ and $\boldsymbol{p}_{{3}}$ . The solid line represents the fitting using the survival probability function of the Weibull distribution. The dashed line denotes values of the spherical case.

6.1. Comparison with existing drag correlations

Concurrently, existing models accounting for the orientation dependence of drag coefficients could be classified into two categories. The first exploits the linearity of the Stokes equations and leads to orientation laws of the Happel–Brenner type for axisymmetric bodies, in which the drag varies between principal orientations through a one-angle closure. Such formulations, together with their subsequent extensions, perform well for oblate and prolate spheroids. Therefore, the present law should not be viewed as a competing asymptotic theory in that setting. Rather, for smooth ellipsoids, it is a compact reformulation of the same creeping-flow structure in terms of the normalised projected area. This explains why the present $C_{\overline{D}}=a\overline{A}_{\!p}^{-0.8}$ relation remains competitive for one-angle-dependent bodies in the Stokes limit.

The distinction becomes more important for bodies for which drag depends on two rotation angles. Triaxial ellipsoids and rough grains cannot be represented naturally by one-angle closures, whereas the projected-area-based formulation remains applicable, because orientation enters through the measurable quantity $\overline{A}_{\!p}$ . By contrast, widely used empirical or semi-empirical correlations, such as Ganser (Reference Ganser1993) and related models for irregular particles, were designed for broader Reynolds number ranges and more general shape classes. That broader generality is valuable outside the present regime, but it also tends to compromise accuracy in the creeping-flow DNS benchmark considered here. Accordingly, the main contribution of the present law is not that it universally outperforms every existing correlation, but that it provides a compact Stokes regime closure that is competitive for spheroids and, more importantly, extends naturally to two-angle-dependent grains as demonstrated in figure 10 in Appendix E.

6.2. Rationale for the projected-area-based power-law closure

The adoption of a power-law closure in $\overline{A}_{\!p}$ is motivated by both theory and parsimony. For smooth ellipsoids in the Stokes limit, combining the resistance-tensor formulation with the projected-area relation reduces the orientation dependence of drag to a simple function of the normalised projected area, which is well approximated by $C_{\overline{D}}\propto \overline{A}_{\!p}^{-0.8}$ . The DNS results show that this exponent remains robust not only for spheroidal and triaxial ellipsoids, but also for rough grains, while morphology primarily modifies the prefactor $a$ . In this sense, the proposed law is not introduced as an ad hoc empirical fit; rather, it is a compact geometry-driven rewriting of the creeping-flow drag dependence.

Alternative correlations for irregular or non-spherical particles were not adopted as the main parametrisation for three reasons. First, many such models were calibrated over substantially wider Reynolds number ranges, and were not developed specifically for the Stokes regime considered here. Second, they often rely on multiple global shape descriptors, for which connection to instantaneous orientation is indirect. Third, several available closures are tailored to one-angle-dependent bodies and therefore do not transfer naturally to triaxial ellipsoids or rough grains. To fill the gap, our objective is to identify a simple closure for isolated grains in creeping flow, in which the instantaneous projected area plays the central role. This makes the model easy to evaluate from imaging or equivalent-ellipsoid reconstruction, and naturally compatible with orientation-distribution-based predictions.

6.3. Applicability, assumptions and limitations of the power law

The present closure is derived and validated for isolated grains immersed in unbounded Newtonian fluids in the creeping-flow limit. Further, during the interaction between fluids and the grain, the grain is held stationary at a fixed orientation. The DNS database is centred on an extremely small Reynolds number, and the supplementary tests indicate that the normalised drag coefficients remain essentially insensitive to ${\textit{Re}}$ up to the Stokes range upper limit examined here. Similar to the high- ${\textit{Re}}$ flow in a wall-bounded rough channel, roughness could reduce drag over a rough grain at creeping flow. The turbulence in fully rough channel flow at high ${\textit{Re}}$ is governed by inertial effects, inter-element interactions, roughness-layer development, and sheltering phenomena absent from our Stokes regime simulations of an isolated grain. Our power law and distributional results are therefore not intended for quantitative prediction in rough-wall flows. They do, however, provide qualitative micromechanistic insight: roughness tends to decrease the friction-drag component more than the pressure-drag component, and orientation controls the variability via projected area. Any quantitative mapping to rough-wall parameters (e.g. $k_{s}$ ) requires dedicated high- ${\textit{Re}}$ studies that include the above mechanisms. In particular, although ${\textit{Re}}\approx 10^{-5}$ is selected for all cases, as confirmed in figure 9(d), within the Stokes range explored there ( ${\textit{Re}}\in [10^{-4},0.1]$ ) for extra simulations of three distinctive grain shapes at a certain rotation, the normalised metrics $C_{\overline{D}}$ , $C_{\overline{D}}^{f}$ and $C_{\overline{D}}^{p}$ are insensitive to ${\textit{Re}}$ , so the present law should be viewed as a Stokes regime closure to be augmented by explicit inertial corrections when extending to finite ${\textit{Re}}\gt 0.1$ . For quantitative applications at higher ${\textit{Re}}$ , or in the presence of walls and particle–particle interactions, the $\overline{A}_{\!p}^{-0.8}$ dependence should be combined with explicit inertial and confinement corrections to the prefactor $\alpha$ ; outside the calibrated ranges of shape ( $\overline{A}_{\!p}\in [0.9,1.8]$ ) and ${\textit{Re}}\lt 0.1$ , our results should therefore be interpreted as qualitative guidance rather than a universal law.

In addition, the present framework does not include near-wall effects, particle–particle interactions, non-Newtonian rheology, or unsteady coupling between translation and rotation. Finally, when the closure is used to predict ensemble drag statistics in settling or transport problems, the orientation distribution must be supplied by experiment, DNS or an auxiliary dynamical model; the projected-area law maps orientation statistics into drag statistics, but does not by itself predict the orientation dynamics.

Figure 9.

Sensitivity study of (a) mesh size, (b) lateral domain size, (c) orientation number, and (d) particle Reynolds number, for (a,b) the roughest case, (c) the rough grain with $D_{\!f}=2.5$ and $R_{r}=0.2$ , and (d) the spherical, ellipsoidal ( $\mathcal{A}_{\mathscr{r}}=0.6023$ ) and roughest ( $D_{\!f}=2.5$ and $R_{r}=0.35$ ) grains at one specific rotation angle.

6.4. Physical implications: settling experiments and micromechanics

Revisiting settling-grain-experiment-based empirical equations to predict $C_{\overline{D}}$ of irregular grains at creeping flows (Michaelides & Feng Reference Michaelides and Feng2023), nearly all of them have the same general formulation: $C_{\overline{D}}=C_{\overline{D}0}(1+\mathcal{F})$ , with $\mathcal{F}\gt 0$ being a function of shape indices. The drag enhancement in experiments is contradictory to the prevalent view – surface roughness equals drag reduction – in both rough channel flow (Choi et al. Reference Chéron, Evrard and van Wachem2008) and simulations of single rough grains with various rotations. Focusing on the single-grain drag, the analytical derivations of Leith (Reference Leith1987) are also contrary to the general view, $C_{\overline{D}}=C_{\overline{D}0}(1+\mathcal{F})$ . As demonstrated in figure 4(a), higher roughness could induce a net decrease of $C_{\!D}$ , although Leith’s model systematically departs from DNS for non-spherical shapes. Appendix E benchmarks this 0.8 power law against widely used correlations for axisymmetric and irregular shapes; within the Stokes limit, it performs competitively for spheroids (one-angle dependence) and, crucially, generalises to two-angle-dependent bodies (triaxial ellipsoids and rough grains), where many one-angle closures deteriorate.

To explain this paradox, we start with examining the settling grain in a static fluid, concentrating on two parameters: the projected area $A_{\!p}$ , and the steady settling down velocity $w_{\!f}$ . The experimental drag coefficient is determined by $C_{\!D}=2(\rho _{s}-\rho _{\!f})gV_{\!p}/(\rho _{\!f}A_{\!p}w_{\!f}^{2})$ , where $\rho _{s}$ is the density of the grain. First, although rotated less or more periodically, the gravity-driven grain tends to fall with its maximum $A_{\!p}$ of higher likelihood to be perpendicular to the relative movement direction (Dietrich Reference Dietrich1982). This manifestation could induce lower $w_{\!f}$ of the grain. Note that although in such a type of experiment $w_{\!f}$ is quantified by the vertically falling down velocity, $w_{f,z}\approx w_{\!f}$ , in fact the grain could dive spirally, and $w_{\!f}=\sqrt{w_{f,x}^{2}+w_{f,y}^{2}+w_{f,z}^{2}}\gt w_{f,z}$ . As demonstrated in a recent study on an experiment on the settling motion of coral grains (Chen et al. Reference Chen2024), the spiral trajectory becomes more tumbling and even chaotic from an elliptical shape in the $x$ – $y$ plane with increasing grain irregularity. Further, in their experiments the lateral velocity $\sqrt{w_{f,x}^{2}+w_{f,y}^{2}}$ could be enlarged without an asymptotic value. This indicates that using the approximation $w_{\!f}\approx w_{f,z}$ could highly underestimate the falling down velocity for grains with higher roughness. Second, quantifying the inconsistent values of $A_{\!p}$ is non-trivial during the falling down process of the grain with quasi-periodic rotations. Considering the ease of measuring grain volume by $V_{\!p}={m_{s}}/{\rho _{s}}$ with the grain mass $m_{s}$ , the projected area $A_{\textit{eq}}$ of volume-equivalent sphere is always implemented instead. Such a value can be lower by up to 80 % of the maximum $A_{\!p}$ , as illustrated in figure 6. The significant underestimations of both $w_{\!f}$ and $A_{\!p}$ are attributed to the overestimation of $C_{\!D}$ in this series of experiments. The present framework yields testable predictions for such experiments: (i) if the orientation distribution during fall is skewed towards large $\overline{A}_{\!p}$ , then the measured $C_{\overline{D}}$ distribution should match the Fréchet-type upper tail implied by (4.2); and (ii) measuring the distribution of $\overline{A}_{\!p}$ (e.g. via imaging or equivalent-ellipsoid approximation) suffices to generate an a priori prediction of the full $C_{\overline{D}}$ distribution for a given morphology via the 0.8 power mapping.

Intuitively, one may ask why the general formula $C_{\overline{D}}=C_{\overline{D}0}(1+\mathcal{F})$ and its inference – roughness means drag enhancement – are so popular. Although the mean value is negatively related to $R_{r}$ , the standard deviation even vanishes at extremely low $R_{r}$ in our DNS with fixed grain rotation; this facilitates the possibility where $C_{\overline{D}}$ of high $R_{r}$ could be larger than that of low $R_{r}$ , resulting in roughness-induced drag enhancement in experiments. One may argue that in the situation of complete drag reduction, i.e. for extremely high $R_{r}$ and $D_{\!f}$ , $C_{\overline{D}}$ is always smaller than $C_{\overline{D}0}$ , and no drag improvement could be encountered. Considering the narrow ranges of natural grains ( $R_{r}\lt 0.1$ and $D_{\!f}\lt 2.3$ ), for which settling experiments are conducted and are mostly from the geoscience community for sedimentary study, complete drag reduction is less likely to happen. To bypass the large uncertainty in grain rotation during settling, we propose a power law $C_{\overline{D}}=a\overline{A}_{\!p}^{-0.8}$ to introduce the rotational dependence of drag coefficient for the single aspherical grain of various rotation angles.

At the grain scale, the DNS results indicate that roughness-induced drag reduction originates mainly from the friction-drag contribution. Local ridges intensify shear only over limited surface patches, whereas dimples and concave regions broaden low-shear zones over a larger fraction of the surface; the net effect is therefore a downward shift of the friction-drag contribution, while the pressure-drag contribution remains comparatively insensitive to morphology. This interpretation also explains why the prefactor $a$ decreases as grains depart from the purely convex limit. Increasing concavity weakens the effective skin-friction contribution without requiring a comparable reduction in pressure drag.