2.1 Definition and pairwise approximation of undisturbed flow

A schematic of a typical Eulerian–Lagrangian (EL) simulation is shown in figure 1. Also shown in the figure is a zoom-in of a single grid cell consisting of many particles inside. If the EL simulation employs a first-order finite volume approach, then the fluid velocity seen by all the particles shown in figure 1(b) is typically taken to be the same, and equal to the cell-centre fluid velocity. Furthermore, if the particles are stationary then the force on all the particles predicted with the standard point-particle models is the same, and there will not be any particle-to-particle force variation predicted within the cell. If the EL simulation employs a finite difference approach or a finite volume approach with second-order interpolation, then the fluid velocity at the particle location will be interpolated from the surrounding grid points or cell-centre values of fluid velocity. In this case, the fluid velocity seen by the different particles within a grid cell can vary and correspondingly there can be particle-to-particle force variation within a grid cell. However, this fluid velocity variation seen by the particles is due to continuous phase velocity variation at the sub-grid scale, which is typically small. The interpolated fluid velocity at a particle, however, does not account for the perturbation fluid velocity due to the influence of the neighbouring particles within the cell. Since the location of all the particles within the cell and their sizes are precisely known, we wish to address the following question: can the undisturbed flow seen by the $i$ th particle (marked blue in figure 1 b) be more accurately evaluated from the knowledge of the cell-centre (or grid point) value of the flow velocity at the macroscale and the location of all the neighbouring particles?

Figure 1.

General schematic for EL simulations at the (a) macroscale and (b) microscale.

Consider a random distribution of spherical particles subjected to a macroscale velocity field, $\boldsymbol{u}_{mac}(x)$ , where $x$ is the spatial coordinate vector, as computed in an EL simulation. Here the macroscale velocity field has been computed at the cell centre or grid points and can be interpolated to obtain the velocity at any point $x$ within the computational domain. In a two-way coupled EL simulation, where the dispersed-to-continuous back coupling is accounted for, this computed macroscale velocity field $\boldsymbol{u}_{mac}(x)$ is, in a sense, influenced by the collective effect of all the particles.

We define $\boldsymbol{u}_{mic}(x_{i})$ as the microscale flow at the location $x_{i}$ occupied by the $i$ th particle due to the perturbation flow created by the presence of all other neighbouring particles. The undisturbed flow at the location of the $i$ th particle can then be expressed as a superposition of the macroscale flow and the microscale perturbation as,

(2.1) $$\begin{eqnarray}\displaystyle \boldsymbol{u}_{un,i}(x_{i})=\boldsymbol{u}_{mac}(x_{i})+\boldsymbol{u}_{mic,i}(x_{i}). & & \displaystyle\end{eqnarray}$$

The above separation of the undisturbed flow into a macroscale and a microscale component can be well defined in the limit where the particles are much smaller than the grid. This is the limit typically considered in the traditional point-particle approach. The Eulerian fluid velocity evaluated at the $i$ th particle location is the macroscale velocity, $\boldsymbol{u}_{mac}$ , which will be taken to be the undisturbed fluid velocity as seen by the $i$ th particle in the traditional point-particle approach. This macroscale velocity does not account for the sub-grid scale perturbation velocity induced by all other particles at the microscale. In (2.1) $\boldsymbol{u}_{mic}(x_{i})$ accounts for this microscale component of the undisturbed flow at the $i$ th particle. It must be pointed out that $\boldsymbol{u}_{mic}(x_{i})$ , thus, accounts for the pseudo-turbulence generated by the individual particles, while $\boldsymbol{u}_{mac}(x)$ can be a laminar flow or a turbulent flow, with flow scales only above the grid size represented.

The separation of the undisturbed flow becomes complicated, and perhaps even non-unique, when the size of the particle approaches or exceeds the grid spacing. Since the continuous phase velocity computed at the grid points is substantially influenced by the presence of such large particles, obtaining the undisturbed flow velocity becomes a challenge, even in the implementation of the point-particle approach. Recent advances by Gualtieri et al. (Reference Gualtieri, Picano, Sardina and Casciola2015), Horwitz & Mani (Reference Horwitz and Mani2016) have presented physical and numerical methodologies for the removal of self-disturbance so that the undisturbed flow velocity can be extracted. Such difficulties will persist with the application of PIEP model. Therefore, here we limit our attention to scenarios where the particles are much smaller than the grid as depicted in figure 1.

In (2.1) and henceforth we do not limit $x_{i}$ to be only the centre of the $i$ th particle, but to extend over the entire volume of the $i$ th particle. Thus, the undisturbed flow at the location of the $i$ th particle implies velocity distribution over the entire volume occupied by the particle. Even when the macroscale flow is spatially uniform, the perturbation flow due the neighbouring particles will be strongly inhomogeneous, and as a result $\boldsymbol{u}_{un,i}(x_{i})$ can have substantial variation over the volume occupied by the particle. This generalized definition of a spatially varying undisturbed flow is important, since later in the Faxén form of the quasi-steady and other forces, we require the undisturbed flow to be averaged over the surface or volume of the particle.

Our goal is then to find a suitable expression for the microscale component of the undisturbed velocity at the $i$ th particle by accounting for the relative position of all the neighbouring particles. In the limit of very dilute distribution of particles, the microscale correction $\boldsymbol{u}_{mic,i}(x_{i})\rightarrow 0$ and the undisturbed flow will simply approach the macroscale flow. With increasing particle volume fraction the mean inter-particle distance goes as $\unicode[STIX]{x1D719}^{-1/3}$ , where $\unicode[STIX]{x1D719}$ is the local volume fraction. So, even at modest particle volume fraction, the likelihood of finding neighbouring particles sufficiently close to contribute to the microscale component of the undisturbed flow increases.

A precise evaluation of the undisturbed flow at the location of the $i$ th particle, $\boldsymbol{u}_{un,i}(x_{i})$ , can be obtained from a fully resolved DNS simulation that includes all the other particles, except the $i$ th particle. But this formal definition is useful only as a conceptual picture, since the undisturbed flow at the $N$ particles will require $N$ fully resolved DNS simulations, each involving $N-1$ particles. Our interest is to obtain a good approximation for the undisturbed flow at the location of each of the particles, by taking into account the leading-order effect of the neighbours.

Here we propose a pairwise interaction approximation, where the perturbation flow of each of the neighbours is considered separately. The total effect of all the neighbours would then be the superposition of the individual effects. Using this pairwise interaction approximation, the undisturbed flow at the location of the $i$ th particle can then be simplified as,

(2.2) $$\begin{eqnarray}\displaystyle \boldsymbol{u}_{un,i}(x_{i})\approx \boldsymbol{u}_{mac}(x_{i})+\mathop{\sum }_{\substack{ j=1 \\ j\neq i}}^{N}\boldsymbol{u}_{j\rightarrow i}(x_{i}), & & \displaystyle\end{eqnarray}$$

where $\boldsymbol{u}_{j\rightarrow i}$ is the velocity perturbation at the $i$ th particle caused by the presence of the $j$ th particle. Here the sum is over all the $N$ particles except $j=i$ . In (2.2) the second term on the right is only an approximation to $\boldsymbol{u}_{mic,i}(x_{i})$ , since it accounts for the effect of the neighbouring particles one at a time and thereby ignores tertiary, $\boldsymbol{u}_{k\rightarrow j\rightarrow i}$ , and higher-order fluid-mediated interactions between the particles, which are implicitly captured in a fully resolved DNS.

An important advantage of the pairwise interaction approximation is that it can be implemented in an EL simulation. According to (2.2) the velocity perturbation at the $i$ th particle caused by the presence of the $j$ th particle is independent of the presence of all other particles in the neighbourhood. It depends only on the macroscale flow $\boldsymbol{u}_{mac}(x)$ and the relative location of the $i$ th particle with respect to the $j$ th particle. For a given macroscale flow, if we place the $j$ th particle at its location, it will create a perturbation flow due to its no-slip and no-penetration boundary conditions. This perturbation flow evaluated at the location of the $i$ th particle is $\boldsymbol{u}_{j\rightarrow i}(x_{i})$ . Since the definition of this perturbation flow excludes the presence of all particles except the $j$ th particle, this can be obtained from the flow over an isolated particle. Further, if we consider the macroscale flow to be uniform at the level of a single grid cell (as in a first-order finite volume approach), $\boldsymbol{u}_{j\rightarrow i}(x_{i})$ can be readily obtained from the finite Reynolds number uniform flow around an isolated particle. This evaluation of the perturbation flow from DNS of an isolated particle will be presented in § 3.

The undisturbed flow at the location of the $i$ th particle will be used to calculate the quasi-steady, pressure-gradient, added-mass, Basset history and vorticity-induced forces on the $i$ th particle. The evaluation of these forces requires not only the undisturbed velocity field at the location of the $i$ th particle, but also other related quantities such as the stress tensor and the vorticity of the undisturbed flow. We separate these quantities as well into a macroscale and a microscale contribution. We again employ the pairwise superposition approximation for the microscale contribution. The stress tensor of the undisturbed flow at the location of the $i$ th particle can be written as,

(2.3) $$\begin{eqnarray}\displaystyle -p_{un,i}(x_{i})\boldsymbol{I}+\unicode[STIX]{x1D749}_{un,i}(x_{i}) & {\approx} & \displaystyle -(p_{mac}(x_{i})+p_{mic}(x_{i}))\boldsymbol{I}+\unicode[STIX]{x1D749}_{mac}(x_{i})+\unicode[STIX]{x1D749}_{mic}(x_{i})\nonumber\\ \displaystyle & {\approx} & \displaystyle -\!\left(p_{mac}(x_{i})+\mathop{\sum }_{\substack{ j=1 \\ j\neq i}}^{N}p_{j\rightarrow i}(x_{i})\right)\boldsymbol{I}+\unicode[STIX]{x1D749}_{mac}(x_{i})+\mathop{\sum }_{\substack{ j=1 \\ j\neq i}}^{N}\unicode[STIX]{x1D749}_{j\rightarrow i}(x_{i}),\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

where $p_{mac}$ and $\unicode[STIX]{x1D749}_{mac}$ are the pressure and viscous stress tensor of the macroscale flow and $p_{j\rightarrow i}(x_{i})$ and $\unicode[STIX]{x1D749}_{j\rightarrow i}(x_{i})$ are the perturbation pressure and perturbation viscous stress tensor due to the presence of the $j$ th particle at the location of the $i$ th particle, ignoring the effect of all other neighbours.

Note that the above definition involves different layers of approximation. Due to the nonlinear dependence of pressure on velocity, if we define macroscale pressure to obey the following Poisson equation: $\unicode[STIX]{x1D6FB}^{2}p_{mac}=-\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\boldsymbol{u}_{mac}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}_{mac})$ , the corresponding Poisson equation for microscale pressure becomes

(2.4) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FB}^{2}p_{mic}=-\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\boldsymbol{u}_{mac}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}_{mic})-\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\boldsymbol{u}_{mic}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}_{mac})-\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\boldsymbol{u}_{mic}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}_{mic}) & & \displaystyle\end{eqnarray}$$

and thus $p_{mic}$ will depend on both micro and macroscale velocity. Here we assume $\boldsymbol{u}_{mac}$ to be nearly uniform on the microscale, and ignore the first two contributions on the right-hand side. In obtaining the final expression in (2.3) we assume pairwise interaction with additive contribution from each neighbour. Here again we ignore the nonlinear interaction between the different perturbations with the assumption

(2.5) $$\begin{eqnarray}\displaystyle \left(\mathop{\sum }_{\substack{ j=1 \\ j\neq i}}^{N}\boldsymbol{u}_{j\rightarrow i}\right)\boldsymbol{\cdot }\unicode[STIX]{x1D735}\left(\mathop{\sum }_{\substack{ j=1 \\ j\neq i}}^{N}\boldsymbol{u}_{j\rightarrow i}\right)\approx \mathop{\sum }_{\substack{ j=1 \\ j\neq i}}^{N}(\boldsymbol{u}_{j\rightarrow i}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}_{j\rightarrow i}). & & \displaystyle\end{eqnarray}$$

The vorticity of the undisturbed flow at the location of the $i$ th particle can be written as,

(2.6) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D74E}_{un,i}(x_{i})\approx \unicode[STIX]{x1D74E}_{mac}(x_{i})+\mathop{\sum }_{\substack{ j=1 \\ j\neq i}}^{N}\unicode[STIX]{x1D74E}_{j\rightarrow i}(x_{i}), & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D74E}_{mac}$ is the vorticity of the macroscale flow and $\unicode[STIX]{x1D74E}_{j\rightarrow i}(x_{i})$ is the perturbation vorticity due to the presence of the $j$ th particle at the location of the $i$ th particle. Pairwise superposition can be employed in a similar manner to evaluate any other undisturbed flow quantity of interest at any of the particle locations, as required in the force formula to be discussed below.

2.2 Force model

Now that we have established the definition of the undisturbed flow, and an approximation to it based on pairwise interaction, in this section we will obtain an expression for the hydrodynamic force on the $i$ th particle. By definition the undisturbed flow $\boldsymbol{u}_{un,i}(x_{i})$ is the flow that would exist in the absence of the $i$ th particle, and the corresponding pressure and viscous stresses at the $i$ th particle are $p_{un,i}(x_{i})$ and $\unicode[STIX]{x1D749}_{un,i}(x_{i})$ . The presence of the $i$ th particle will introduce a perturbation flow, since the $i$ th particle’s undisturbed flow is now required to satisfy no-slip and no-penetration conditions on the $i$ th particle. We will term this perturbation flow self-perturbation, in order to distinguish it from the perturbation flows due to the presence of the neighbouring particles that contributed to microscale component of $\boldsymbol{u}_{un,i}(x_{i})$ in (2.2). Let us denote this self-perturbation flow and its pressure and viscous stress fields to be $\boldsymbol{u}_{i}^{\prime }(x_{i})$ , $p_{i}^{\prime }(x_{i})$ and $\unicode[STIX]{x1D749}_{i}^{\prime }(x_{i})$ . Following Maxey & Riley (Reference Maxey and Riley1983) the hydrodynamic force on the $i$ th particle can be expressed as

(2.7) $$\begin{eqnarray}\displaystyle \boldsymbol{F}_{i}=\oint _{S_{p}}(-p_{un,i}(x_{i})\boldsymbol{I}+\unicode[STIX]{x1D749}_{un,i}(x_{i}))\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}A+\oint _{S_{p}}(-p_{i}^{\prime }(x_{i})\boldsymbol{I}+\unicode[STIX]{x1D749}_{i}^{\prime }(x_{i}))\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}A, & & \displaystyle\end{eqnarray}$$

where the integrals are over the surface $S_{p}$ of the $i$ th particle. The first term on the right is the force due to the undisturbed flow and will be denoted as $\boldsymbol{F}_{un,i}$ , and the second term on the right is the force due to the self-perturbation flow, which will be denoted as $\boldsymbol{F}_{i}^{\prime }$ . Using Gauss’s theorem the surface integral can be replaced by the volume integral to obtain

(2.8) $$\begin{eqnarray}\displaystyle \boldsymbol{F}_{i}=\oint _{V_{p}}(-\unicode[STIX]{x1D735}p_{un,i}(x_{i})+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}_{un,i}(x_{i}))\,\text{d}V+\oint _{V_{p}}(-\unicode[STIX]{x1D735}p_{i}^{\prime }(x_{i})+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}_{i}^{\prime }(x_{i}))\,\text{d}V. & & \displaystyle\end{eqnarray}$$

In the limit when all the neighbours are far away, the boundary conditions satisfied by the self-perturbation flow are (i) on the surface of the $i$ th particle the self-perturbation flow should counteract the undisturbed flow to satisfy no-slip and no-penetration boundary conditions and (ii) in the far field the self-perturbation flow must decay to zero. With these requirements $\boldsymbol{F}_{i}^{\prime }$ is traditionally written as the sum of quasi-steady, added-mass, Basset history and vorticity-induced force contributions. Expressions for these contributions have been rigorously obtained by Gatignol (Reference Gatignol1983) and Maxey & Riley (Reference Maxey and Riley1983) in the Stokes regime, with empirical corrections that extend their applicability to finite Reynolds numbers (Magnaudet & Eames Reference Magnaudet and Eames2000; Balachandar & Eaton Reference Balachandar and Eaton2010; Crowe et al. Reference Crowe, Schwarzkopf, Sommerfeld and Tsuji2011). It is important to note that these force expressions are based on the undisturbed flow at the location of the $i$ th particle, although the force $\boldsymbol{F}_{i}^{\prime }$ itself arises due to self-perturbation.

Here are we are interested in situations where the neighbours are close by. In which case, in addition to the boundary conditions (i) and (ii) stated above, the self-perturbation flow must also satisfy (iii) no-slip and no-penetration boundary conditions on the surface of all the neighbours. The self-perturbation flow $\boldsymbol{u}_{i}^{\prime }(x_{i})$ and the corresponding pressure and viscous stress fields will be affected by the presence of neighbours and therefore will deviate from what they would be for an isolated $i$ th particle. Thus, the neighbours play a dual role, first in perturbing the macroscale flow and defining the undisturbed flow, and second in modifying the self-perturbation flow induced by the $i$ th particle.

Based on the above line of reasoning, we can conceptually separate $\boldsymbol{F}_{i}^{\prime }$ into

(2.9) $$\begin{eqnarray}\displaystyle \boldsymbol{F}_{i}^{\prime }=\unicode[STIX]{x1D74D}_{i}^{\prime }F_{i0}^{\prime }, & & \displaystyle\end{eqnarray}$$

where $\boldsymbol{F}_{i0}^{\prime }$ is the force due to the self-perturbation flow, ignoring the presence of the neighbours (only in the evaluation of the self-perturbation flow) and the correction $\unicode[STIX]{x1D74D}_{i}^{\prime }$ accounts for the force modification due to changes in the self-perturbation flow because of the presence of neighbours. In Stokesian suspensions of finite numbers of particles, the corrections can be precisely obtained using method of reflections (Guazzelli & Morris Reference Guazzelli and Morris2011; Kim & Karrila Reference Kim and Karrila2013). Given the undisturbed flow at the $i$ th particle $\boldsymbol{F}_{i0}^{\prime }$ can be readily evaluated using the force expressions of Gatignol (Reference Gatignol1983) and Maxey & Riley (Reference Maxey and Riley1983) and their finite Reynolds number extensions. This is made possible due to the fact that $\boldsymbol{F}_{i0}^{\prime }$ has been defined to be force on an isolated particle. Unfortunately, Stokesian dynamics and related rigorous theoretical developments do not apply at the finite Reynolds numbers of present interest. Empirical correction functions that depend on local volume fraction must instead be used to account for the collective effect of neighbouring particles. For example, in the case of quasi-steady force, $\unicode[STIX]{x1D74D}_{i}^{\prime }$ would correspond to finite volume fraction correction to quasi-steady force, and so on. Such correction functions account for the presence of neighbours only in the mean sense, in terms of local volume fraction.

In summary, the force on a particle can be expressed by the following sum

(2.10) $$\begin{eqnarray}\displaystyle \boldsymbol{F}=\boldsymbol{F}_{un}+\boldsymbol{F}_{qs}+\boldsymbol{F}_{am}+\boldsymbol{F}_{bh}+\boldsymbol{F}_{L}, & & \displaystyle\end{eqnarray}$$

where we have suppressed the subscript $i$ . The terms on the right-hand side are respectively the undisturbed, quasi-steady, added-mass, Basset history (or viscous unsteady) and rotation-induced lift forces. All these forces can be expressed entirely in terms of the undisturbed flow characterized by $\boldsymbol{u}_{un,i}(x_{i})$ , $p_{un,i}(x_{i})$ and $\unicode[STIX]{x1D749}_{un,i}(x_{i})$ , without the knowledge of the precise location of the neighbours. Thus, at the present level of approximation, the PIEP model uses the precise location of the neighbouring particles only in obtaining the undisturbed flow, but the neighbours contribute to the calculation of the force due to the undisturbed flow only in an average sense. The rest of our discussion will focus attention to steady flows over an array of particles, and this allows us to further simplify the force expressions. In what follows we discuss the modelling of each of these force contributions.

2.2.1 Undisturbed flow force

The undisturbed flow force, also referred to as the pressure-gradient or the Archimedes force, is due to the pressure and viscous stress variation of the undisturbed flow over the location of the particle. Here the undisturbed flow over the $i$ th particle is defined as in (2.1) and thus includes the perturbation flow due to all the neighbours of the $i$ th particle. It is undisturbed only in the sense that it is the flow that would exist in the absence of the $i$ th particle. By definition, the undisturbed flow force is experienced even in the absence of the $i$ th particle, by the fluid volume that resides at the location of the $i$ th particle. From (2.8) the undisturbed flow force can be expressed as

(2.11) $$\begin{eqnarray}\displaystyle \boldsymbol{F}_{un}=V_{p}(-\unicode[STIX]{x1D735}p_{un,i}(\boldsymbol{x}_{i})+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}_{un,i}(\boldsymbol{x}_{i})), & & \displaystyle\end{eqnarray}$$

where in evaluating the volume integral, the spatial variation of pressure-gradient and viscous stress over the volume have been ignored for now. In § 2.3 we will use Faxén’s law to systematically account for such spatial variation and obtain the more general version of the undisturbed flow force and the other force contributions to be discussed below. Separating the undisturbed flow into the macroscale and microscale contributions we obtain

(2.12) $$\begin{eqnarray}\displaystyle V_{p}(-\unicode[STIX]{x1D735}p_{mac}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}_{mac})+V_{p}(-\unicode[STIX]{x1D735}p_{mic,i}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}_{mic,i}). & & \displaystyle\end{eqnarray}$$

Finally we apply the pairwise interaction approximation (2.3) for the microscale pressure and viscous stress tensor at the $i$ th particle to obtain

(2.13) $$\begin{eqnarray}\displaystyle \boldsymbol{F}_{un}\approx V_{p}(-\unicode[STIX]{x1D735}p_{mac}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}_{mac})+V_{p}\left[-\mathop{\sum }_{\substack{ j=1 \\ j\neq i}}^{N}\unicode[STIX]{x1D735}p_{j\rightarrow i}+\mathop{\sum }_{\substack{ j=1 \\ j\neq i}}^{N}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}_{j\rightarrow i}\right]. & & \displaystyle\end{eqnarray}$$

2.2.2 Quasi-steady force

In (2.10) all the force contributions other than $\boldsymbol{F}_{un}$ are due to the self-perturbation flow created by the no-slip and no-penetration boundary conditions imposed by the $i$ th particle. The quasi-steady force on the $i$ th particle can be expressed using the standard drag law (Schiller & Naumann Reference Schiller and Naumann1933) as

(2.14) $$\begin{eqnarray}\displaystyle \boldsymbol{F}_{qs}=\boldsymbol{F}_{ra}(\unicode[STIX]{x1D719},Re)+\{3\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}d(\boldsymbol{u}_{un,i}(x_{i})-v_{i})[1+0.15(Re)^{0.687}]-\boldsymbol{F}_{std}\}, & & \displaystyle\end{eqnarray}$$

where $\boldsymbol{F}_{ra}(\unicode[STIX]{x1D719},Re)$ is a suitable expression for the mean drag as a function of mean Reynolds number and mean volume fraction. In the above $v_{i}$ is the velocity of the $i$ th particle and the Reynolds number of the $i$ th particle is defined as,

(2.15) $$\begin{eqnarray}\displaystyle Re=\frac{d|\boldsymbol{u}_{mac}(\boldsymbol{x}_{i})-v_{i}|}{\unicode[STIX]{x1D708}}, & & \displaystyle\end{eqnarray}$$

with $\unicode[STIX]{x1D708}$ being the viscosity of the fluid and $d$ the particle diameter. In the above $\boldsymbol{F}_{std}$ is the standard drag on the $i$ th particle ignoring the influence of the neighbours at the microscale. It then depends only the macroscale flow and is given by

(2.16) $$\begin{eqnarray}\displaystyle \boldsymbol{F}_{std}=3\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}d(\boldsymbol{u}_{mac}(x_{i})-v_{i})[1+0.15(Re)^{0.687}]. & & \displaystyle\end{eqnarray}$$

Thus, the second term in (2.14) within the curly parenthesis accounts for the quasi-steady effect of the perturbation flow created by the neighbouring particles at the microscale (i.e. quasi-steady force due to $\boldsymbol{u}_{un,i}(x_{i})-\boldsymbol{u}_{mac}(x_{i})$ ). It should be noted that the standard drag law by Schiller and Naumann is valid for up to $Re_{p}\approx 1000$ . Any other finite Reynolds number drag law (Clift, Grace & Weber Reference Clift, Grace and Weber1978; Khan & Richardson Reference Khan and Richardson1987; Brown & Lawler Reference Brown and Lawler2003) can be used instead in (2.14).

The standard drag correlation used in both terms within the curly parenthesis of (2.14) is for an isolated particle. It has now been well established that mean drag force on a particle within a random bed of particles is substantially higher than what will be predicted by the standard drag relation (Hill, Koch & Ladd Reference Hill, Koch and Ladd2001; Tenneti et al. Reference Tenneti, Garg and Subramaniam2011). The increase in quasi-steady mean drag is due to the fact that the self-perturbation flow around a particle, even under steady conditions, is strongly influenced by the presence of other particles, and differs from that of an isolated particle. The average influence of finite volume fraction has been incorporated into enhanced drag law relations which are dependent on both the macroscale Reynolds number $(Re)$ and the macroscale volume fraction $(\unicode[STIX]{x1D719})$ . A few such drag laws, denoted as $\boldsymbol{F}_{ra}(\unicode[STIX]{x1D719},Re)$ , obtained from the literature (Beetstra et al. Reference Beetstra, Van der Hoef and Kuipers2007; Tenneti et al. Reference Tenneti, Garg and Subramaniam2011; Zaidi et al. Reference Zaidi, Tsuji and Tanaka2014; Bogner et al. Reference Bogner, Mohanty and Rüde2015; Tang et al. Reference Tang, Peters, Kuipers, Kriebitzsch and Hoef2015) are plotted alongside (2.16) in figure 2 for volume fractions of 0.1 and 0.21. (The subscript ‘ra’ denotes the mean drag correlation for a random distribution of particles. The results for structured arrays of particles, such as body-centred cubic or face-centred cubic arrangements, will have different mean drag correlation.) Although there is significant difference between the different finite volume fraction drag laws, all the correlations suggest a higher drag than the standard drag law due to the intrinsic effects within the array, which are amplified with increasing volume fraction. In the results to be presented we will use the mean drag correlation by Tenneti et al. (Reference Tenneti, Garg and Subramaniam2011) for $\boldsymbol{F}_{ra}(\unicode[STIX]{x1D719},Re)$ but qualitatively similar results are obtained with the use of other finite volume fraction drag laws.

Figure 2.

Drag laws for finite volume fractions versus the superficial Reynolds number for (a) $\unicode[STIX]{x1D719}=0.1$ and (b) $\unicode[STIX]{x1D719}=0.21$ . The standard drag law (SDL) for a single sphere is plotted alongside for comparison.

The quasi-steady drag relation (2.14) can now be interpreted in the following way. Given the macroscale information (i.e. local volume fraction at the macroscale, Reynolds number based on the macroscale flow) the best estimate of the mean quasi-steady drag is given by $\boldsymbol{F}_{ra}(\unicode[STIX]{x1D719},Re)$ . The purpose of the second term is to account for the perturbation to the quasi-steady force due to changes in the undisturbed flow at the location of the $i$ th particle induced by the local arrangement of neighbouring particles. When (2.2) and (2.16) are substituted for $\boldsymbol{u}_{un,i}$ and $\boldsymbol{F}_{std}$ , the second term within the curly parenthesis in (2.14) becomes,

(2.17) $$\begin{eqnarray}\displaystyle 3\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}d\left(\mathop{\sum }_{\substack{ j=1 \\ j\neq i}}^{N}\boldsymbol{u}_{j\rightarrow i}(x_{i})\right)[1+0.15(Re)^{0.687}]. & & \displaystyle\end{eqnarray}$$

The above perturbation force will approach zero when averaged over a large number of particles and over an ensemble of realizations. From the definition of perturbation velocity we expect

(2.18) $$\begin{eqnarray}\displaystyle \left\langle \mathop{\sum }_{\substack{ j=1 \\ j\neq i}}^{N}\boldsymbol{u}_{j\rightarrow i}(x_{i})\right\rangle \rightarrow 0, & & \displaystyle\end{eqnarray}$$

where the angle brackets represent an average over an ensemble of realizations. This will ensure that the mean drag force for a steady macroscale flow will correctly reduce to $\boldsymbol{F}_{ra}(\unicode[STIX]{x1D719},Re)$ .

Finally, the expression (2.14) for the quasi-steady force is not the only acceptable form. For example, the definition of the Reynolds number of the $i$ th particle can be based on undisturbed fluid velocity, as opposed to the macroscale velocity. In other words, $Re_{i}$ in (2.14) can be replaced by $d|\boldsymbol{u}_{un,i}(\boldsymbol{x}_{i})-v_{i}|/\unicode[STIX]{x1D708}$ . With this modification, however, it is not clear if the perturbation force within the curly parenthesis on the right-hand side of (2.14) summed over an ensemble of realizations will reduce to zero.

2.2.3 Added-mass force

The added-mass force arises from the relative acceleration of the undisturbed flow at the particle location with respect to the acceleration of the particle. This force is expressed as

(2.19) $$\begin{eqnarray}\displaystyle \boldsymbol{F}_{am} & = & \displaystyle C_{m}m_{f}\left(\frac{\text{D}\boldsymbol{u}_{un,i}(x_{i})}{\text{D}t}-\frac{\text{d}v_{i}}{\text{d}t}\right)\nonumber\\ \displaystyle & = & \displaystyle C_{m}m_{f}\left(\frac{\text{D}\boldsymbol{u}_{mac}(x_{i})}{\text{D}t}-\frac{\text{d}v_{i}}{\text{d}t}\right)+C_{m}V_{p}(-\unicode[STIX]{x1D735}p_{mic,i}(x_{i})+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}_{mic,i}(x_{i}))\nonumber\\ \displaystyle & {\approx} & \displaystyle C_{m}m_{f}\left[\frac{\text{D}\boldsymbol{u}_{mac}(x_{i})}{\text{D}t}-\frac{\text{d}v_{i}}{\text{d}t}\right]+C_{m}V_{p}\left[-\mathop{\sum }_{\substack{ j=1 \\ j\neq i}}^{N}\unicode[STIX]{x1D735}p_{j\rightarrow i}(x_{i})+\mathop{\sum }_{\substack{ j=1 \\ j\neq i}}^{N}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}_{j\rightarrow i}(x_{i})\right].\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

Note that both the flow acceleration as well as the acceleration of the $i$ th particle contribute to the added-mass force. Again, the flow acceleration has been separated into macroscale and microscale contributions. We have used the pairwise interaction approximation (2.2) to simplify the microscale contribution. In the present problem of a steady flow over a stationary bed of particles, $\text{d}v_{i}/\text{d}t=0$ and $\unicode[STIX]{x2202}\boldsymbol{u}_{un,i}/\unicode[STIX]{x2202}t=0$ . But the total acceleration of the undisturbed flow at the location of the $i$ th particle may not be zero due to convective acceleration (i.e.  $\text{D}\boldsymbol{u}_{un,i}(x_{i})/\text{D}t$ is equal to $\unicode[STIX]{x2202}\boldsymbol{u}_{un,i}(x_{i})\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}_{un,i}(x_{i})$ which may be non-zero). Thus, $\boldsymbol{F}_{am}$ can be non-zero even in a steady flow through an array of particles. This force arises from the self-perturbation flow that results from the imposition of no-penetration boundary condition on the surface of the $i$ th particle. In the above equation $C_{m}$ is the added-mass coefficient and is equal to $1/2$ for an isolated sphere. The added-mass coefficient has been shown to increase with particle volume fraction. Zuber (Reference Zuber1964) and Sangani, Zhang & Prosperetti (Reference Sangani, Zhang and Prosperetti1991) have shown the added-mass coefficient to scale as $(1+2\unicode[STIX]{x1D719}_{mac})/2$ for small values of volume fraction.

2.2.4 Basset history force

The Basset history force arises from the viscous effect of relative acceleration of the undisturbed flow at the particle location with respect to the acceleration of the particle. This force is also known as the viscous unsteady force and it can be expressed as a convolution integral (Maxey & Riley Reference Maxey and Riley1983; Magnaudet & Eames Reference Magnaudet and Eames2000; Balachandar & Eaton Reference Balachandar and Eaton2010)

(2.20) $$\begin{eqnarray}\displaystyle \boldsymbol{F}_{bh}=\frac{3}{2}d^{2}\sqrt{\unicode[STIX]{x03C0}\unicode[STIX]{x1D70C}\unicode[STIX]{x1D707}}\int _{-\infty }^{t}K_{bh}(t-\unicode[STIX]{x1D70F})\left(\frac{\text{d}\boldsymbol{u}_{un,i}}{\text{d}t}-\frac{\text{d}v_{i}}{\text{d}t}\right)_{@\unicode[STIX]{x1D70F}}\,\text{d}\unicode[STIX]{x1D70F}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D70C}$ and $\unicode[STIX]{x1D707}$ are the density and viscosity of the fluid, respectively, and the subscript $@\unicode[STIX]{x1D70F}$ indicates that the relative acceleration is evaluated at time equal to $\unicode[STIX]{x1D70F}$ . The Basset history kernel $K_{bh}$ has a $1/\sqrt{(t-\unicode[STIX]{x1D70F})}$ behaviour for small values of the argument, but for larger values of the argument the kernel decays more rapidly as $1/(t-\unicode[STIX]{x1D70F})^{2}$ or faster (Mei & Adrian Reference Mei and Adrian1992; Lovalenti & Brady Reference Lovalenti and Brady1993; Kim, Elghobashi & Sirignano Reference Kim, Elghobashi and Sirignano1998). Unfortunately, the precise form of the kernel somewhat depends on the nature of the problem. Nevertheless in the present problem, according to the definition given above, the Basset history force is identically zero for all the particles. This is because $\text{d}\boldsymbol{u}_{un,i}/\text{d}t$ is the total derivative of the undisturbed flow at the $i$ th particle following the particle, which is identically zero since we consider a steady flow over a stationary distribution of particles.

Here we have followed Maxey & Riley (Reference Maxey and Riley1983) and Balachandar & Eaton (Reference Balachandar and Eaton2010) and defined fluid acceleration to be different in the added-mass and Basset history forces (i.e.  $\text{D}\boldsymbol{u}_{un,i}/\text{D}t$ in added mass versus $\text{d}\boldsymbol{u}_{un,i}/\text{d}t$ in Basset history force). As pointed out by Maxey & Riley (Reference Maxey and Riley1983) their difference can be ignored in the linear limit, but not in the finite Reynolds number regime. Solution of the linearized Navier–Stokes equations for the unsteady Stokes flow over an arbitrarily moving particle by Gatignol (Reference Gatignol1983) and Maxey & Riley (Reference Maxey and Riley1983) yielded $\text{d}\boldsymbol{u}_{un,i}/\text{d}t$ for both the unsteady forces. Based on the inviscid analysis of Auton (Reference Auton1987) and Auton, Hunt & Prud’Homme (Reference Auton, Hunt and Prud’Homme1988) the added-mass force in (2.19) is written in terms of $\text{D}\boldsymbol{u}_{un,i}/\text{D}t$ . For now, we will consider the Basset history force to be zero in the present steady flow, but later in the discussion section we will consider the effect of alternate forms of the Basset history force.

2.2.5 Rotational force

Even when the ambient flow at the macroscale $\boldsymbol{u}_{mac}$ is irrotational, the $i$ th particle will be subjected to a vortical flow, since it is immersed in the perturbation flow of the neighbouring particles. A simple example is the $i$ th particle sitting in the vortical wake of an upstream neighbour. The hydrodynamic force on a particle due to ambient vorticity can be separated into an inviscid contribution and a viscous contribution. A simple model of the rotational force is

(2.21) $$\begin{eqnarray}\displaystyle F_{L}\approx C_{L}m_{f}[(\boldsymbol{u}_{un,i}(\boldsymbol{x}_{i})-v_{i})\times (\unicode[STIX]{x1D714}_{un,i}(\boldsymbol{x}_{i})+2\unicode[STIX]{x1D734}_{i})], & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D734}_{i}$ is the angular velocity of the particle and $C_{L}$ is the lift coefficient. In the above equation the undisturbed flow velocity and vorticity can be approximated with the pairwise interaction model given in (2.2) and (2.6). From the analysis of Auton et al. (Reference Auton, Hunt and Prud’Homme1988) $C_{L}$ is $1/2$ for an isolated stationary sphere in a steady weak ambient vorticity field. In the more general case of a translating and rotating sphere in an ambient vortical field, the effects of ambient vorticity and particle rotation superpose, and alternate expressions for the lift coefficient have been proposed (Rubinow & Keller Reference Rubinow and Keller1961; Saffman Reference Saffman1965). The effects of finite Reynolds number, finite volume fraction, or time dependence further complicates the rotational or lift force model (Bagchi & Balachandar Reference Bagchi and Balachandar2002a ,Reference Bagchi and Balachandar b ; Wakaba & Balachandar Reference Wakaba and Balachandar2005). Here we have chosen to use the simple expression (2.21) with $C_{L}=1/2$ to model the rotational force. As we will see in the results section, the rotational force does not make a significant contribution to the overall force model.

2.3 Faxén form of the forces

The undisturbed flow at the $i$ th particle $\boldsymbol{u}_{un,i}(x_{i})$ can be expected to substantially vary on the scale of the $i$ th particle. However, the different force contributions given in (2.13), (2.14), (2.19) and (2.21) implicitly assume the undisturbed flow quantities, such as $\boldsymbol{u}_{un,i}(x_{i})$ , $p_{un,i}(x_{i})$ , $\unicode[STIX]{x1D749}_{un,i}(x_{i})$ and $\unicode[STIX]{x1D74E}_{un,i}(x_{i})$ , to be defined at the centre of the $i$ th particle. For the present non-uniform undisturbed flow we will consider the generalization of the Faxén’s theorem (Reference Faxén1922) by Gatignol (Reference Gatignol1983) and Maxey & Riley (Reference Maxey and Riley1983). For a spatially varying ambient flow the undisturbed flow force can be generalized as

(2.22) $$\begin{eqnarray}\displaystyle \boldsymbol{F}_{un}=V_{p}(-\overline{\unicode[STIX]{x1D735}p_{mac}}^{V}+\overline{\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}_{mac}}^{V})+V_{p}\left(\mathop{\sum }_{\substack{ j=1 \\ j\neq i}}^{N}[-\overline{\unicode[STIX]{x1D735}p_{j\rightarrow i}}^{V}+\overline{\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}_{j\rightarrow i}}^{V}]\right), & & \displaystyle\end{eqnarray}$$

where $\overline{()}^{V}$ indicates an average over the volume of the $i$ th particle. By Gauss’s theorem, $\boldsymbol{F}_{un}$ correctly reduces to the integral of the undisturbed pressure and viscous stresses over the surface of the particle. In obtaining the last equation we have commuted the summation operation over the neighbours and the volume integral. In the final form, the undisturbed flow force due to the macroscale flow and that due to the perturbation of each neighbouring particle has been separated.

The Faxén form of the quasi steady force can be expressed as

(2.23) $$\begin{eqnarray}\displaystyle \boldsymbol{F}_{qs}=\boldsymbol{F}_{ra}(\unicode[STIX]{x1D719},Re)+\left\{3\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}d\left(\mathop{\sum }_{\substack{ j=1 \\ j\neq i}}^{N}\overline{\boldsymbol{u}_{j\rightarrow i}(x_{i})}^{S}\right)[1+0.15(Re)^{0.687}]\right\}, & & \displaystyle\end{eqnarray}$$

where $\overline{()}^{S}$ indicates an average over the surface of the $i$ th particle and $S_{p}=4\unicode[STIX]{x03C0}R^{2}$ is the surface area of the $i$ th particle. Following the derivation of the undisturbed fluid force presented above, the Faxén form of the added-mass force can be expressed as

(2.24) $$\begin{eqnarray}\displaystyle \boldsymbol{F}_{am}=C_{m}m_{f}\left[\frac{\overline{\text{D}\boldsymbol{u}_{mac}(x_{i})}^{V}}{\text{D}t}-\frac{\text{d}\boldsymbol{v}_{i}}{\text{d}t}\right]+C_{m}V_{p}\left(\mathop{\sum }_{\substack{ j=1 \\ j\neq i}}^{N}[-\overline{\unicode[STIX]{x1D735}p_{j\rightarrow i}}^{V}+\overline{\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}_{j\rightarrow i}}^{V}]\right), & & \displaystyle\end{eqnarray}$$

where we have used the fact that the translational acceleration of the particle remains the same over the entire particle.

Finally, due to its nonlinear dependence on velocity, the rotational force can be extended in several different ways in the case of spatially varying undisturbed flow. The form that has been used in the results to be presented below is

(2.25) $$\begin{eqnarray}\displaystyle \boldsymbol{F}_{L}\approx C_{ir}m_{f}\left[\overline{\boldsymbol{u}_{mac}(x_{i})}^{V}+\mathop{\sum }_{\substack{ j=1 \\ j\neq i}}^{N}\overline{\boldsymbol{u}_{j\rightarrow i}(x_{i})}^{V}-v_{i}\right]\times \left[\overline{\unicode[STIX]{x1D74E}_{mac}(x_{i})}^{V}+\mathop{\sum }_{\substack{ j=1 \\ j\neq i}}^{N}\overline{\unicode[STIX]{x1D714}_{j\rightarrow i}(x_{i})}^{V}+2\unicode[STIX]{x1D734}_{i}\right]. & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

Note that in the above definition (i) we first average velocity and vorticity perturbations over the volume of the $i$ th particle, (ii) the volume averages are then summed over all the neighbours and added to the corresponding macroscale contribution and (iii) finally the cross-product of the net velocity and vorticity contribute to the rotational force. But the order of these operations can be rearranged to get a different definition of the rotational force. Also the above definition using volume averages is somewhat arbitrary. Instead, the volume averages can be replaced by averages over the surface of the $i$ th particle. Unfortunately, the rigorous theory of Gatignol (Reference Gatignol1983) and Maxey & Riley (Reference Maxey and Riley1983) is based on linear approximation and does not provide the precise form of the nonlinear rotational force. We have tried different forms of the inviscid rotational force, and since the contribution of $\boldsymbol{F}_{L}$ is generally small, the results to be presented are not sensitive to the precise form.

In summary, the PIEP model force is given by,

(2.26) $$\begin{eqnarray}\displaystyle \boldsymbol{F}=\boldsymbol{F}_{un}+\boldsymbol{F}_{qs}+\boldsymbol{F}_{am}+\boldsymbol{F}_{L}, & & \displaystyle\end{eqnarray}$$

where the four force contributions on the right-hand side are given by (2.22), (2.23), (2.24) and (2.25), respectively.

3.1 Numerical methodology

An isolated sphere in a uniform ambient flow presents a relatively simple geometry. At Reynolds numbers below 200 the flow remains steady and axisymmetric, thus allowing highly resolved simulations with a two-dimensional polar grid fitted around the sphere. Here we use a fully spectral simulation with Fourier representation in the circumferential direction and Chebyshev representation along the radial direction (Bagchi & Balachandar Reference Bagchi and Balachandar2002a ,Reference Bagchi and Balachandar b ). The computational domain is chosen to be a larger sphere of radius 15 diameters, which has been shown to be sufficiently large in prior investigations of flow over a sphere. A Chebyshev–Gauss–Lobatto distribution of grid points is used in the radial direction, starting from the definition

(3.1) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D709}_{i}=-\cos \left(\frac{\unicode[STIX]{x03C0}(i-1)}{N_{r}-1}\right), & & \displaystyle\end{eqnarray}$$

where $N_{r}$ is the number of grid points in the radial direction. A transformation is then implemented to stretch the grid and allocate more points in the region close to the sphere where velocity gradients are more prominent. The transformation is given by,

(3.2) $$\begin{eqnarray}\displaystyle \hat{\unicode[STIX]{x1D709}}_{i}=C_{0}+C_{1}\unicode[STIX]{x1D709}_{i}-C_{0}\unicode[STIX]{x1D709}_{i}^{2}+C_{3}\unicode[STIX]{x1D709}_{i}^{3}. & & \displaystyle\end{eqnarray}$$

By requiring $\hat{\unicode[STIX]{x1D709}}(\unicode[STIX]{x1D709}=1)=1$ and $\hat{\unicode[STIX]{x1D709}}(\unicode[STIX]{x1D709}=-1)=-1$ , we constrain the constants $C_{1}$ and $C_{3}$ to,

(3.3) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle C_{1}=0.5(-\unicode[STIX]{x1D6FE}_{1}+2C_{0}+3),\\ \displaystyle C_{3}=0.5(\unicode[STIX]{x1D6FE}_{1}-2C_{0}-1).\end{array}\right\} & & \displaystyle\end{eqnarray}$$

The extent of the stretch is thus determined by the constants $C_{0}$ and $\unicode[STIX]{x1D6FE}_{1}$ , which are chosen to be 0.6 and 0.2 in the present simulations. The grid points in the radial direction are then scaled to cover the whole domain from the radius of the particle $(r_{in})$ to the radius of the outer spherical domain $(r_{out})$ ,

(3.4) $$\begin{eqnarray}\displaystyle r_{i}={\textstyle \frac{1}{2}}\hat{\unicode[STIX]{x1D709}}_{i}[r_{in}-r_{out}]+{\textstyle \frac{1}{2}}[r_{in}+r_{out}],\quad \text{for }i=1,2,\ldots ,N_{r}. & & \displaystyle\end{eqnarray}$$

By symmetry the grid points in the tangential direction are defined from 0 to $\unicode[STIX]{x03C0}$ . Here again the grid points are stretched and clustered in the wake region to provide increased resolution of the wake vortical structures:

(3.5) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D703}(j)=\tan ^{-1}\left(\frac{\sin \displaystyle \left(\frac{\unicode[STIX]{x03C0}}{N_{\unicode[STIX]{x1D703}}}(j-1/2)\right)(1-h^{2})}{\cos \displaystyle \left(\frac{\unicode[STIX]{x03C0}}{N_{\unicode[STIX]{x1D703}}}(j-1/2)\right)(1+h^{2})-2h}\right),\quad \text{for }j=1,2,\ldots ,N_{\unicode[STIX]{x1D703}}, & & \displaystyle\end{eqnarray}$$

where $N_{\unicode[STIX]{x1D703}}$ is the number of grid points in the tangential direction. The grid is stretched towards the downstream direction by specifying a negative parameter $h$ . In the present simulations $h$ is set to $-0.35$ which provides a higher resolution in the wake of the sphere. A schematic of the grid is shown in figure 3.

Figure 3.

Body fitted spherical mesh used for the single-particle simulations. $N_{r}=121$ , $N_{\unicode[STIX]{x1D703}}=100$ , $C_{0}=0.6$ , $\unicode[STIX]{x1D6FE}_{1}=0.2$ and $h=-0.35$ . The arrows indicate the direction of the flow.

A Dirichlet boundary condition for a uniform flow of unit non-dimensional magnitude is imposed at the inlet. A non-reflecting boundary condition is used at the outflow by parabolizing the governing equations (Mittal & Balachandar Reference Mittal and Balachandar1996). The Fourier–Chebyshev spectral method is used to solve the governing equations given by the non-dimensional incompressible Navier–Stokes equations,

(3.6) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$

(3.7) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}=\unicode[STIX]{x1D735}p+\frac{1}{Re}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}. & \displaystyle\end{eqnarray}$$

A mixed second-order Adams–Bashforth and Crank–Nicolson scheme is used for the time discretization. A resolution of $N_{r}=121$ and $N_{\unicode[STIX]{x1D703}}=100$ is used. A more detailed description on the numerical scheme can be found in Bagchi & Balachandar (Reference Bagchi and Balachandar2002a ). In the above equations all the variables are non-dimensional, with the particle diameter as the length scale and the uniform inflow velocity as the velocity scale. Correspondingly, pressure is scaled by density times the square of uniform inflow velocity. Henceforth we will consider only non-dimensional quantities (unless explicitly stated otherwise).

3.2 Pairwise interaction force maps

Using the above described Fourier–Chebyshev spectral method, flow around an isolated sphere was simulated for a range of Reynolds numbers from 10 to 215. A sample result for $Re=87$ is shown in figure 4, where frames (a–d) present contour plots of streamwise $(x)$ velocity, cross-stream $(y)$ velocity, pressure and circumferential vorticity. Note that the mean value of non-dimensional streamwise velocity is unity, while the mean value of all other quantities are zero. The macroscale flow being considered is $\boldsymbol{u}_{mac}=1$ and $\unicode[STIX]{x1D74E}_{mac}=0$ . The steady velocity fields, pressure and vorticity shown in figure 4 represent well known results. Here they have been obtained at high resolution, so that they can be used in the pairwise interaction model.

Figure 4.

Contour plots of the (a) streamwise velocity, (b) spanwise velocity, (c) pressure and (d) vorticity obtained for a uniform ambient flow over a sphere (shaded black) at $Re_{i}=87$ . The plots are on the $x$ – $y$ plane passing through the centre of the sphere. These results are used in the evaluation of the perturbation flow. The black shaded sphere is the $j$ th neighbouring sphere (source of the perturbation) while the open circle represents the $i$ th sphere, where the effect of the perturbation will be evaluated. The open circle represents a virtual sphere which does not affect the contour plots shown.

The central sphere in figure 4 represents the $j$ th particle and the quantities plotted are the perturbation flow due to the presence of the $j$ th particle in the uniform macroscale flow. As an example, now consider a particular location of the $i$ th particle as represented by the thin black circle in figure 4. With respect to the centre of this $i$ th particle the $j$ th particle is 1.5 non-dimensional units upstream and 0.75 units shifted down along the $y$ axis. The volume and surface integrals of the various perturbation quantities defined in the Faxén form of the force expressions (2.22)–(2.25) can now be evaluated. For the specific location of the $i$ th particle with respect to the $j$ th particle shown in figure 4 we obtain

(3.8) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle \overline{\boldsymbol{u}_{j\rightarrow i}}^{S}=-0.19\boldsymbol{e}_{x}-0.04\boldsymbol{e}_{y},\quad \overline{\boldsymbol{u}_{j\rightarrow i}}^{V}=-0.18\boldsymbol{e}_{x}-0.05\boldsymbol{e}_{y},\\ \displaystyle \overline{\unicode[STIX]{x1D735}p_{j\rightarrow i}}^{V}=-0.04\boldsymbol{e}_{x}-(1.9\times 10^{-3})\boldsymbol{e}_{y},\\ \displaystyle \overline{\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}_{j\rightarrow i}}^{V}=-0.002\boldsymbol{e}_{x}+0.0027\boldsymbol{e}_{y},\quad \overline{\unicode[STIX]{x1D74E}_{j\rightarrow i}}^{V}=-0.74\boldsymbol{e}_{z}.\end{array}\right\} & & \displaystyle\end{eqnarray}$$

The position of the $i$ th particle with respect to the $j$ th particle is systematically varied and the requisite volume or surface integrals are computed from the single-particle DNS perturbation field. The contours of the resulting streamwise $(x)$ and lateral $(y)$ components of surface-averaged velocity $\overline{\boldsymbol{u}_{j\rightarrow i}}^{S}$ are shown in figure 5. The surface average can be considered as a filter operation and therefore the contour plots in figure 5 are similar to those in figure 4. In figure 5 the value of the surface average is plotted at the centre of the $i$ th particle. Thus, the contour plots are not shown in the region one diameter from the centre of the $j$ th particle, since the centres of the two identical spheres cannot come closer than one diameter.

Figure 5.

Mapping of the average velocity perturbation $\overline{\boldsymbol{u}_{j\rightarrow i}}^{S}$ on the surface of a virtual sphere in the (a) streamwise and (b) spanwise direction. $Re_{i}=87$ .

As can be expected the streamwise velocity perturbation at the $i$ th particle is negative if it is either in the wake or upstream of the $j$ th particle, while the perturbation is weak and positive if it is located laterally. The lateral velocity perturbations at the $i$ th particle are stronger from a particle located immediately downstream, but shifted off-axis with respect to the macroscale flow direction. The lateral velocity perturbations due to a downstream located particle is relatively weak. These maps are dependent only on the Reynolds number. Here we show the map for $Re=87$ as an example, and the qualitative behaviour remains the same over the Reynolds number range considered. The $x$ and $y$ extents of the contour plots are chosen such that the surface-averaged perturbation velocity falls to 0.5 % of the incoming uniform ambient flow. Neighbouring spheres located outside these limits would have a negligible effect.

The contour plots of the streamwise and lateral components of the pressure-gradient force integrated over the surface of the particle, $-S_{p}\overline{p_{j\rightarrow i}\boldsymbol{n}}^{S}=-V_{p}\overline{\unicode[STIX]{x1D735}\boldsymbol{\cdot }p_{j\rightarrow i}}^{V}$ , are shown in figure 6 as a function of streamwise and lateral separation of the $i$ th particle from the $j$ th particle (i.e. as a function of $r_{x}$ and $r_{y}$ ). This figure is similar to figure 5 in that the value of the contour for any value of $r_{x}$ and $r_{y}$ corresponds to pressure force integrated over the surface of the $i$ th particle whose centre is located at $(r_{x},r_{y})$ , where the perturbation flow is due to the $j$ th particle whose centre is located at the origin. These figures can be easily interpreted in conjunction with the perturbation pressure contours shown in figure 4(c). The large stagnation pressure at the front of a downstream neighbour will tend to push the $i$ th particle upstream. Similarly the low pressure region in the wake of a downstream neighbour will also tend to push the $i$ th particle upstream. The pressure effect of a neighbour located laterally will tend to push the $i$ th particle downstream.

Figure 6.

Mapping of the pressure-gradient force $-V_{p}\overline{\unicode[STIX]{x1D735}\boldsymbol{\cdot }p_{j\rightarrow i}}^{V}$ on a virtual sphere in the (a) streamwise and (b) spanwise direction. $Re_{i}=87$ .

In terms of lateral forces, the effect of perturbation pressure due a downstream neighbour is to push the $i$ th particle off axis. While the effect of perturbation pressure due an upstream neighbour is to push the $i$ th particle towards the axis. The contour plots of the streamwise and lateral components of the viscous force integrated over the surface of the particle, $S_{p}\overline{\unicode[STIX]{x1D749}_{j\rightarrow i}\boldsymbol{\cdot }\boldsymbol{n}}^{S}=V_{p}\overline{\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}_{j\rightarrow i}}^{V}$ , are shown in figure 7 as a function of $r_{x}$ and $r_{y}$ . In general, these force contributions are smaller in magnitude than the pressure contribution. From figures 6 and 7 it can be seen that the undisturbed flow force due to perturbation flow remains significant only for neighbours located within two diameters.

Figure 7.

Mapping of the stress gradient force $V_{p}\overline{\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}_{j\rightarrow i}}^{V}$ on a virtual sphere in the (a) streamwise and (b) spanwise direction. $Re_{i}=87$ .

Finally in figure 8 we present contours of volume-averaged vorticity, $\overline{\unicode[STIX]{x1D714}_{j\rightarrow i}(x_{i})}^{V}$ , and the streamwise and lateral components of rotational force due to the perturbation flow given by $[\overline{\boldsymbol{u}_{mac}}^{V}+\overline{\boldsymbol{u}_{j\rightarrow i}}^{V}]\times \overline{\unicode[STIX]{x1D74E}_{j\rightarrow i}}^{V}$ . As can be expected, rotational force is significant only in case when the $i$ th particle is in the vortical wake of an upstream neighbour. Also, the lateral component of the inviscid rotational force is much stronger than the streamwise component. The streamwise component tends to push the $i$ th particle downstream, while the lateral force tends to move the particle away from the axis.

Figure 8.

Mapping of the (a) vorticity volume average, (b) inviscid rotational force volume average in the streamwise direction and (c) spanwise direction. $Re_{i}=87$ .

All the force maps shown in figures 5–8 are for a Reynolds number of 87. The general qualitative description of these maps hold for other Reynolds numbers, but their magnitude and spatial extent somewhat vary. As an example, contours of the streamwise and spanwise component of the undisturbed flow force, $S_{p}(-\overline{p_{j\rightarrow i}\boldsymbol{n}}^{S}+\overline{\unicode[STIX]{x1D749}_{j\rightarrow i}\boldsymbol{\cdot }\boldsymbol{n}}^{S})$ , as a function of $r_{x}$ and are shown in figure 9 for four different Reynolds numbers. The variation can be seen to be continuous, but complex, for varying Reynolds number. Because of the axisymmetric nature of the simulations it is relatively inexpensive to create these maps for different Reynolds numbers. With small increments of $Re$ , a library of the maps shown in figures 5–9 can be created and stored for future application of pairwise interaction model.

Figure 9.

Mappings of the inertial force on a virtual sphere in the (a) streamwise and (b) lateral direction for different Reynolds number.

The force maps shown in figures 5–9 are due to one neighbouring particle. For undisturbed flow force, quasi-steady force and added-mass force, from their definitions given in (2.22), (2.23), and (2.24), contributions from each neighbour can be simply added. However, according to its definition given in (2.25), the rotational force contributions from the individual neighbours cannot be simply added. Velocity and vorticity perturbations from individual neighbours must be first added before taking their cross-product.