2.1. Governing equations

Conservation of momentum and mass for the fluid phase is expressed in terms of incompressible Navier–Stokes and continuity equations for a Newtonian fluid as follows:

(2.1)\begin{gather} \frac{\partial{\boldsymbol{u}}}{\partial{t}} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla}{\boldsymbol{u}} = -\boldsymbol{\nabla}{p} + \frac{1}{\textit{Re}} \nabla^2 \boldsymbol{u} , \end{gather}

(2.2)\begin{gather} \boldsymbol{\nabla}{\boldsymbol{\cdot}}\, \boldsymbol{u} = 0, \end{gather}

where $\boldsymbol {u}$, $p$ and $\textit {Re}$ respectively denote the fluid velocity vector, pressure and particle Reynolds number, defined as

(2.3)\begin{equation} \textit{Re} = \frac{\rho U d}{\mu} = \frac{\rho (1 - \phi) u_s d}{\mu}, \end{equation}

where $\rho$, $\mu$ are the fluid phase density and dynamic viscosity. The particle Reynolds number is defined based on the superficial velocity $U = (1 - \phi ) u_s$, and $u_s$ represents the average interstitial fluid velocity. In (2.1) and (2.2) and what follows, all variables are non-dimensionalised with respect to the particle diameter $d$ as the length scale, $U$ as the velocity scale, $\rho U^2$ as the pressure scale and $\rho U^2 d^2$ as the force scale. The hydrodynamic force and torque exerted on each particle denoted by $\boldsymbol {F}$ and $\boldsymbol {T}$ are given as

(2.4a)\begin{gather} \boldsymbol{F}= \int_{S}\left[ -p\boldsymbol{\mathsf{I}}+\frac{1}{\textit{Re}}(\boldsymbol{\nabla}{\boldsymbol{u}}+\boldsymbol{\nabla}{\boldsymbol{u}^{\textrm{T}}}) \right] \boldsymbol{\cdot} \boldsymbol{n} \,\textrm{d}S, \end{gather}

(2.4b)\begin{gather} \boldsymbol{T} = \int_{S}\boldsymbol{r}\times\left[ -p\boldsymbol{\mathsf{I}}+\frac{1}{\textit{Re}}(\boldsymbol{\nabla}{\boldsymbol{u}}+\boldsymbol{\nabla}{\boldsymbol{u}^{\textrm{T}}}) \right] \boldsymbol{\cdot} \boldsymbol{n} \,\textrm{d}S, \end{gather}

with $S$ denoting the surface enclosing the solid body, $\boldsymbol{\mathsf{I}}$ the identity matrix, $(\cdot )^\textrm {T}$ the matrix transpose, $\boldsymbol {n}$ the unit vector normal to the boundary of the solid body and $\boldsymbol {r}$ the position vector relative to the particle centre of mass.

2.2. Numerical method

As our PR-DNS tool, we use PeliGRIFF (Parallel Efficient Library for GRains in Fluid Flow) which is a multiphase flow solver based on the distributed Lagrange multiplier-fictitious domain (DLM-FD) formulation proposed by Glowinski et al. (Reference Glowinski, Pan, Hesla and Joseph1999). In our implementation, we employ a finite-volume staggered-grid scheme for the fluid conservation equations (Wachs et al. Reference Wachs, Hammouti, Vinay and Rahmani2015). Similar to the IBM, particles are immersed in the fluid domain in the DLM-FD method and rigid-body constraints on the fictitious fluid inside the solid region are enforced through a set of Lagrange multipliers collocated in the particle domain. For a fixed array of particles, the combined momentum and continuity equations in a non-variational form are given as

(2.5a)\begin{gather} \frac{\partial{\boldsymbol{u}}}{\partial{t}} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla}{\boldsymbol{u}} = -\boldsymbol{\nabla}{p} + \frac{1}{\textit{Re}} \nabla^2\boldsymbol{u} - \boldsymbol{\lambda} \quad \text{in} \ \mathbb{D}, \end{gather}

(2.5b)\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{u}} = 0 \quad\text{in} \ \mathbb{D}, \end{gather}

(2.5c)\begin{gather} \boldsymbol{u} = \boldsymbol{0} \quad\text{in} \ \mathbb{P}, \end{gather}

where the solid domain is denoted by $\mathbb {P}$ and the fluid–particle domain by $\mathbb {D}$. Furthermore, $\boldsymbol {\lambda }$ shows the distributed Lagrange multiplier vector which is used to enforce the rigid-body motion constraint (2.5c). For the temporal discretisation, we employ a first-order Marchuk–Yanenko operator-splitting algorithm. At each time $t^{n+1}$, we solve:

1. (i) A classical projection scheme for the solution of the Navier–Stokes problem: find $\boldsymbol {u}^{n+1/2}$ and $p^{n+1}$ such that

   (2.6a)\begin{align} &\frac{\tilde{\boldsymbol{u}}^{n+1/2} - \boldsymbol{u}^{n}}{{\rm \Delta} t} - \frac{1}{2 \textit{Re}}\boldsymbol{\nabla}^{2}{\tilde{\boldsymbol{u}}^{n+1/2}}\nonumber\\ &\quad =-{\boldsymbol{\nabla} p^{n}} + \frac{1}{2 \textit{Re}}{\nabla^{2} \boldsymbol{u}^{n}} -\frac{1}{2} \left( 3\boldsymbol{u}^{n}\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}^{n} - \boldsymbol{u}^{n-1}\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}^{n-1} \right) - \boldsymbol{\lambda}^{n},\end{align}
   (2.6b)\begin{gather} \left.\begin{array}{c@{}} \nabla^{2}{\psi^{n+1}} = \dfrac{1}{{\rm \Delta} t}\boldsymbol{\nabla}\boldsymbol{\cdot}{\tilde{\boldsymbol{u}}^{n+1/2}},\\ \dfrac{\partial{\psi^{n+1}}}{\partial{n}} =0\quad \text{on} \ \partial \mathbb{D}, \end{array}\right\} \end{gather}
   (2.6c)\begin{gather} \boldsymbol{u}^{n+1/2} = \tilde{\boldsymbol{u}}^{n+1/2} - {\rm \Delta} t \boldsymbol{\nabla}\psi^{n+1}, \end{gather}
   (2.6d)\begin{gather} p^{n+1} = p^n + \psi^{n+1} - \frac{{\rm \Delta} t}{2 \textit{Re}}\nabla^{2}{\psi^{n+1}}. \end{gather}

2. (ii) A fictitious domain problem: find $\boldsymbol {u}^{n+1}$ and $\boldsymbol {\lambda }^{n+1}$ such that

   (2.7a)\begin{gather} \frac{\boldsymbol{u}^{n+1}-\boldsymbol{u}^{n+1/2}}{{\rm \Delta} t} +\boldsymbol{\lambda}^{n+1} = \boldsymbol{\lambda}^{n}, \end{gather}
   (2.7b)\begin{gather} \boldsymbol{u}^{n+1} = \boldsymbol{0} \quad \text{in} \ \mathbb{P}, \end{gather}

where ${\rm \Delta} t$ denotes the time step, $\psi$ the pseudo-pressure and $\partial \mathbb {D}$ the domain boundary. In (2.6a), second order in time Crank–Nicolson and Adams–Bashforth schemes are used to discretise the viscous and advective terms, respectively, and the saddle-point problem in step (ii) is handled by an Uzawa algorithm (Wachs et al. Reference Wachs, Hammouti, Vinay and Rahmani2015). Considering the high-order correction of the pressure, the projection scheme in step (i) is also second-order accurate in time. However, the first-order time discretisation of the fictitious domain sub-problem in step (ii) and the first-order Marchuk–Yanenko method reduce the global time accuracy of our algorithm to first order only. Equations presented in step (i) are spatially discretised with a second-order central scheme for the diffusion term, whereas the advective term is treated with a total variation diminishing scheme combined with a Superbee flux limiter. Despite the second-order discretisation of the flow equations, the accuracy of our method is between first and second order due to the presence of rigid bodies immersed within the domain. It can be shown that with the DLM-FD method, the Lagrange multiplier $\boldsymbol {\lambda }$ can be directly integrated over the volume of each particle to obtain the hydrodynamic force and torque acting on the particle $\mathbb {P}$

(2.8a)\begin{gather} \boldsymbol{F} = \int_{\mathbb{P}} \boldsymbol{\lambda} \,\textrm{d}\boldsymbol{x}, \end{gather}

(2.8b)\begin{gather} \boldsymbol{T} = \int_{\mathbb{P}} \boldsymbol{r} \times \boldsymbol{\lambda} \,\textrm{d}\boldsymbol{x}. \end{gather}

2.3. Simulation set-up

A summary of all considered cases is given in table 1. For our PR-DNS simulations, we consider triply periodic cubic domains of edge length $L$ containing $N_p$ spherical particles each taking up a volume of $v_p = \pi / 6$, corresponding to a solid volume fraction of $\phi = N v_p/ L^3$. A constant flow rate so as to attain the desired Reynolds number is imposed in the $x$ direction using a dynamically adjusted pressure drop. The $x$ direction hence corresponds here to the streamwise direction. As noted by Akiki et al. (Reference Akiki, Jackson and Balachandar2016), Tenneti et al. (Reference Tenneti, Garg and Subramaniam2011) demonstrated that using a domain size of only $2.4 d$ guarantees the decorrelation of fluid velocities for $(\textit {Re}, \phi ) = (20, 0.2)$ and $(\textit {Re}, \phi ) = (300, 0.2)$. Our computational domains given in table 1 all extend far beyond $2.4 d$ containing $\approx$2500–3000 particles each, ensuring both the decay of fluid correlations and statistical reliability. Initialisation of particle locations for cases where $\phi \in \{0.1, 0.2\}$ is performed by distributing spheres randomly in the domain according to a random number generation algorithm without allowing any overlap. For the highly dense cases with $\phi = 0.4$, we start with a structured array of particles where each particle is given a random translational and rotational velocity, and we let the system reach an asymptotic motionless state through dissipative collisions. For all solid volume fractions, pair correlation functions have been obtained and verified to be satisfactorily matching theoretical radial distribution functions of hard spheres (Percus & Yevick Reference Percus and Yevick1958; Wertheim Reference Wertheim1963). Visualisations of two sample cases in the present study are shown in figure 2. The time step is taken to be smaller than ${\rm \Delta} t = 2 \times 10^{-3}$ in all cases, ensuring time accuracy of the simulations along with satisfying the Courant–Friedrichs–Lewy (CFL) condition; namely, $\mathrm {CFL} < 0.4$ for the spatial resolutions presented in table 1. The simulations have been run until steady state is achieved, and force and torque data are all collected from the steady part. For $\textit {Re} = 150$ at all volume fractions, the flow becomes time dependent and oscillating. This occurs because particles in close proximity act as a single obstacle and hence increase the effective length scale and the Reynolds number. Data from these cases are collected after a statistically pseudo-stationary state is established by averaging over sufficiently extended time intervals. In terms of computational demand, the mesh resolution and the domain size in each case along with the load per CPU core dictate resource requirements. With each processor core handling $512 \times 10^3$ grid cells, the smallest simulations with $\approx \! 110 \times 10^6$ grid cells were run on 192 cores, whereas the largest ones with $512 \times 10^6$ grid cells required 960 cores. All computations were carried out on groups of 48-core nodes each equipped with $\approx$187 GB of memory, provided by the Cedar supercomputer as a part of Compute Canada's advanced research computing infrastructure (http://www.computecanada.ca).

Figure 2.

Fixed beds of spherical particles at $\textit {Re}=40$ and solid volume fractions of $\phi = 0.1$ (a) and $\phi = 0.4$ (b). The streamlines shown are coloured with respect to the fluid velocity magnitude.

Table 1.

Summary of the parameters used for PR-DNS of fixed beds of spheres.

2.4. Validation

The aim in this section is to provide a comparison of the results of our PR-DNS simulations with the reported drag correlations in the literature (Beetstra et al. Reference Beetstra, van der Hoef and Kuipers2007; Tenneti et al. Reference Tenneti, Garg and Subramaniam2011; Bogner et al. Reference Bogner, Mohanty and Rüde2015; Tang et al. Reference Tang, Peters, Kuipers, Kriebitzsch and van der Hoef2015). In terms of spatial resolution, our PR-DNS is comparable with the work of Bogner et al. (Reference Bogner, Mohanty and Rüde2015) and Tenneti et al. (Reference Tenneti, Garg and Subramaniam2011) for the range of $(\textit {Re}, \phi )$ considered here, but generally better resolved than simulations done by Tang et al. (Reference Tang, Peters, Kuipers, Kriebitzsch and van der Hoef2015) and Beetstra et al. (Reference Beetstra, van der Hoef and Kuipers2007). Nevertheless, validations are presented merely to establish the accuracy of our code rather than intending to offer benchmarking data. The ensemble-average drag is obtained by summing the drag force acting on each particles in the bed and dividing by the total number of particles $N_p$

(2.9)\begin{equation} \langle F_d \rangle = \frac{1}{N_p} \sum_{i=1}^{N_p} \boldsymbol{F}_i \boldsymbol{\cdot} \hat {\boldsymbol{e}}_{x}. \end{equation}

Obviously, the normal components of the mean force $\langle F_y \rangle = \langle \boldsymbol {F} \boldsymbol {\cdot } \hat {\boldsymbol{e}}_{y} \rangle$ and $\langle F_z \rangle = \langle \boldsymbol {F} \boldsymbol {\cdot } \hat {\boldsymbol {e}_z} \rangle$ are both expected to be vanishingly small. In order to compare our data with the available drag correlations, normalisation of forces is done with respect to the Stokes drag given as

(2.10)\begin{equation} F_{st} = 3 \pi \mu \,\textrm{d} U. \end{equation}

As shown in figure 3, the drag computed for our cases listed in table 1 all lie within the range of existing correlations and generally indicate good accordance. In particular, our data seem to agree well with correlations proposed by Bogner et al. (Reference Bogner, Mohanty and Rüde2015) and Tang et al. (Reference Tang, Lau, Deen, Peters and Kuipers2016) at $\textit {Re} = 150$ where significant discrepancies are observed between different correlations. In addition, we have also examined force distributions and the extent of data scatter for each case in table 2. The magnitude of the ensemble-average lift force given as $\langle F_L \rangle = \langle F_y \rangle$ or $\langle F_z \rangle$ is practically zero for all cases, as expected. However, the ratio of the standard deviation of drag and lift shown by $\sigma _{F_d}/\langle F_d \rangle$ and $\sigma _{F_L}/\langle F_d \rangle$ are both significant and $\sigma _{F_d}/\langle F_d \rangle$ turns out to be always greater than $\sigma _{F_L}/\langle F_d \rangle$. In cases where a rough comparison is possible, our results for $\sigma _{F_d}/\langle F_d \rangle$ and $\sigma _{F_L}/\langle F_d \rangle$ in table 2 show a difference of $\approx$2 %–4 % with those obtained by Akiki et al. (Reference Akiki, Jackson and Balachandar2016). A contributing feature that might explain the difference is that our sample sizes for each case are significantly larger by a factor of $5\hbox{--}6$, rendering the results statistically more converged. This is obvious from the $\langle F_L \rangle$ values in table 2 compared with those reported by Akiki et al. (Reference Akiki, Jackson and Balachandar2016), as their mean lift data (while still quite close to zero) are orders of magnitude greater than present results. According to Akiki et al. (Reference Akiki, Jackson and Balachandar2016), realisation dependence alone could cause variance in $\sigma _{F_d}/\langle F_d \rangle$ and $\sigma _{F_L}/\langle F_d \rangle$ of up to $2.7\,\%$ and $1.7\,\%$, respectively.

Figure 3.

Validation of the drag data obtained from PR-DNS of the present work with correlations of Beetstra et al. (Reference Beetstra, van der Hoef and Kuipers2007), Tenneti et al. (Reference Tenneti, Garg and Subramaniam2011), Tang et al. (Reference Tang, Peters, Kuipers, Kriebitzsch and van der Hoef2015) and Bogner et al. (Reference Bogner, Mohanty and Rüde2015).

Table 2.

Statistics of the drag and lift data for cases presented in table 1.

2.5. Dataset construction

For the purpose of analysis and development of our model, we first need to construct datasets corresponding to each simulation. Each dataset contains as many rows as the number of particles $N_p$ (which we call samples), while columns represent the input and output variables. For each particle $i$, the set of position vectors of the neighbouring spheres denoted by $\{ \boldsymbol {r}_{j=1}, \ldots , \boldsymbol {r}_{j=M} \}$ along with the average fluid velocity $\overline {\boldsymbol {u}}_i$ around particle $i$ constitute the inputs, whereas the hydrodynamic forces and torques experienced by particle $i$ are the outputs that we aim to model. The position vectors of neighbours are expressed relative to the location of particle $i$. The first $30$ nearest neighbours are identified for each particle $i$ by looping over all other particles in the simulation. These neighbours are then numbered from $j=1$ to $j=30$ depending on their relative distance, with $j=1$ being the closest. It is worthwhile noting that, since our simulations are performed in tri-periodic domains, forces and torques on particles near boundaries could be affected by periodic images of particles on the opposite side of the domain whose positions do not exist in the simulation data. For this reason, prior to performing the neighbour identification loop, we create the periodic image positions manually and explicitly insert them in the data so that the periodic effects are correctly accounted for.

3.1. Motivation

As pointed out in § 1, we seek to predict the hydrodynamic forces and torques on each particle as a function of the flow conditions, and more importantly, of the unique neighbourhood of each particle as highlighted by (1.2). The hydrodynamic force and torque exerted on each sphere may be decomposed and expressed as follows:

(3.1a)\begin{gather} \boldsymbol{F}_i = \langle \boldsymbol{F}_i \rangle (\textit{Re}, \phi) + {\rm \Delta} \boldsymbol{F}_i(\textit{Re}, \phi, \{ \boldsymbol{r}_{j=1}, \ldots, \boldsymbol{r}_{j=M} \}), \end{gather}

(3.1b)\begin{gather} \boldsymbol{T}_i = {\rm \Delta} \boldsymbol{T}_i(\textit{Re}, \phi, \{ \boldsymbol{r}_{j=1}, \ldots, \boldsymbol{r}_{j=M} \}), \end{gather}

where ${\rm \Delta} \boldsymbol {F}_i$ and ${\rm \Delta} \boldsymbol {T}_i$ show the fluctuating contribution to the force and torque due to the specific arrangement of neighbours surrounding particle $i$. Statistically, the lateral components of $\langle \boldsymbol {F}_i \rangle$ and all components of $\langle \boldsymbol {T}_i \rangle$ are close to zero for a sufficiently large number of particles in each simulation. In the streamwise direction, the ensemble-average force $\langle \boldsymbol {F}_i \rangle$ is identical for all particles by definition and therefore only depends on macroscopic variables, i.e. the Reynolds number and solid volume fraction. Many correlations exist that describe such a dependence as discussed in § 1. For this reason, we henceforth focus on ${\rm \Delta} \boldsymbol {F}_i$ and ${\rm \Delta} \boldsymbol {T}_i$ in (3.1); namely, the deviations from the mean arising due to the unique local neighbourhood of each particle.

The key question posed here may be expressed as follows: How does changing the location of each neighbour affect the force and torque experienced by a particle? The most obvious and naive strategy would be attempting multiple simulations where the location of a single neighbour particle is changed systematically and the resulting effect on forces and torques are recorded. The outcome would be similar to, for example, the two-dimensional interaction maps given by Zhou & Alam (Reference Zhou and Alam2016) for two neighbouring cylinders characterising forces for various locations of a second cylinder. Assuming a three-dimensional system, adding only one or two more neighbours to the scenario makes the number of simulations needed for populating the modelling dataset an almost intractable computational task. For their physics-based PIEP model, Akiki et al. (Reference Akiki, Jackson and Balachandar2017a) needed to account for 15–40 closest neighbours to obtain good accuracy, which in three dimensions translates into 45–120 input variables, respectively. It is well known that in a data-driven framework, such systems are afflicted with the ‘curse of dimensionality’ (Hastie, Tibshirani & Friedman Reference Hastie, Tibshirani and Friedman2009). It means that constructing a reliable dataset for a regression-type model requires millions of samples for each case, while we are only able to model a few thousand particles with PR-DNS at best. The pairwise interactions and order-invariance approximations utilised by Moore et al. (Reference Moore, Balachandar and Akiki2019) for their data-driven model is a good demonstration of a strategy for reducing the number of input variables to a manageable level.

3.2. Theory

Consider a monodisperse random array of spherical particles with a volume fraction $\phi$ subject to fluid flow with a Reynolds number $\textit {Re}$. Let $\boldsymbol {X}$ and $\boldsymbol {Y}$ be two continuous random variable vectors in a $d$-dimensional space ($\boldsymbol {X}, \boldsymbol {Y} \in \boldsymbol {R}^d$) with marginal probability density functions (PDFs) $f_{\boldsymbol {X}}(\boldsymbol {x})$ and $f_{\boldsymbol {Y}}(\boldsymbol {y})$ where $f_{\boldsymbol {X}}(\boldsymbol {x}), f_{\boldsymbol {Y}}(\boldsymbol {Y}) : \boldsymbol {R}^d \rightarrow \boldsymbol {R}$. By definition, PDFs in general are required to satisfy the following conditions:

(3.2)\begin{equation} \left.\begin{array}{c@{}} f_{\boldsymbol{X}}(\boldsymbol{x}) \geq 0,\\ \displaystyle\int f_{\boldsymbol{X}}(\boldsymbol{x}) \,\textrm{d}\boldsymbol{x} = 1, \end{array}\right\} \end{equation}

where the integration is performed on all possible values of $\boldsymbol {x}$. The value $f_{\boldsymbol {X}}(\boldsymbol {x}) \,\textrm {d}\boldsymbol {x}$ represents the probability of $\boldsymbol {X}$ falling in the infinitesimal volume $\textrm {d}\boldsymbol {x}$ about $\boldsymbol {x}$. Therefore, the probability that the random variable $\boldsymbol {X}$ will happen to be in $\mathcal {A}$ is given by

(3.3)\begin{equation} P \left[ \boldsymbol{X} \in \mathcal{A} \right] = \int_{ \mathcal{A} } f_{\boldsymbol{X}}(\boldsymbol{x}) \,\textrm{d}\boldsymbol{x}. \end{equation}

The same properties similarly hold for $f_{\boldsymbol {Y}}(\boldsymbol {y})$. The expected value of the random variable $\boldsymbol {X}$ is the weighted average of all possible values of $\boldsymbol {X}$, each value being weighted according to its probability of occurrence (Ross Reference Ross2010). The expected value of $\boldsymbol {X}$ is thus governed by its probability density distribution $f_{\boldsymbol {X}}(\boldsymbol {x})$, and is given as

(3.4)\begin{equation} E \left[ \boldsymbol{X} \right] = \langle \boldsymbol{X} \rangle = \int \boldsymbol{x} f_{\boldsymbol{X}}(\boldsymbol{x}) \,\textrm{d}\boldsymbol{x}. \end{equation}

The probability associated with a random variable may also be conditioned on another random variable. For this case, $f_{\boldsymbol {X} \mid \boldsymbol {Y}}(\boldsymbol {x} \mid \boldsymbol {y})$ may be defined as the conditional PDF of $\boldsymbol {X}$ given that $\boldsymbol {Y}$ assumes a particular fixed value $\boldsymbol {y}$.

(3.5)\begin{equation} E \left[\boldsymbol{X} \: \vert \: \boldsymbol{Y}=\boldsymbol{y} \right] = \langle \boldsymbol{X} \mid \boldsymbol{Y}=\boldsymbol{y} \rangle = \int \boldsymbol{x} f_{\boldsymbol{x} \mid \boldsymbol{y}} \left(\boldsymbol{x} \: \vert \: \boldsymbol{y} \right)\,\textrm{d}\boldsymbol{x}. \end{equation}

Referring to (3.1), the three components of the force/torque fluctuations pertaining to particle $i$ are shown by ${\rm \Delta} \boldsymbol {F}_i = \left ({\rm \Delta} F_{i,x}, {\rm \Delta} F_{i,y} , {\rm \Delta} F_{i,z}\right ) ^\textrm {T}$ or ${\rm \Delta} \boldsymbol {T}_i = \left ({\rm \Delta} T_{i,x}, {\rm \Delta} T_{i,y} , {\rm \Delta} T_{i,z}\right ) ^\textrm {T}$. A particular configuration of the neighbourhood surrounding particle $i$ is denoted by $\boldsymbol {R}_i = \{ \boldsymbol {r}_{j=1}, \ldots , \boldsymbol {r}_{j=M} \}$. Here, we present the analysis only for force fluctuations; namely, ${\rm \Delta} \boldsymbol {F}$, since the derivation of similar equations is similarly done for torque fluctuations ${\rm \Delta} \boldsymbol {T}$. From this point onwards, the subscript $i$ is dropped from all variables in order to minimise cluttering in the notations. Upon setting $\boldsymbol {X} = {\rm \Delta} F$ and $\boldsymbol {Y} = \boldsymbol {R}$, the following PDF

(3.6)\begin{equation} p_{{\rm \Delta} F \mid \boldsymbol{R}} \left({\rm \Delta} F \: \vert \: \boldsymbol{R} \right) \end{equation}

represents the probability distribution of the force/torque fluctuations given the locations of $M$ neighbours. In the above equations, $\boldsymbol {R}$ shows the set of position vectors $\{\boldsymbol {r}_1, \boldsymbol {r}_2, \boldsymbol {r}_3, \ldots , \boldsymbol {r}_M \}$ belonging to the neighbours numbered according to their distance from a reference particle, starting with $j=1$ being the closest. Each distribution $p \left ({\rm \Delta} F \: \vert \: \boldsymbol {R} \right )$ is to be obtained from collecting ${\rm \Delta} F$ values experienced by particles that happen to have the same neighbour configuration $\boldsymbol {R}$, while no constraints are imposed on neighbours $\{ \boldsymbol {r}_j : j > M \}$ located further away. The expected value of the PDF in (3.6) gives a prediction of each fluctuating force/torque component as follows:

(3.7)\begin{equation} E \left[{\rm \Delta} F \: \vert \: \boldsymbol{R} \right] = \langle {\rm \Delta} F \mid \boldsymbol{R} \rangle = \int \xi p_{{\rm \Delta} F \mid \boldsymbol{R}} \left({\rm \Delta} F \: \vert \: \boldsymbol{R} \right) \,\textrm{d} \xi . \end{equation}

Two extreme cases of $M$, the number of included neighbours, are worth elaborating on.

1. (a) $M=0$: this means that no conditions are imposed on $p$, hence $E[{\rm \Delta} F] = 0$ and the predicted force/torque is the same for all particles in the array. Such a situation corresponds to, for instance, the conventional microstructure-ignorant correlations of the form $\boldsymbol {F} \cdot \hat {\boldsymbol {e}_x} = \langle F_x \rangle (\textit {Re}, \phi )$ in case of the drag force.

2. (b) $M \to N_p-1$: the positions of all neighbours in the particle array are imposed as a condition on $p$, hence the variance of $p$ tends to zero and the value of ${\rm \Delta} F$ is uniquely determined.

Since the problem is fully constrained in the latter case, the predicted value of ${\rm \Delta} F$ would be free of uncertainty. Put another way, $p$ would no longer be a distribution per se, but rather resembles a Dirac delta function instead:

(3.8)\begin{equation} \lim_{M \to N_{p-1}} E \left[{\rm \Delta} F \: \vert \: \boldsymbol{R} \right] = \int_{-\infty}^{+\infty} \xi \delta_{{\rm \Delta} F \mid \boldsymbol{R}} \left( \xi - {\rm \Delta} F_{\textit{DNS}} \: \vert \: \boldsymbol{R} \right) \,\textrm{d} \xi = {\rm \Delta} F_{\textit{DNS}}, \end{equation}

which is equivalent of having a dataset populated with an infinite number of samples. In other words, this ‘ideal’ dataset would contain samples for each and every combination of neighbour locations. Such a dataset would conceivably provide a full and exact description of the problem, thus representing a solution to the Navier–Stokes equations. Consequently, (3.6) can be viewed as a data-driven description of forces/torques exerted on each particle within an array derived from probabilistic arguments, which is potentially capable of predicting fluctuations with varying degrees of accuracy, depending on how constrained we force $p$ to be. It is crucial to realise that unless $M = N_p - 1$, our best estimate of the fluctuating component $E \left [{\rm \Delta} F \: \vert \: \boldsymbol {R} \right ]$ would be an average (i.e. the central tendency) of the distribution, implying that for $M < N_p-1$, influence of neighbours whose positions are not constrained is only accounted for in a statistical sense. Having established the notion and properties of the conditioned force/torque PDFs, the probability-driven prediction of a fluctuating component based on its corresponding distribution is given as

(3.9)\begin{equation} {\rm \Delta} F_{\textit{MPP}} = E \left[{\rm \Delta} F \: \vert \: \boldsymbol{R} \right] = \int_{-\infty}^{+\infty} \xi p_{{\rm \Delta} F \mid \boldsymbol{R}} \left( \xi \: \vert \: \boldsymbol{R} \right) \,\textrm{d} \xi . \end{equation}

The integral above may be split to cover the positive and negative force/torque contributions separately, yielding

(3.10)\begin{equation} {\rm \Delta} F_{\textit{MPP}} = \int_{-\infty}^{0} \xi p_{{\rm \Delta} F \mid \boldsymbol{R}} \left( \xi \: \vert \: \boldsymbol{R} \right) \,\textrm{d} \xi + \int_{0}^{+\infty} \xi p_{{\rm \Delta} F \mid \boldsymbol{R}} \left( \xi \: \vert \: \boldsymbol{R} \right) \,\textrm{d} \xi . \end{equation}

The impeding hurdle with (3.10) is that computing ${\rm \Delta} F_{\textit{MPP}}$ requires a priori knowledge of the force/torque distributions for each and every possible configuration $\boldsymbol {R}$. Such an approach is clearly not feasible and certainly would not fulfil the goal of constructing a model based on a limited dataset. An elegant alternative is to take advantage of Bayes’ theorem for probabilities, which serves to convert a conditional probability problem to its reverse case. According to the Bayes’ formula for distributions of random variables (Papoulis, Pillai & Pillai Reference Papoulis, Pillai and Pillai2002), $f_{\boldsymbol {X} | \boldsymbol {Y}}(\boldsymbol {x} \mid \boldsymbol {y}) = f_{\boldsymbol {Y} | \boldsymbol {X}}(\boldsymbol {y} \mid \boldsymbol {x}) f_{\boldsymbol {X}}(\boldsymbol {x}) / f_{\boldsymbol {Y}}(\boldsymbol {y})$. Therefore, it follows that

(3.11)\begin{equation} p_{{\rm \Delta} F \mid \boldsymbol{R}} \left({\rm \Delta} F \: \vert \: \boldsymbol{R} \right) = \frac{ p_{\boldsymbol{R} \mid {\rm \Delta} F} \left( \boldsymbol{R} \: \vert \: {\rm \Delta} F \right) p_{{\rm \Delta} F}\left({\rm \Delta} F \right) }{ p_{\boldsymbol{R}}\left( \boldsymbol{R} \right) }. \end{equation}

Substituting for $p_{{\rm \Delta} F \mid \boldsymbol {R}} \left ( \xi \: \vert \: \boldsymbol {R} \right )$ in first integral in (3.10) yields

(3.12)\begin{equation} \frac{1}{p_{\boldsymbol{R}}\left( \boldsymbol{R} \right)} \int_{-\infty}^{0} \xi p_{\boldsymbol{R} \mid {\rm \Delta} F} \left( \boldsymbol{R} \: \vert \: \xi \right) p_{{\rm \Delta} F}\left( \xi \right) \textrm{d} \xi . \end{equation}

Assuming the continuity of both $p_{\boldsymbol {R} \mid {\rm \Delta} F}$ and $p_{{\rm \Delta} F}$ and realising that $p_{\boldsymbol {R} \mid {\rm \Delta} F} , p_{{\rm \Delta} F} \geq 0$, the generalised mean value theorem for integrals can be applied to the above equation to give

(3.13)\begin{equation} \frac{\alpha p_{{\rm \Delta} F}(\alpha) }{p_{\boldsymbol{R}}\left( \boldsymbol{R} \right)} \int_{-\infty}^{0} p_{\boldsymbol{R} \mid {\rm \Delta} F} \left( \boldsymbol{R} \: \vert \: \xi \right) \,\textrm{d} \xi = \left[ \frac{\alpha p_{{\rm \Delta} F}(\alpha) }{p_{\boldsymbol{R}}\left( \boldsymbol{R} \right)} \right] p_{\boldsymbol{R} \mid {\rm \Delta} F} \left( \boldsymbol{R} \: \vert \: {\rm \Delta} F < 0 \right), \end{equation}

where $\xi = \alpha$ is a constant value that belongs to the interval $(-\infty , 0]$. Following the same operations for the second integral of (3.10), there exists a value $\xi = \beta$ in $[0, +\infty )$ for which $\xi p_{{\rm \Delta} F}( \xi )$ can be brought out of the integral. Substituting results back in (3.10), we arrive at the following:

(3.14)\begin{equation} {\rm \Delta} F_{\textit{MPP}} = \left[ \frac{\alpha p_{{\rm \Delta} F}(\alpha) }{p_{\boldsymbol{R}}\left( \boldsymbol{R} \right)} \right] p_{\boldsymbol{R} \mid {\rm \Delta} F} \left( \boldsymbol{R} \: \vert \: {\rm \Delta} F < 0 \right) + \left[ \frac{\beta p_{{\rm \Delta} F}(\beta) }{p_{\boldsymbol{R}}\left( \boldsymbol{R} \right)} \right] p_{\boldsymbol{R} \mid {\rm \Delta} F} \left( \boldsymbol{R} \: \vert \: {\rm \Delta} F > 0 \right). \end{equation}

Note that, according to (3.5), $p_{\boldsymbol {R} \mid {\rm \Delta} F} \left ( \boldsymbol {R} \: \vert \: {\rm \Delta} F < 0 \right )$ gives the PDF of particle positions provided that the force/torque fluctuations obey ${\rm \Delta} F < 0$ (same applies to $p_{\boldsymbol {R} \mid {\rm \Delta} F} \left ( \boldsymbol {R} \: \vert \: {\rm \Delta} F > 0 \right )$ when ${\rm \Delta} F > 0$). Since particles are assumed to be randomly, yet uniformly distributed, $p_{\boldsymbol {R}}( \boldsymbol {R} )$ is expected to have an approximately constant distribution. This can be verified by inspection of unconditioned particle location PDFs for various neighbours. Therefore, both terms inside the square brackets in (3.14) are constant values and can be expressed as two coefficients $c_{\alpha }$ and $c_{\beta }$. The final form of the force/torque model reads

(3.15)\begin{equation} {\rm \Delta} F_{\textit{MPP}}(\boldsymbol{R}, \textit{Re}, \phi) = c_{\alpha}\,\tilde{p} \left(\boldsymbol{R} \: \vert \: {\rm \Delta} F < 0 \right) + c_{\beta}\,\tilde{p} \left(\boldsymbol{R} \: \vert \: {\rm \Delta} F > 0 \right), \end{equation}

where $\tilde {p}$ represents the conditional PDF $p_{\boldsymbol {R} \mid {\rm \Delta} F}$. It should be noted that the $\tilde {p}$ probability maps depend on the Reynolds number and solid volume fraction, which is why ${\rm \Delta} F_{\textit{MPP}}$ is shown as ${\rm \Delta} F_{\textit{MPP}}(\boldsymbol {R}, \textit {Re}, \phi )$. Thus far, no approximation was involved in the derivation of ${\rm \Delta} F_{\textit{MPP}}$, which means that in case a full knowledge of the aforementioned distributions exists, (3.15) would be theoretically exact. Practically, however, for a dataset with a few thousand samples the probability distributions $\tilde {p} \left (\boldsymbol {R} \: \vert \: {\rm \Delta} F < 0 \right )$ or $\tilde {p} \left (\boldsymbol {R} \: \vert \: {\rm \Delta} F > 0 \right )$ would be extremely sparse, as alluded to in § 3.1. Taking $N$ and $D$ to denote sample size and number of input dimensions, the data density would be proportional to $N^{1/D}$. As $D$ increases, it is straightforward to see that maintaining a constant sample density requires exponentially more data points. In our case, each neighbour adds three inputs to the problem, each representing one component of the position vector $\boldsymbol {r}_j = (x_j, y_j, z_j)$. If, for instance, $N_1 = 1000$ is deemed a sufficiently dense dataset to form reliable local averages when $D_1 = M \times 3 = 1 \times 3$ (i.e. a single neighbour considered), in case of 15 neighbours $D_2 = 15 \times 3 = 45$, thus $N_2 = N_1^{(D_2/D_1)} = 1000^{15}$. Ordinary PR-DNS-generated simulations of particle arrays are far from being sufficiently dense for creating meaningful, reliable functional forms for $\tilde {p}$. As a result, the curse of dimensionality precludes the construction of PDFs in (3.15) in their current form. Any such effort is in fact bound to produce a severely over-fitted model with almost no generalisable prediction capability. Consequently, instead of attempting to obtain $\tilde {p} \left (\boldsymbol {R} \: \vert \: {\rm \Delta} F < 0 \right )$ and $\tilde {p} \left (\boldsymbol {R} \: \vert \: {\rm \Delta} F > 0 \right )$ directly, we suggest that

(3.16)\begin{equation} \left.\begin{array}{c@{}} \tilde{p} \left(\boldsymbol{R} \: \vert \: {\rm \Delta} F < 0 \right) \approx c_{\alpha, 1}\,\tilde{p}_1 \left(\boldsymbol{r}_1 \: \vert \: {\rm \Delta} F < 0 \right) + c_{\alpha, 2}\,\tilde{p}_2 \left(\boldsymbol{r}_2 \: \vert \: {\rm \Delta} F < 0 \right) + \cdots \\ = \displaystyle \sum_{j=1}^M c_{\alpha, j} \, \tilde{p}_j \left(\boldsymbol{r}_j \: \vert \: {\rm \Delta} F < 0 \right),\\ \tilde{p} \left(\boldsymbol{R} \: \vert \: {\rm \Delta} F > 0 \right) \approx c_{\beta, 1}\,\tilde{p}_1 \left(\boldsymbol{r}_1 \: \vert \: {\rm \Delta} F > 0 \right) + c_{\beta, 2}\,\tilde{p}_2 \left(\boldsymbol{r}_2 \: \vert \: {\rm \Delta} F > 0 \right) + \cdots\\ = \displaystyle \sum_{j=1}^M c_{\beta, j} \, \tilde{p}_j \left(\boldsymbol{r}_j \: \vert \: {\rm \Delta} F > 0 \right). \end{array}\right\} \end{equation}

In the above equation, each function $\tilde {p}_j \left (\boldsymbol {r}_j \: \vert \: \cdot \right )$ is equivalent to a marginal probability distribution which is obtained by integrating out the position $\boldsymbol {r}_k$ of each and every neighbour where $k \ne j$ in $\tilde {p} \left (\boldsymbol {R} = \{ \boldsymbol {r}_{k=1}, \ldots , \boldsymbol {r}_{k=M} \} \: \vert \: \cdot \right )$. In other words, this approximation only considers the effect of a single neighbour on the functional form of $\tilde {p}$, while accounting for the presence of all other surrounding particles in an average manner. This marginal distribution can be expressed as

(3.17)\begin{equation} \tilde{p}_j \left(\boldsymbol{r}_j \: \vert \: \cdot \right) = \int \tilde{p} \left(\boldsymbol{R} \: \vert \: \cdot \right) \,\textrm{d}\boldsymbol{r}_1 \ldots \textrm{d}\boldsymbol{r}_{j-1}\,\textrm{d}\boldsymbol{r}_{j+1} \ldots \textrm{d}\boldsymbol{r}_M. \end{equation}

Equation (3.16) is reminiscent of the pairwise approximation employed by Moore et al. (Reference Moore, Balachandar and Akiki2019) for the purpose of reducing the number of independent variables by accounting for the influence of only one neighbour at a time. In their data-driven model, the authors also invoked the ‘order-invariance’ approximation, which removes the dependence of the model functions on the neighbour number. This was done by weighting the functionals according to the probability of a particular neighbour being the $j$th closest neighbour. In the MPP model, the neighbours first need to be sorted based on their distance to the reference particle. Each neighbour's position can be then passed to the appropriate PDF.

Substitution of (3.16) in (3.15) leads to

(3.18)\begin{equation} {\rm \Delta} F_{\textit{MPP}}(\boldsymbol{R}, \textit{Re}, \phi) = \sum_{j=1}^M c_{\alpha, j} \, \tilde{p}_j \left(\boldsymbol{r}_j \: \vert \: {\rm \Delta} F < 0 \right) + \sum_{j=1}^M c_{\beta, j} \, \tilde{p}_j \left(\boldsymbol{r}_j \: \vert \: {\rm \Delta} F > 0 \right). \end{equation}

In (3.18), $c_{\alpha , j}$ and $c_{\beta , j}$ are constant unknown coefficients. With the final form of ${\rm \Delta} F_{\textit{MPP}}$, we have assumed an additive nature for the effect that each neighbour might have on the force/torque fluctuation experienced by a reference particle. Notably, both distributions $\tilde {p}_j \left (\boldsymbol {r}_j \: \vert \: \cdot \right )$ (this notation is used to refer to a PDF with any arbitrary force/torque conditioning) are now functions of three variables only; that is, the three components of the position vector $\boldsymbol {r}_j = (r_{j,x}, r_{j,y}, r_{j,z})$. This means that the functional form of the PDFs in (3.18) can be conveniently inferred with a data-driven approach, as the high-dimensionality of the input space has been evaded. To this end, a dataset can be constructed by running PR-DNS of fixed beds for desired sets of parameters. The discrete estimation of the distribution $\tilde {p}_j$ can be extracted by filtering neighbour positions according to the conditions $\{ {\rm \Delta} F < 0, {\rm \Delta} F > 0\}$. The discrete form of the distribution is then fitted with a multivariate kernel density estimation (KDE) function (Chacón & Duong Reference Chacón and Duong2018) of the form

(3.19)\begin{equation} \tilde{p} \left(\boldsymbol{r} \: \vert \: \cdot \right) = \frac{1}{m} \sum_{q=1}^{m} K_{\boldsymbol{\mathsf{H}}}(\boldsymbol{r} - \boldsymbol{r}_q), \end{equation}

where

(3.20)\begin{equation} K_{\boldsymbol{\mathsf{H}}}(\boldsymbol{r}) = {|\boldsymbol{\mathsf{H}}|}^{-1/2} K(\boldsymbol{\mathsf{H}}^{-1/2} \boldsymbol{r}). \end{equation}

In the above equations, $m$ shows the total number of samples after filtering with the aforementioned conditions, $\boldsymbol{\mathsf{H}}$ the bandwidth matrix and $K$ is a symmetric density function of choice. We select a Gaussian distribution as the kernel function defined as

(3.21)\begin{equation} K(\boldsymbol{r}) = \frac{1}{(2 \pi)^{D/2}} \exp \left\{ -\frac{1}{2} \boldsymbol{r}^{\textrm{T}} \boldsymbol{r} \right\}. \end{equation}

We will demonstrate and discuss particular examples of these functions in § 3.3. Now that the appropriate functional forms are identified and known, it remains to decide the values of the constant coefficients $c_{\alpha , j}$ and $c_{\beta , j}$ in (3.18). The values are obtained using the ordinary least-squares method for linear regression

(3.22)\begin{equation} \{ c_{\alpha, j}, c_{\beta, j} \}_{j=1}^{M} = \text{argmin} \sum_{i=1}^{N_p} \left( {\rm \Delta} F_{\textit{MPP}, i} - {\rm \Delta} F_{\textit{DNS}, i} \right)^{2}, \end{equation}

which minimises the residual sum of squares between values of force/torque fluctuations predicted by the MPP model and the true values from PR-DNS. After optimal values for $c_{\alpha , j}$ and $c_{\beta , j}$ are found, (3.18) can be used to make predictions for each particle based on the particular microstructure of the surrounding particles. It is of great importance to realise that in the presented approach, the model relies solely on the unique configuration of neighbours surrounding each particle, and the PR-DNS-generated force/torque data (as opposed to velocity and pressure fields) in order to make predictions. Consequently, the MPP model accounts for the combined effect of various force contributions such as the undisturbed flow, quasi-steady and added-mass forces. This can be contrasted with the PIEP physics-driven model (Akiki et al. Reference Akiki, Jackson and Balachandar2017a) which uses perturbed velocity and pressure fields together with the Faxén's theorem to obtain hydrodynamic forces experienced by each particle.

The above analysis applies identically to all force and torque fluctuation components; namely, $\left ({\rm \Delta} F_{x}, {\rm \Delta} F_{y} , {\rm \Delta} F_{z}\right )$ and $\left ({\rm \Delta} T_{x}, {\rm \Delta} T_{y} , {\rm \Delta} T_{z}\right )$. Let $\boldsymbol {C}_{\alpha }$ denote a $1 \times M$ vector that contains the $c_{\alpha , j}$ coefficients in (3.18)

(3.23)\begin{equation} \boldsymbol{C}_{\alpha} =\begin{pmatrix} c_{\alpha, 1} & \cdots & c_{\alpha, j} & \cdots & c_{\alpha, M} \end{pmatrix} \end{equation}

and let $\boldsymbol {P}_{{\rm \Delta} F < 0}$ denote an $M \times 1$ vector that contains the PDFs $\tilde {p}_j \left (\boldsymbol {r}_j \: \vert \: {\rm \Delta} F < 0 \right )$ in (3.18)

(3.24)\begin{equation} \boldsymbol{P}_{{\rm \Delta} F < 0} =\begin{pmatrix} \tilde{p}_1 \left(\boldsymbol{r}_1 \: \vert \: {\rm \Delta} F < 0 \right) \\ \vdots \\ \tilde{p}_j \left(\boldsymbol{r}_j \: \vert \: {\rm \Delta} F < 0 \right) \\ \vdots \\ \,\tilde{p}_M \left(\boldsymbol{r}_M \: \vert \: {\rm \Delta} F < 0 \right) \end{pmatrix}. \end{equation}

Together with their counter-parts, i.e. $\boldsymbol {C}_{\beta }$ and $\boldsymbol {P}_{{\rm \Delta} F > 0}$, the above vectors can be ultimately defined for all force components. The complete model equations for force fluctuations can thus be expressed in the following form:

(3.25)\begin{equation} {\rm \Delta} \boldsymbol{F}_{\textit{MPP}}(\boldsymbol{R}, \textit{Re}, \phi) =\begin{pmatrix} {\rm \Delta} F_x \\ {\rm \Delta} F_y \\ {\rm \Delta} F_z \end{pmatrix}=\begin{pmatrix} \boldsymbol{C}^x_{\alpha} \boldsymbol{P}_{{\rm \Delta} F_x < 0} + \boldsymbol{C}^x_{\beta} \boldsymbol{P}_{{\rm \Delta} F_x > 0} \\ \boldsymbol{C}^y_{\alpha} \boldsymbol{P}_{{\rm \Delta} F_y < 0} + \boldsymbol{C}^y_{\beta} \boldsymbol{P}_{{\rm \Delta} F_y > 0} \\ \boldsymbol{C}^z_{\alpha} \boldsymbol{P}_{{\rm \Delta} F_z < 0} + \boldsymbol{C}^z_{\beta} \boldsymbol{P}_{{\rm \Delta} F_z > 0} \end{pmatrix}, \end{equation}

where the superscripts of the coefficients indicate the corresponding component. Following the above definitions, torque fluctuations are given as

(3.26)\begin{equation} {\rm \Delta} \boldsymbol{T}_{\textit{MPP}}(\boldsymbol{R}, \textit{Re}, \phi) =\begin{pmatrix} {\rm \Delta} T_x \\ {\rm \Delta} T_y \\ {\rm \Delta} T_z \end{pmatrix}=\begin{pmatrix} \boldsymbol{D}^x_{\alpha} \boldsymbol{P}_{{\rm \Delta} T_x < 0} + \boldsymbol{D}^x_{\beta} \boldsymbol{P}_{{\rm \Delta} T_x > 0} \\ \boldsymbol{D}^y_{\alpha} \boldsymbol{P}_{{\rm \Delta} T_y < 0} + \boldsymbol{D}^y_{\beta} \boldsymbol{P}_{{\rm \Delta} T_y > 0} \\ \boldsymbol{D}^z_{\alpha} \boldsymbol{P}_{{\rm \Delta} T_z < 0} + \boldsymbol{D}^z_{\beta} \boldsymbol{P}_{{\rm \Delta} T_z > 0} \end{pmatrix}, \end{equation}

where $\boldsymbol {D}_{\alpha }$, $\boldsymbol {D}_{\beta }$, $\boldsymbol {P}_{{\rm \Delta} T > 0}$ and $\boldsymbol {P}_{{\rm \Delta} T < 0}$ for each torque component are equivalent to $\boldsymbol {C}_{\alpha }$, $\boldsymbol {C}_{\beta }$, $\boldsymbol {P}_{{\rm \Delta} F > 0}$ and $\boldsymbol {P}_{{\rm \Delta} F < 0}$ in (3.25). Note that the vectors $\boldsymbol {C}_\alpha$, $\boldsymbol {C}_\beta$, $\boldsymbol {D}_\alpha$, $\boldsymbol {D}_\beta$, $\boldsymbol {P}_{{\rm \Delta} F < 0}$, $\boldsymbol {P}_{{\rm \Delta} F > 0}$, $\boldsymbol {P}_{{\rm \Delta} T < 0}$, $\boldsymbol {P}_{{\rm \Delta} T > 0}$ each contain as many coefficients and PDFs as the number of included neighbours $M$, the values of which depend on $\textit {Re}$ and $\phi$ of each case represented in table 1.

3.3. Probability distribution maps

In this section, we scrutinise a few of the probability distribution maps obtained for different force/torque components and various neighbours. Before proceeding further, we first describe the necessary steps taken for generating the foregoing distributions. In order to construct the maps, those particles that satisfy the desired condition (e.g. ${\rm \Delta} F_x < 0$ or $|{\rm \Delta} F_y| > \sigma$, with $\sigma$ being the standard deviation of the corresponding variable) are identified and filtered by looping over the entire array of spheres. In this filtered subset of the original data, the positions of the chosen neighbours are recorded. These positions constitute data points to which KDE functions in (3.19)–(3.21) are then fitted, thus giving $\tilde {p}_j \left (\boldsymbol {r}_j \: \vert \: \cdot \right )$. The generated PDFs are functions of the three spatial coordinates, which makes it impossible to visualise these functions directly. Inevitably, we resort to contour plots of the PDFs; however, three-dimensional contour surfaces are also neither efficient nor convenient to visualise. Consequently, we have opted for depicting their projections on two-dimensional planes. In doing so, the contour surfaces are entirely compressed onto an $x$–$y$ plane. This is why the reference particle will be shown as a blank circle in contrast to a solid disk, since the neighbouring spheres can happen to be located even at $(x,y) = (0,0)$ when $|z| \geq 0.5$.

In the case of streamwise force, all maps are in fact axisymmetric about the $x$ axis; meaning that the choice of the plane on which the maps are shown does not make any difference, as long as it is parallel to the $x$ axis. As for the lateral forces, the maps are plane symmetric with respect to the same plane whose normal has the same direction as the force component (e.g. $\tilde {p} \left (\boldsymbol {r}_j \: \vert \: |{\rm \Delta} F_y| > \sigma \right )$ would be symmetric about the $y$ axis). For the hydrodynamic torques, the maps are also plane symmetric but with respect to the plane that contains the torque component. Owing to the symmetry about the streamwise direction, the probability distribution map for $|{\rm \Delta} F_y | > \sigma$ shown on the $x$–$y$ plane is essentially almost identical to that for $|{\rm \Delta} F_z | > \sigma$ on the $y$–$z$ plane. The same is also true for $|{\rm \Delta} T_z | > \sigma$ on $x$–$y$ and $|{\rm \Delta} T_y | > \sigma$ on $y$–$z$. Noting the existence of this similarity, only projections on the $x$–$y$ plane will be shown in this section. With an ideal densely populated dataset, the PDFs are expected to exhibit perfect symmetry. The limited number of samples in our dataset, however, causes some PDFs to show minor yet noticeable deviations from symmetry. This would be particularly evident in PDFs generated for farther neighbours due to the lower data density compared with closer neighbours. In what follows, statistical symmetry is hence enforced in all PDF depictions.

Examination of the probability maps in all cases remarkably reveals quite distinct, non-uniform and physically meaningful functional forms. Figure 4 shows the distribution of the first (i.e. closest) neighbour positions, conditioned on two different ranges of the streamwise force, while the blank circle in the middle of represents the reference particle which experiences the hydrodynamic force or torque. According to these probability maps, it is considerably more likely for a reference particle to experience higher than the mean drag force when the closest neighbour is located laterally, and slightly downstream of the particle. The right plot in figure 4 appears more intuitive: drag force diminishes if the reference particle is shielded by the closest neighbour directly in front of it, as expected. A somewhat less intuitive observation is that even if the neighbour happens to be immediately downstream of the reference particle, drag force would be lower than average, due to the suppression of the low-pressure region behind the reference particle. It is crucial to stress this point once more that these maps demonstrate the higher likelihood of coming across neighbours at particular regions depending on the condition imposed on force/torque fluctuation. This means that, for instance, even if the closest neighbour occurs to be located laterally, the drag force might be lower than average because the second neighbour happens to directly shield the reference particle. It signifies the fact that while accounting for the influence of one neighbour provides a great deal of valuable information about the fluctuations, effects of other neighbours are indispensable in explaining the force/torque fluctuation accurately.

Figure 4.

(a) Unconditioned PDF of the first closest neighbour position and PDF of the first closest neighbour position when the reference particle experiences a (b) higher than average or (c) lower than average drag force. The PDFs are obtained for the case of $Re = 40$ and $\phi = 0.1$.

Figures 5 and 6 show PDFs for the distribution of closest neighbours when conditions are imposed on the lateral force in the $y$ direction and on the lateral torque in the $z$ direction, respectively. Note that for demonstration purposes, we have chosen more extreme conditions in figures 5 and 6, e.g. ${\rm \Delta} F_y > \sigma$ and ${\rm \Delta} F_y < -\sigma$, and ${\rm \Delta} T_z > \sigma$ and ${\rm \Delta} T_z < - \sigma$, respectively, in order to make the PDFs stand out more prominently. However, the MPP model ((3.25) and (3.26)) is always constructed with the original positive and negative conditioning (e.g. ${\rm \Delta} F_y > 0$ and ${\rm \Delta} F_y < 0$ for the lateral force in the $y$ direction). On the leftmost plot of figure 5, the reference particle is seen to experience a positive lateral force in the $y$ direction when the first neighbour is located within the region that lies mostly in the bottom right quadrant where $x > 0$ and $y < 0$, whereas a negative lateral force results when the first neighbour happens to be on the opposite side where $x, y > 0$ depicted in the middle plot. As pointed out by Prahl et al. (Reference Prahl, Hölzer, Arlov, Revstedt, Sommerfeld and Fuchs2007) for a staggered arrangement of only two spheres, a higher than normal pressure field develops in the gap between the two particles, which causes a repulsive force acting mainly on the upstream particle. Upon increasing the distance between the particles, the fluid can penetrate the gap more easily. The nozzle effect hence creates a low-pressure region in the gap which gives rise to an attraction force. The latter does not appear in figure 5 as the first neighbour is quite close to the reference particle, but regions corresponding to such a condition show up remarkably on the PDFs for the lateral force in case we consider more distant neighbours. The right most plot shows the PDFs conditioned with $|{\rm \Delta} F| > \sigma$, which is a combination of the foregoing plots with $\{ {\rm \Delta} F < -\sigma , {\rm \Delta} F > \sigma \}$. These are regions for which the reference particle experiences a non-zero lateral force in the $y$ direction.

Figure 5.

The PDF of the first closest neighbour position when the reference particle experiences a (a) higher than $\sigma$ or (b) lower than $\sigma$ lift, with $\sigma$ being the standard deviation of the data. (c) The PDF of the first closest neighbour position when the lift is either higher or lower than $\sigma$. The PDFs are obtained for the case of $Re = 40$ and $\phi = 0.1$.

Figure 6.

The PDF of the first closest neighbour position when the reference particle experiences a (a) higher than $\sigma$ and (b) lower than $\sigma$ lateral torque, with $\sigma$ being the standard deviation of the data. (c) The PDF of the first closest neighbour position when the lateral torque is either higher or lower than $\sigma$. The PDFs are obtained for the case of $Re = 40$ and $\phi = 0.1$.

In figure 6, we can see essentially the same kind of information as in figure 5. The most critical regions for the modification of the torque are located where the highest amount of vorticity is generated in the flow over a single sphere. As with the lateral force, the higher pressure in the gap decreases the fluid velocity and the resulting vorticity. The vorticity imbalance between the top and bottom regions of the reference particle brings about a net torque exerted on the particle. The interpretation of the PDFs emerging from our analysis is consistent with the observations made with a binary system of two spheres exhibiting the same dependency of the forces on relative positions of the spheres (Prahl et al. Reference Prahl, Hölzer, Arlov, Revstedt, Sommerfeld and Fuchs2007; Yoon & Yang Reference Yoon and Yang2007). Such modification of the forces and torques were also elucidated by Akiki et al. (Reference Akiki, Jackson and Balachandar2016). Our results imply that even when the system being modelled is a dense multi-particle system, the forces and torques acting on each particle still vary in the same manner on average, in response to the positioning of the closest neighbour. For instance, when a particle happens to experience less drag force, it is most probably shielded by an upstream neighbour (or supported by a neighbour immediately downstream with a comparatively less probability) as shown on the right plot in figure 4. Although this may not be always true, it is the most likely reason. Similar arguments apply in regards to all other conditional PDFs involving lateral forces and torques as well.

Interestingly, similar patterns are found for neighbours located further away. An example of the PDFs for farther neighbours is shown in figure 7. While the distributions are spread out and have become fainter for surrounding particles located farther away from the reference particle, regions with higher probability are discernible and turn out to be generally in accordance with those obtained for $j=1$, i.e. the closest neighbour. The reason for the fact that the densities are smaller for $j > 1$ is twofold. First, farther particles are distributed over a larger volume around the reference particle. The second and more important reason is that the more distant particles are comparatively less likely to alter the force/torque fluctuations. The functional forms of $\tilde {p}$ in (3.16) has been kept distinct for each neighbour, meaning that $\tilde {p}_1$ is a different function compared with $\tilde {p}_2$, and they both are different from $\tilde {p_3}$ and so on. This is also reflected in the probability maps shown in figure 7. A possible simplification can be made by combining all functions $\tilde {p}_j$ into a single function $\tilde {p}$ that still contains the contribution of all included neighbours, except for not making distinction based on the ordering of the neighbours. In this manner, instead of filtering the locations of the $j$th neighbour, we filter the location of all the included neighbours (up to $j = M$) based on the desired condition. It is to be noted that we still retain the contribution of each neighbour to the net force/torque fluctuation by allowing separate terms for every included neighbour in (3.16). That is to say, we construct only one PDF, namely, $\tilde {p}$ using the entire neighbourhood, while each term $\tilde {p} \left (\boldsymbol {r}_j \: \vert \: \cdot \right )$ is assigned a separate unknown coefficient and influences the prediction of fluctuations differently. An example of such a combined PDF is shown for the streamwise force in figure 8(a) for $\phi = 0.1$. Comparing with figure 4, it is immediately obvious that the significant region around the reference particle has expanded due to the influence of more distant neighbours. Nevertheless, the general form of the maps that were obtained for $j = 1$ in figure 4 still prevails even when all the neighbours up to $j = M$ are included all at once. Another remarkable observation can be made from figure 8(b), where PDFs similar to those in figure 8(a) have been shown for the highest solid volume fraction of $\phi = 0.4$. The PDFs in figure 8(b) have gained fore-aft symmetry and are shrunk towards the reference particle, reminiscent of the force maps obtained with the data-driven model of Moore et al. (Reference Moore, Balachandar and Akiki2019) at $\phi = 0.45$.

Figure 7.

The PDFs of various neighbour positions for (a) ${\rm \Delta} F_x < - \sigma$, and (b) ${\rm \Delta} F_x > \sigma$, with $\sigma$ being the standard deviation of the data. Note that the numbering represented by $j$ is based on proximity to the reference particle, where $j=1$ shows the closest neighbour. The PDFs are obtained for $Re = 40$ and $\phi = 0.1$.

Figure 8.

PDFs of positions of $15$ closest neighbours for two different solid volume fractions. (a) $\textit {Re} = 40$ and $\phi = 0.1$, (b) $\textit {Re} = 40$ and $\phi = 0.4$.