A numerical approach for particle-vortex interactions based on volume-averaged equations

TL;DR

This paper introduces a volume-averaged numerical model for particle-vortex interactions in turbulent flows, accurately capturing history effects and flow disturbances with lower computational cost than traditional methods.

Contribution

The authors extend previous volume-averaged fluid equations to better estimate fluid forces, validated against detailed interface-resolved simulations, improving accuracy for particles comparable to flow eddies.

Findings

The model accurately captures history effects and lift forces.

VA simulation results closely match fully-resolved simulations for high density ratios.

Flow disturbances and particle trajectories are well represented, with reduced computational cost.

Abstract

To study the dynamics of particles in turbulence when their sizes are comparable to the smallest eddies in the flow, the Kolmogorov length scale, efficient and accurate numerical models for the particle-fluid interaction are still missing. Therefore, we here extend the treatment of the particle feedback on the fluid based on the volume-averaged fluid equations (VA simulation) in the previous study of the present authors, by estimating the fluid force correlated with the disturbed flow. We validate the model against interface-resolved simulations using the immersed-boundary method. Simulations of single particles show that the history effect is well captured by the present estimation method based on the disturbed flow. Similarly, the simulation of the flow around a rotating particle demonstrates that the lift force is also well captured by the proposed method. We also consider the interaction between non-negligible size particles and an array of Taylor-Green vortices. For density ratios $\rho_d/\rho_c\geq$ 10, the results show that the particle motion captured by the VA approach is closer to that of the fully-resolved simulations than that obtained with a traditional two-way coupling simulation. The flow disturbance is also well represented by the VA simulation. In particular, it is found that history effects enhance the curvature of the trajectory in vortices and this enhancement increases with the particle size. Furthermore, the flow field generated by a neighboring particle at distances of around ten particle diameters significantly influences particle trajectories. The computational cost of the VA simulation proposed here is considerably lower than that of the interface-resolved simulation.

Tables

Table 1. Table 1: Acceleration obtained by the present estimation method. The subscripts 1 and 2 correspond to x 1 subscript 𝑥 1 x_{1} and x 2 subscript 𝑥 2 x_{2} -components.

${\rho }_d/{\rho }_c$	$a_1^{sv}$	$a_2^{sv}$	$a_1^{sv+pg}$	$a_2^{sv+pg}$	$a_1^{FR}$	$a_2^{FR}$
1	$2.50$	$-0.122$	$2.35$	$-18.4$	$0.477$	$-19.2$
10	$-1.25$	$-4.85$	$-1.62$	$-7.24$	$-1.52$	$-7.69$
100	$0.292$	$-0.314$	$0.464$	$-0.568$	$0.547$	$-0.628$
1000	$-0.174$	$0.174$	$-0.176$	$0.176$	$-0.196$	$0.196$

Equations

$${\alpha }_k$$

$${⟨B⟩}_k$$

$$\nabla \cdot \mathbf{u}=0,$$

$$\frac{\partial \mathbf{u}}{\partial t}+\nabla \cdot (\mathbf{u}\mathbf{u})=-\frac{1}{{\rho }_c}\nabla p+\nu {\nabla }^2\mathbf{u}+\mathbf{h},$$

$$\nabla \cdot ({\alpha }_c{⟨\mathbf{u}⟩}_c+{\alpha }_d{⟨\mathbf{v}⟩}_d)=0,$$

$$\frac{\partial ({\alpha }_c{⟨\mathbf{u}⟩}_c)}{\partial t}+\nabla \cdot ({\alpha }_c{⟨\mathbf{u}⟩}_c{⟨\mathbf{u}⟩}_c)=-\frac{1}{{\rho }_c}\nabla P+\nu {\nabla }^2({\alpha }_c{⟨\mathbf{u}⟩}_c+{\alpha }_d{⟨\mathbf{v}⟩}_d)+{\alpha }_c{⟨\mathbf{h}⟩}_c-\nabla \cdot \mathbf{\tau }+\frac{\mathbf{f}}{V},$$

$$\mathbf{\tau }={\alpha }_c{⟨\delta \mathbf{u}\delta \mathbf{u}⟩}_c,$$

$$\mathbf{\tau }$$

$$\mathbf{f}=-{\int }_{S_d}\left\{-\frac{1}{{\rho }_c}\delta p\mathbf{n}+\nu \left(\nabla \mathbf{u}+(\nabla \mathbf{u}{)}^T\right)\cdot \mathbf{n}\right\}dS$$

$$\mathbf{f}={\mathbf{f}}_{\mathrm{unif}}+{\mathbf{f}}_{\mathrm{vg}}+{\mathbf{f}}_{\mathrm{pg}}+{\mathbf{f}}_{\mathrm{rot}},$$

$${\mathrm{Re}}_p=\frac{|{\mathbf{U}}_r|D}{\nu },$$

$${\mathbf{f}}_{\mathrm{unif}}(\mathbf{x})=\left\{-\frac{F_{\mathrm{drag}}}{{\rho }_c}\xi +3\chi \xi (1-\xi )\frac{\mathbf{x}-{\mathbf{x}}_p}{|\mathbf{x}-{\mathbf{x}}_p|}\cdot \mathbf{m}\right\}\mathbf{m}-\chi \xi (1-\xi )\frac{\mathbf{x}-{\mathbf{x}}_p}{|\mathbf{x}-{\mathbf{x}}_p|}.$$

$$F_{\mathrm{drag}}$$

$$\chi$$

$${\mathbf{f}}_{\mathrm{vg}\gamma }(\mathbf{x})=-\pi \nu \gamma D^2\xi (1-\xi )\left(4\frac{x_2-x_{p2}}{|\mathbf{x}-{\mathbf{x}}_p|}{\mathbf{e}}_1+\frac{x_1-x_{p1}}{|\mathbf{x}-{\mathbf{x}}_p|}{\mathbf{e}}_2\right),$$

$${\mathbf{f}}_{\mathrm{vg}}(\mathbf{x})=-\pi \nu D^2\xi (1-\xi )\left\{4\left(\frac{\mathbf{x}-{\mathbf{x}}_p}{|\mathbf{x}-{\mathbf{x}}_p|}\cdot \nabla \right){\mathbf{U}}_{\mathrm{ud}}+(\nabla {\mathbf{U}}_{\mathrm{ud}})\cdot \frac{\mathbf{x}-{\mathbf{x}}_p}{|\mathbf{x}-{\mathbf{x}}_p|}\right\}.$$

$${\mathbf{F}}_{\mathrm{pg}}=-\frac{3m_c}{2{\rho }_c}\nabla P_{\mathrm{ud}}-\frac{m_c}{2}\frac{d{\mathbf{v}}_p}{dt},$$

$${\mathbf{f}}_{\mathrm{pg}}(\mathbf{x})=-\xi \left\{1-(1-\xi )(1-2\xi )\right\}\frac{{\mathbf{F}}_{\mathrm{pg}}}{{\rho }_c}-3\xi (1-\xi )(1-2\xi )\left(\frac{{\mathbf{F}}_{\mathrm{pg}}}{{\rho }_c}\cdot \frac{\mathbf{x}-{\mathbf{x}}_p}{|\mathbf{x}-{\mathbf{x}}_p|}\right)\frac{\mathbf{x}-{\mathbf{x}}_p}{|\mathbf{x}-{\mathbf{x}}_p|}.$$

$${\mathbf{f}}_{\mathrm{rot}}(\mathbf{x})=3\pi \nu D^2\xi (1-\xi ){\mathbf{\Omega }}_p\times \frac{\mathbf{x}-{\mathbf{x}}_p}{|\mathbf{x}-{\mathbf{x}}_p|}.$$

$${\mathrm{Re}}_p=4.64{\left(\frac{{\alpha }_c|{⟨\mathbf{u}⟩}_c-{\mathbf{v}}_p|D}{\nu }\right)}^{0.81},$$

$$\delta (f,l{\mathbf{e}}_i)=\frac{f({\mathbf{x}}_p+l{\mathbf{e}}_i)-f({\mathbf{x}}_p-l{\mathbf{e}}_i)}{2l},$$

$$\frac{\partial {\mathbf{U}}_{\mathrm{ud}}}{\partial x_i}\approx \delta ({\alpha }_c{⟨\mathbf{u}⟩}_c+{\alpha }_d{⟨\mathbf{v}⟩}_d,l{\mathbf{e}}_i),\frac{\partial P_{\mathrm{ud}}}{\partial x_i}\approx \delta (P,l{\mathbf{e}}_i).$$

$${\delta }_{urr}({\mathrm{Re}}_p)$$

$${\delta }_{pr}({\mathrm{Re}}_p)$$

$$\nabla {\mathbf{U}}_{\mathrm{ud}}$$

$$\nabla P_{\mathrm{ud}}$$

$$\mathbf{F}=F_{\mathrm{drag}}\mathbf{m}+m_c\frac{D{\mathbf{U}}_{\mathrm{ud}}}{Dt}+\frac{m_c}{2}\left(\frac{D{\mathbf{U}}_{\mathrm{ud}}}{Dt}-\frac{d{\mathbf{v}}_p}{dt}\right)+{\mathbf{F}}_{\mathrm{h}}+(m_d-m_c)\mathbf{g}+{\mathbf{F}}_{\mathrm{ext}}$$

$$\frac{d{\mathbf{v}}_p}{dt}=\frac{2}{2m_d+m_c}\left\{F_{\mathrm{drag}}\mathbf{m}-\frac{\pi D^3}{4}\nabla P_{\mathrm{ud}}+(m_d-m_c)\mathbf{g}+{\rho }_c{\int }_{V_p}\mathbf{h}dV\right\}.$$

$$\frac{d{\mathbf{v}}_p}{dt}=\frac{1}{m_d}\left\{F_{\mathrm{drag}}\mathbf{m}+(m_d-m_c)\mathbf{g}+{\rho }_c{\int }_{V_p}\mathbf{h}dV\right\}.$$

$$\frac{d{\mathbf{x}}_p}{dt}={\mathbf{v}}_p.$$

$$\frac{d{\mathbf{\Omega }}_p}{dt}=\frac{\pi {\rho }_c\nu D^3}{I_d}\left(\frac{1}{2}\nabla \times {\mathbf{U}}_{\mathrm{ud}}-{\mathbf{\Omega }}_p\right),$$

$${\alpha }_d(\mathbf{x})=\{\begin{array}{l}\frac{r_d^3}{R^3}\mathrm{for}0\le |\mathbf{x}-{\mathbf{x}}_p|<R-r_d\\ \frac{1}{16R^3}\left\{|\mathbf{x}-{\mathbf{x}}_p{|}^3-6(r_d^2+R^2)|\mathbf{x}-{\mathbf{x}}_p|-\frac{3(r_d^2-R^2{)}^2}{|\mathbf{x}-{\mathbf{x}}_p|}+8(R^3+r_d^3)\right\}\\ \mathrm{for}R-r_d\le |\mathbf{x}-{\mathbf{x}}_p|\le R+r_d\\ 0\mathrm{for}|\mathbf{x}-{\mathbf{x}}_p|>R+r_d\end{array},$$

Full text

Keywords:

particle-laden flow, volume-averaged equation, particle-vortex interaction, history effect

1 Introduction

Interactions between particles and a turbulent flow are important in many industrial processes like cyclone separators and pulverised coal combustion. Many factors determine the fluid-particle interaction such as the flow configuration, the particle relaxation time (the Stokes number), the role of flow inertia (Reynolds number), the importance of gravity (Froude number), the solid volume fraction and the mass fraction (the latter two related by the density ratio). One critical factor is the ratio between the particle size and a typical length scale of the flow. In pipe flows and free jets laden with particles (including bubbles and droplets), for example, the turbulence intensity increases when the particle diameter $D$ is larger than one-tenth of the integral length scale (Gore and Crowe, 1989). Experimental works with dilute suspensions, on the other hand, report significant reductions of the turbulence intensity when the particle diameter is comparable to the turbulence Kolmogorov length scale $\eta$ (Kulick et al., 1994; Paris and Eaton, 2001; Hwang and Eaton, 2006). To understand the mechanisms of the interaction between the phases, numerical simulations can be used to capture both the turbulence structures and the particle motion. In many numerical studies, however, the force on the particle is approximated and the feedback force on the fluid is either ignored (i.e., one-way coupling) or simplified (i.e., traditional two-way coupling) to a point-source. Thus, the turbulence attenuation by particles of $D\sim\eta$ is not reasonably reproduced by traditional two-way coupling simulations (Eaton, 2009; Schneiders et al., 2017). On the other hand, fully-resolved simulations like those in Kempe and Fröhlich (2014), Picano et al. (2015), Fornari et al. (2016) and Santarelli and Fröhlich (2016) are still too expensive for configurations of practical interest, which justifies the need for better models.

Focusing on the particle motion, the mean particle settling velocity is influenced by the background turbulence (Nielsen, 1993). As recently shown in the experimental study of Good et al. (2014), the settling velocity of particles can both increase and decrease when the particle diameter is slightly smaller than the Kolmogorov length scale. Although one-way coupling simulations capture this trend qualitatively, a quantitative difference is recognised even for very dilute cases. For direct numerical simulations of turbulent flows, a grid width $\Delta x$ of the order of the Kolmogorov length scale is necessary. Therefore, for a simulation of a particle-laden flow of particle size $D\sim\eta$, an appropriate two-way interaction model between the flow and the particle is required.

The most reliable numerical approach is resolving the particle boundary, in which case the fluid force is directly computed. The immersed boundary method (IBM) is one of the possible approaches of this type as shown in several studies (Kajishima and Takiguchi, 2002; Lucci et al., 2010; Tenneti and Subramaniam., 2014; Fornari et al., 2016). As the particle diameter needs to be resolved by ten or more grid points, fully-resolved simulations are practical when the particle is sufficiently larger than the Kolmogorov length scale. In other words, for the case $D\sim\eta$, fully-resolved simulations are not feasible because of the prohibitive computational costs. Therefore, the effect of the particle on the fluid has to be modelled without capturing the boundary layer. In the traditional two-way coupling simulations, the drag force model is based on the undisturbed flow velocity and the particle is assumed to be much smaller than the Kolmogorov length and the grid width. In the implementation, the local flow disturbance around the particle (see Fig. 1) is neglected and the disturbed velocity interpolated at the particle position is used as the undisturbed velocity in the expressions for the force (Squires and Eaton, 1990; Boivin et al., 1998; Sundaram and Collins, 1999; Li et al., 2001; Rani et al., 2004). However, as Gualtieri et al. (2015) pointed out, the effect of the disturbance around the particle itself cannot be ignored even when $D\ll\eta$. These authors proposed, therefore, an estimation of the fluid force based on the Stokes flow around the particle, still considered to be smaller than the grid size. In the case of $D\sim\eta$, the disturbance flow around the particle becomes more important and the assumption of Stokes flow is questionable. On the fluid side, moreover, the point-source feedback force on the momentum equation is numerically distributed in space. Since the particle size is ignored, the distribution does not consider the effect of the physical surface position.

To overcome these limitations, one possibility is volume averaging of the momentum equation that enable us to distribute physically-meaningful feedback force. This force is referred to as interaction force in this paper. Fukada et al. (2016) recently developed a distribution model of the interaction force for particle of diameters slightly larger than the grid size, $\Delta x$. The interaction forces for uniform and simple shear flows around a sphere are modelled for particle Reynolds numbers ${\rm Re}_{p}=\mathcal{O}(10)$ and shear Reynolds number based on the particle diameter ${\rm Re}_{\gamma}=\mathcal{O}(1)$. The asymmetric distribution of the interaction force resulted in qualitatively and quantitatively reasonable flow fields consistent with the fully-resolved results. The energy transfer on the volume-averaged field was also captured, something which is not considered in traditional two-way coupling models. However, the simulations in this previous work were limited to the case of a fixed particle and known steady undisturbed flow.

In the present study, we therefore propose a novel estimation method of the fluid force based on the disturbed flow around the particle. This approach is suited for the volume-averaged framework unlike a conventional two-way coupling approach. The study aims to show the applicability of the volume-average framework for the flow including moving particles.

We will initially consider the history effect on the particle motion, an effect whose importance is increasingly recognised (Olivieri et al., 2014; Daitche, 2015). The history effects are highly influenced by the background flow and the modelling is therefore difficult (Bagchi and Balachandar, 2003). The traditional Basset history model based on the assumption of Stokes flow (Maxey and Riley, 1983) is not applicable for a long physical time since the model overestimates the past effects (Mei and Adrian, 1992). Some models developed for finite Reynolds number are, on the other hand, limited to specific and relatively simple flows (Mei and Adrian, 1992; Wakaba and Balachandar, 2005). The high computational cost of the integration of the history effect is also a factor to consider. However, in an appropriate two-way coupling simulation, the history effects are included in the force estimation if the effect of unsteady disturbances is correctly captured (Gualtieri et al., 2015). In a similar way, the lift force can be also represented by an appropriate two-way coupling algorithm that captures the flow disturbance, again reducing the dependence on a specific model. To investigate how the history and the lift forces appear in the present simulation framework based on the volume-averaged equation (referred to as VA simulation), the settling of a particle in a fluid at rest and the flow around a rotating particle will be examined. For comparison and validation, we refer to the result from a fully-resolved IBM simulation, which is carried out in this study, and previous results (Rubinow and Keller, 1961; Kurose and Komori, 1999, Bagchi and Balachandar, 2002; Bluemink et al., 2010).

We will then focus on the interaction between particles and a cellular vortical flow, the Taylor-Green vortex. The particle diameter is $\mathcal{O}(10)$ times smaller than the vortex and its size is therefore non-negligible. The particle trajectory in the Taylor-Green vortex has been first investigated by Maxey (1987) in the one-way coupling regime. For particles of non-negligible size, however, the particle-vortex interaction leads to the flow disturbance at scales larger than the particle size as well as local disturbances around the particle. Bergougnoux et al. (2014) showed in their experimental study that weak disturbances of the vortex influence the particle trajectory significantly. Therefore, a two-way coupling investigation is necessary to correctly capture the particle motion in vortices. Flow disturbances also induce and modify the interactions between two different particles. The hydrodynamic forces on particles at distances of the order of $D$ in a uniform flow have been the main objectives of previous studies (Tsuji et al., 2003; Yoon and Yang, 2007; Ozgoren, 2013). However, for a better understanding of the particle motion in turbulence, inter-particle interactions at larger distances in a vortical flow should also be considered. In this study, we simulate the behaviour of the particle in the Taylor-Green vortices and show advantage of the VA approach for the particle motion in comparison to the point-model based method. The fully-resolved IBM simulations are also carried out as reference. Finally, the importance of two-way coupling on the particle trajectories is confirmed for cases with different initial positions and inter-particle distance of around $10D$.

2 Governing equations

2.1 Volume-averaged equations of the fluid phase

The volume-averaged mass and momentum equations for dispersed multiphase flows are derived by Anderson and Jackson (1967) under averaging length scales much larger than the inter-particle spacing. The derivation is also detailed in Crowe et al. (1997). On the other hand, the treatments of the residual stress and the interaction force terms for a case of averaging length scale comparable to the particle size have been developed by Fukada et al. (2016). The volume-averaged equations and the models of these terms are briefly described in the following.

The multiphase flow consists of the continuous (fluid) phase ($c$-phase) and the dispersed phase ($d$-phase). Only rigid spherical particles are considered for the dispersed phase. Considering the spherical averaging volume $V$ as shown in Fig. 2, the volume fraction and the phase average of a physical property $B$ are defined as

[TABLE]

where $V_{k}$ is the volume occupied by $k$-phase (either $c$ or $d$) inside $V$. For a quantity $f$ defined through the averaging volume $V$, we use the notation $f(\bm{x})$ where $\bm{x}$ denotes the centre of $V$.

The basic mass and momentum equations for an incompressible Newtonian fluid are written as follows:

[TABLE]

[TABLE]

where $\bm{u}$ is the velocity, $t$ time, $\rho_{c}$ the fluid density, $p$ pressure, $\nu$ viscosity and $\bm{h}$ an external forcing. In the present study, $\bm{h}$ is given to keep the background flow steady. As shown in Fukada et al. (2016), through volume-averaging these equations, we obtain

[TABLE]

[TABLE]

where $\bm{v}$ is the velocity inside a particle, $P$ the scalar function corresponding to pressure, $\bm{\tau}$ the residual stress and $\bm{f}$ the interaction force. The fluid variables used in the simulations are $\alpha_{c}\left<\bm{u}\right>_{c}$ and $P$. Appendix A shows the differentiability of the volume-averaged quantities. The form of the viscous term is different from that in other volume-averaged equations (Anderson and Jackson, 1967; Wachem et al., 2001), and the decomposition into $\nu\nabla^{2}(\alpha_{c}\left<\bm{u}\right>_{c})+\nu\nabla^{2}(\alpha_{d}\left<\bm{v}\right>_{d})$ is not allowed.

The scalar function $P$ can be decomposed as $\alpha_{c}\left<p\right>_{c}+\sum_{n}\alpha_{d,n}\hat{p}_{n}$ where $\alpha_{d,n}$ and $\hat{p}_{n}$ are the volume fraction and the surface-mean pressure (over the entire surface) of the $n$th particle. In a numerical simulation, $P$ can be obtained without considering the decomposition.

The residual stress is defined as

[TABLE]

where $\delta\bm{u}=\bm{u}-\left<\bm{u}\right>_{c}$. The original model by Fukada et al. (2016) is

[TABLE]

where $x_{m}$ is the $m$th component of the Cartesian coordinates and the summation convention is applied for the subscript $m$. The velocity gradient inside the rigid particle $\partial\bm{v}/\partial x_{m}$ corresponds to the angular velocity of the particle.

The interaction force is defined as

[TABLE]

where $\delta p$ is the deviation from the surface-mean pressure $\hat{p}_{n}$ of the corresponding particle, $\bm{n}$ the unit normal vector on the particle surface directed to the fluid-phase and $S_{d}$ the partitioned particle surface area inside $V$ (see Fig. 2).

2.2 Interaction force models

According to Appendix B of our previous study (Fukada et al., 2016), the independent interaction force models can be successfully superimposed for steady flows with ${\rm Re}_{p}\leq 40$ and $\gamma D^{2}/\nu=\mathcal{O}(1)$, where $\gamma$ is the shear rate. In the present study, we thus assume superposition of the different contributions:

[TABLE]

where the contributions due to the relative velocity, the undisturbed velocity gradient, the undisturbed pressure gradient and the particle rotation are represented as $\bm{f}_{\rm unif}$, $\bm{f}_{\rm vg}$, $\bm{f}_{\rm pg}$ and $\bm{f}_{\rm rot}$. Each component is described in the following.

Introducing the particle Reynolds number as

[TABLE]

where $\bm{U}_{r}$ is the relative velocity based on the undisturbed flow and $D$ is the particle diameter, the interaction force corresponding to the relative velocity is modelled as (Fukada et al., 2016)

[TABLE]

In the above expression, $\xi=S_{d}/\pi D^{2}$ denotes the normalised surface area (hereafter, referred to as surface fraction), $\bm{m}$ the unit vector in the direction of the relative velocity $\bm{U}_{r}$, $\bm{x}_{p}$ the particle centre position, $F_{\rm drag}$ the contribution of the relative velocity to the fluid force

[TABLE]

and $\chi$ is a fitting coefficient

[TABLE]

Note that $V^{-1}\int\bm{f}_{\rm unif}d\bm{x}=-\bm{F_{\rm drag}}/\rho_{c}$ holds, which guarantees the total momentum conservation between the phases. The interaction force (12) is based on the theoretical result of the flow around the sphere for low ${\rm Re}_{p}$ (Proudman and Pearson, 1965). The fitting coefficient $\chi$ is introduced for extending to higher ${\rm Re}_{p}$ according to our numerical result (Fukada et al., 2016).

For a uniform shear flow of undisturbed velocity $\bm{U}_{\rm ud}=\gamma(x_{2}-x_{p2})\bm{e}_{1}$ with ${\rm Re}_{p}=0$, the shear-induced interaction force is modelled as (Fukada et al., 2016)

[TABLE]

where $\bm{e}_{i}$ is the unit basis vector in the $i$th direction. The model for a general undisturbed velocity gradient can be written as the superposition

[TABLE]

Appendix B of this paper briefly summarises the derivation, and for more details the reader is referred to Fukada et al. (2016).

The pressure gradient and the added-mass forces on the particle, $\bm{F}_{\rm pg}$, are modelled as

[TABLE]

where $P_{\rm ud}$ is the undisturbed pressure and $\bm{v}_{p}$ the translating velocity of the particle, $m_{c}=\pi\rho_{c}D^{3}/6$ the mass of the displaced fluid. Assuming that this force corresponds to the additional surface pressure $-(3\bm{F}_{\rm pg}\cdot\bm{n}/\pi D^{2})$ obtained for an inviscid uniform flow, the interaction force can be shown to be (see Appendix B)

[TABLE]

Using the Stokes solution for a rotating sphere with angular velocity $\bm{\Omega}_{p}$, the interaction force due to the particle rotation is obtained as

[TABLE]

2.3 Estimation of undisturbed flow at the particle position

To represent the interaction between the fluid and the particles, the relative velocity $\bm{U}_{r}$ and the undisturbed gradients of $\bm{U}_{\rm ud}$ and $P_{\rm ud}$ at the particle position need to be estimated. In traditional two-way coupling simulations, the disturbed fluid velocity interpolated at the particle position is often regarded as the undisturbed flow when computing the interaction. However, this treatment is justified only when the particle is much smaller than the grid spacing and not appropriate for the present cases ($D\sim\Delta x$). Therefore, in this paper, we propose new estimation methods for the relative velocity $\bm{U}_{r}$ and the undisturbed gradients of $\bm{U}_{\rm ud}$ and $P_{\rm ud}$. The radius $R$ of the averaging volume $V$ is kept to be $0.75D$ throughout this paper as the effectiveness of this value was confirmed (Fukada et al., 2016).

2.3.1 Relative velocity

The contribution of the relative velocity to the fluid force (13) is described in terms of ${\rm Re}_{p}$. As the volume-averaged disturbed velocity $\alpha_{c}\left<\bm{u}\right>_{c}$ is obtained from the volume-averaged equations, the correlation for the particle Reynolds number ${\rm Re}_{p}$ is necessarily based on $(\alpha_{c}\left<\bm{u}\right>_{c})(\bm{x}_{p})$. For the fixed value $R=0.75D$, the correlation equation,

[TABLE]

is obtained by a curve fitting based on the numerical data of our previous study (Fukada et al., 2016) for the steady uniform flow around a single particle. Figure 3 shows that Eq. (20) is valid for ${\rm Re}_{p}\leq 40$, which is sufficient for the present study. The effect of unsteady velocity on the force is qualitatively reflected through $|\left<\bm{u}\right>_{c}-\bm{v}_{p}|$. The estimation method for $\alpha_{c}(\left<\bm{u}\right>_{c}-\bm{v}_{p})$ is in Sec. 3.1. In the numerical implementation, the direction of the relative velocity is assumed to be the same as that of $\alpha_{c}(\left<\bm{u}\right>_{c}-\bm{v}_{p})$.

2.3.2 Velocity and pressure gradients

A simple estimation method for the undisturbed gradients based on the volume-averaged variables is proposed. We define the differentiation operator in the direction $\bm{e}_{i}$ at the particle centre:

[TABLE]

where $f$ is a function of $\bm{x}$, and $l$ is an appropriately determined distance. In this study, $l=D$ is adopted as supported by the tests for the pressure gradient shown in Appendix C. One of the simplest possible estimates of the undisturbed gradients are

[TABLE]

However, even in a uniform flow, non-zero gradients would be estimated around a particle due to the particle relative motion. To remove the effect of the relative motion, the following equations are constructed for the uniform flow case by curve fittings:

[TABLE]

where the direction of the relative velocity $\bm{m}$ is $(\left<\bm{u}\right>_{c}-\bm{v}_{p})/|\left<\bm{u}\right>_{c}-\bm{v}_{p}|$. Figures 4 and 5 compare Eqs. (23) and (24) with the fully-resolved numerical data obtained on a body-fit coordinate system (Fukada et al., 2016). According to these correlations, the undisturbed gradients are approximated by

[TABLE]

where $\bm{I}$ is the identity tensor and the summation convention is applied for the subscript $i$.

2.4 Equation of motion for the dispersed phase

For finite Reynolds numbers, the force $\bm{F}$ on a particle can be modelled as

[TABLE]

(see Crowe et al., 1997), where $\bm{F}_{\rm h}$ is the history force, $\bm{g}$ the gravitational acceleration, $\bm{F}_{\rm ext}$ the external force, $m_{d}=\pi\rho_{d}D^{3}/6$ the particle mass. The first three terms on the right-hand side are the steady viscous force, the pressure gradient force and the added-mass force. The term $\rho_{c}(D\bm{U}_{\rm ud}/Dt)$ is replaced by $-\nabla P_{\rm ud}$ because the viscous force $\nu\nabla^{2}\bm{U}_{\rm ud}$ cancels with the external force $\bm{h}$ for the undisturbed flows considered in the present study. The external force on the particle $\bm{F}_{\rm ext}$ is the same as that on the displaced fluid $\rho_{c}\int_{V_{p}}\bm{h}dV$ where $V_{p}$ is the volume of the particle. For the history force at a finite Reynolds number, a reliable model for general flows is not available so far. However, as discussed by Gualtieri et al. (2015), the history effect can be partly reproduced in accurate two-way coupling simulations without any specific model. Therefore, in the present model, we will numerically solve the following equation

[TABLE]

To investigate the effects of the pressure gradient and the added mass, we will also compare the result by solving the following equation

[TABLE]

The particle position is given by the following equation:

[TABLE]

To include a first-order approximation of the effect of the flow on the particle rotation, the following equation for the angular velocity based on the Stokes solution is considered

[TABLE]

where $I_{d}=m_{d}D^{2}/10$ is the moment of inertia of the spherical particle. This equation is consistent with Eqs. (16) and (19) as both equations assume the same Stokes solution.

3 Numerical methods

For simplicity, the simulation with the volume-averaged equations (Secs. 2.1 and 2.2) is referred to as VA simulation. To investigate the history effect, a one-way coupling simulation is also attempted. Moreover, we will compare the results of the VA simulation to those of a traditional two-way coupling simulation and the fully-resolved simulation with the immersed boundary method (IBM). In the following, these numerical procedures are described briefly. In all the simulations, the 2nd-order central-difference scheme is used for the spatial derivatives with a staggered arrangement for the fluid variables. The fractional step method (Kim and Moin, 1985) is used as the pressure-velocity coupling algorithm for the fluid-phase. The computational cell is a cube of side length $\Delta x$.

3.1 VA simulation

The 2nd-order Runge-Kutta method is used for the time evolution of both phases. As the averaging volume $V$ is larger than the particle ($R=0.75D$), the contribution of each particle to the volume fraction $\alpha_{d}$ is calculated as

[TABLE]

where $r_{d}=D/2$ is the radius of the particle. The surface fraction $\xi$ is calculated as

[TABLE]

The averaged velocity of the solid-phase $\alpha_{d}\left<\bm{v}\right>_{d}$ is calculated as

[TABLE]

where the function $K_{d}(y)$ is

[TABLE]

The derivation of these geometrical functions (32)–(34) is summarised in Appendix D.

To keep the total external force on the system $\int\bm{h}dV$ constant, the external forces are approximated as

[TABLE]

for fluid, in Eq. (6), and

[TABLE]

for solid, in Eq. (28), where the subscript $ijk$ represents the spatial point of grid index $(i,j,k)$.

The estimation of $\{\alpha_{c}(\left<\bm{u}\right>_{c}-\bm{v}_{p})\}(\bm{x}_{p})$ in Eq. (20) is particularly important to predict the drag force. As the distribution of $\alpha_{c}\left<\bm{u}\right>_{c}$ has a local minimum near the particle centre, a linear interpolation is not sufficient. Therefore, the following interpolation steps are used for the velocity:

[TABLE]

[TABLE]

where $l=1,\cdots,8$ indicates the eight definition points of the velocity around $\bm{x}_{p}$. The effect of the second-order derivative is considered in (39) and the linear interpolation (40) guarantees the continuity about $\bm{x}_{p}$. For the pressure, the linear interpolation

[TABLE]

is used ($l$ indicates the definition points of the pressure).

To adjust the flow field far away from the particle to the non-averaged field, the residual stress term $\nabla\cdot\bm{\tau}$ is replaced by $C(\nabla\cdot\bm{\tau})$ where $C$ is given as follows:

[TABLE]

3.2 One-way coupling simulation

In the one-way coupling simulation, the flow field ($\bm{U}_{\rm ud}$, $P_{\rm ud}$) is theoretically or numerically given and only the particle motion is solved for. To consider the history effect, the acceleration including the Basset history term is considered. The equation of motion is given by

[TABLE]

Eqs. (30) and (43) are solved with the efficient implicit method proposed by van Hinsberg et al. (2011). The linear drag force model ($F_{\rm drag}=3\pi\rho_{c}\nu D|\bm{U}_{\rm ud}-\bm{v}_{p}|$) is used as they did. Simulations are also performed without the Basset term for comparison.

To quantify the relevance of the nonlinear drag force model, we solve Eqs. (28) and (30) without the external force, $\bm{h}$, with a 2nd-order Runge-Kutta method. Note that the volume force $\bm{h}$ is entirely distributed on the fluid. The particle Reynolds number used here is given by

[TABLE]

instead of Eq. (20); the direction defined by the unit vector $\bm{m}$ is parallel to $\bm{U}_{\rm ud}-\bm{v}_{p}$.

In the present study, to examine the history effect in time when ${\rm Re}_{p}$ is low, the following three cases of one-way coupling simulations are performed and those are denoted as:

O-LB:

including the linear drag, the added mass, the pressure gradient and the Basset terms, see Eq. (43),

O-L:**  **

same as O-LB but without the Basset term,

O-NL:

including the nonlinear drag, the added mass and the pressure gradient, see Eq. (28).

3.3 Traditional two-way coupling simulation

In the traditional two-way coupling simulation, Eq. (3) and the momentum equation

[TABLE]

are solved without volume-averaging, where $\bm{f}_{b}$ is the feedback force from the particle. Eq. (29) with $\bm{h}=0$ (as $\bm{h}$ is entirely distributed on the fluid) and Eq. (30) are solved for the particle motion. This traditional two-way coupling simulation is referred to as TT below.

To distribute the feedback force $\bm{f}_{b}$, we use the following equation (regularised Dirac delta function)

[TABLE]

where $K$ is the normalisation factor computed as

[TABLE]

The direction $\bm{m}$, used in eq. (46), is determined by $\bm{m}=(\bm{u}(\bm{x}_{p})-\bm{v}_{p})/|\bm{u}(\bm{x}_{p})-\bm{v}_{p}|$.

The fluid velocity at the particle centre is interpolated with

[TABLE]

To compute the drag force by Eq. (13), the particle Reynolds number is estimated as

[TABLE]

The numerical procedure with the 2nd-order Runge-Kutta method is employed for both phases.

3.4 Fully-resolved simulation

The immersed boundary code originally developed by Breugem (2012) is used for the fully-resolved simulation. The continuity equation (3) and the following momentum equation are solved in the whole domain including the regions occupied by the particles:

[TABLE]

where $\bm{f}_{\rm IB}$ is the body force used to impose the no-slip condition on the particle surface. The particle translational and rotational equations are

[TABLE]

where $S$ is the particle surface, $\bm{s}=-p\bm{I}+\rho_{c}\nu(\nabla\bm{u}+\nabla\bm{u}^{T})$ the stress tensor. The force exchange is considered on a set of $N$ Lagrangian points around each particle surface. The force $\bm{F}_{j}$ at the $j$th Lagrangian point is distributed on the fluid as

[TABLE]

where $\delta_{d}$ is a regularised Dirac delta function and $\Delta V_{j}$ the volume of the Lagrangian grid cell. In the simulation, Eqs. (52) and (53) are converted to

[TABLE]

The three-step Runge-Kutta method is used for the time integration. More details can be found in Breugem (2012) and Lambert et al. (2013).

4 Numerical results

In the VA simulation, the history and lift forces are captured without any specific treatment as this method intrinsically incorporates the effect of the flow perturbation through the volume-averaging. In Sec 4.1, the particle settling problem in a stationary fluid is simulated to show how the history force is represented. The lift force on a rotating particle is instead the focus of Sec 4.2. For the study of the turbulence modulation by particles of comparable size to the Kolmogorov length scale, the interaction between the particle and a vortex element should be precisely represented. In Sec 4.3, therefore, the applicability of the VA simulation for the Taylor-Green vortex is investigated. As the fundamental validation for different density ratios, the particle motion in the smallest periodic unit is simulated without gravity. Finally, in Sec 4.4, to highlight the importance of the two-way coupling simulation for a vortical flow, we show the particle trajectory and the inter-particle interaction in an array of Taylor-Green units with gravity.

Throughout this study, grid resolutions of $D/\Delta x=24$ for the fully-resolved simulations and $D/\Delta x=2$ for the VA simulations are commonly employed. Note that the number of grid points is $12^{3}$ times lower for the VA simulation with respect to the fully-resolved simulation and consequently the time step $\Delta t$ is $20$ times larger. Therefore, the total computational cost is $\mathcal{O}(10^{4})$ times lower with the VA model. The computational domain is rectangular of lengths $l_{1},l_{2}$ and $l_{3}$ in the $x_{1},x_{2}$ and $x_{3}$-directions. Periodic boundary conditions are applied in all the directions. The motion of the particle is confined in the $x_{1}$-$x_{2}$ plane due to the symmetry of the flows studied. The choice of $R=0.75D$ and $D/\Delta x=2$ is reasonable according to our previous study (Fulada et al., 2016). The effects of these parameters are further discussed in Appendix E.

4.1 History effect on the settling particle

A single particle settling in a stationary fluid is studied by one-way coupling simulations, the VA approach and the IBM simulation. The fluid and particle velocities are initially set to [math]. Gravity acts in the negative $x_{2}$-direction and the external force, $\bm{h}$, is neglected. The importance of gravity is characterised by the Galileo number defined as

[TABLE]

In the following, using Eq. (13), ${\rm Ga}$ is related to the particle Reynolds number based on the particle terminal velocity as

[TABLE]

The following set of parameters are used for the simulations: $(l_{1}/D,l_{2}/D,l_{3}/D)=(16,32,16)$, $\rho_{d}/\rho_{c}=100$ and ${\rm Ga}=8.44$ (corresponding to ${\rm Re}_{\rm term}=3$). The time step is $(\nu/D^{2})\Delta t=1.19\times 10^{-4}$ for the fully-resolved simulation and $(\nu/D^{2})\Delta t=2.39\times 10^{-3}$ for the other cases. The number of grid points is $384\times 768\times 384$ for the fully-resolved simulation and $32\times 64\times 32$ for the VA simulation.

Figure 6 shows the time evolution of the dimensionless particle settling velocity $v_{p2}D/\nu$. First we note that the solid and dashed lines almost overlap with each other and the proposed VA and the fully-resolved simulations show good agreement. Around $(\nu/D^{2})t=12$, the result of the O-NL simulation, including drag, added-mass and pressure gradient, shows better agreement with that of the fully-resolved simulation due to the nonlinear drag model (13). On the other hand, the O-LB simulation, including the history effect, shows better agreement with the fully-resolved case only for the earlier stage. Therefore, the history force is essential to correctly model the initial transient stage, which is captured in the VA simulation. In Fig. 6(b), focusing on the initial stages of the particle motions, the difference between the results of the two one-way coupling simulations without the Basset term (O-L and O-NL) is small because the nonlinear effect in the drag force is not significant at the initial stage when the particle Reynolds number is low.

As the boundary layer thickness at the beginning of the settling is smaller than that in the steady flow, the friction drag of the unsteady flow becomes larger in the developing stage. In the VA simulation, smaller boundary layer thickness corresponds to larger $|\alpha_{c}(\left<\bm{u}\right>_{c}-\bm{v}_{p})|$ and the history effect is qualitatively reflected in the drag force. This also explains why the result of the VA simulation shows quantitatively good agreement with the fully-resolved simulation.

4.2 Lift force induced by particle rotation

The VA simulation of the flow around a rotating particle is carried out to test the capability of capturing the transversal forces. The uniform velocity $\bm{u}=(U_{\rm init},0,0)$ is given as the initial condition for the fluid flow. The particle centre is fixed in space and the angular velocity is kept constant to $\bm{\Omega}_{p}=(0,\Omega_{\rm const},0)$, thus the particle motion, Eqs. (28), (30) and (31), does not need to be solved. The Reynolds numbers and the angular velocities are varied in the following range: $U_{\rm init}D/\nu=1,5,10,20$ and $\Omega_{\rm const}D/U_{\rm init}=0.196,0.393$, with gravity and the external forces set to zero. The size of the computational domain is $(l_{1}/D,l_{2}/D,l_{3}/D)=(64,32,16)$ and the number of grid points is $128\times 64\times 32$. The time step is $(U_{\rm init}/D)\Delta t=2.81\times 10^{-2}$. The wake of the particle reaches the particle position around $(U_{\rm init}/D)t=64$ due to the periodic boundary condition. The force is thus examined at $(U_{\rm init}/D)t=28.1$ so that effects from the re-entering wake are avoided.

Based on the components of the estimated fluid force $\bm{F}$, the drag and lift coefficients, $C_{D}$ and $C_{L}$, are defined as follows:

[TABLE]

The drag and lift coefficients obtained in the VA simulation are plotted in Fig. 7. Note that the magnitude of the angular velocity does not influence the two coefficients. The drag coefficient estimated in the present simulation shows good agreement with that based on Eq. (12) (solid line). Therefore, the effect of the rotation on the drag force is small as supported by previous researches (Sridhar and Katz, 1995; Bagchi and Balachandar, 2002). As for the lift force, the signs of $\alpha_{c}\left<u_{2}\right>_{c}$ at the particle centre and $F_{2}$ should be the same according to the result that the contribution of the friction lift is in the same direction as $F_{2}$ (Kurose and Komori, 1999). Therefore, the present force estimation (Sec. 2.3.1), the directions of the drag force and $\alpha_{c}(\left<\bm{u}\right>_{c}-\bm{v}_{p})$ being the same, is capable of capturing the direction of the lift force. In the VA simulation, the interaction force model for the particle rotation, Eq. (19), induces $\alpha_{c}\left<u_{2}\right>_{c}<0$ and thus the lift force $F_{2}<0$. According to the theoretical study by Rubinow and Keller (1961), the lift coefficient is $C_{L}=1$ for ${\rm Re}_{p}\ll 1$. On the other hand, numerical studies at ${\rm Re}_{p}=\mathcal{O}(10)$ have shown an estimate of $C_{L}\approx 0.5$ (Bagchi and Balachandar, 2002; Bluemink et al., 2010). The VA simulation captures the direction of the lift force generated by the particle rotation, and the magnitude quantitatively agrees with the previous results for ${\rm Re}_{p}<10$ under the setting of this study.

4.3 Vortical flow without gravity for three density ratios

To study the interaction between a particle and a vortex, the Taylor-Green vortex is used as the background undisturbed flow. The smallest unit structure of the Taylor-Green vortex is considered to compare the results of relatively simple particle motions from different simulations. The directions of the Cartesian coordinates ($x_{1},x_{2},x_{3}$) are determined so that the velocity components of the undisturbed flow are

[TABLE]

where the velocity $A$ and the length $L$ define the vortex intensity and size. The period in the $x_{1}$ and $x_{2}$-directions is $2\pi L$ and the Reynolds number ${\rm Re}=AL/\nu=18$. According to Jiménetz et al. (1993), the intensity of a typical vortex in isotropic turbulence is correlated as $\Gamma/\nu\approx 18\sqrt{{\rm Re}_{\lambda}}$, where $\Gamma$ is the circulation of the vortex and ${\rm Re}_{\lambda}$ the Reynolds number based on the Taylor length scale. The present case (${\rm Re}=18$), where the circulation of one vortex is $\Gamma=16AL$, corresponds therefore to ${\rm Re}_{\lambda}\approx 250$. The size of the computational domain is $l_{1}/L=l_{2}/L=l_{3}/L=2\pi$ and the particle diameter $D/L=2\pi/16$. As the Kolmogorov length scale $\eta$ is around eight times smaller than the diameter of the most intense vortices in turbulence (Jiménetz et al., 1993), the present particle diameter is considered as a model of the case $D\sim\eta$. The flow is maintained by the external force

[TABLE]

The number of grid points is $384\times 384\times 384$ for the fully-resolved simulation and $32\times 32\times 32$ for the VA and the TT simulations. The time step is $(A/L)\Delta t=3.31\times 10^{-4}$ for the fully-resolved simulation and $(A/L)\Delta t=6.63\times 10^{-3}$ for the other two methods. Three different density ratios ($\rho_{d}/\rho_{c}=1,10,1000$) and two different initial particle positions are examined in the following. The corresponding Stokes numbers ${\rm St}=\rho_{d}D^{2}A/(18\rho_{c}\nu L)$ are 0.154, 1.54 and 154. The initial velocity of the particle is set to be the same as the undisturbed fluid velocity at the particle centre and the initial angular velocity is 0.

4.3.1 $\rho_{d}/\rho_{c}=1000$

When the initial particle position is $(x_{1}/L,x_{2}/L,x_{3}/L)=(\pi/2,\pi/2,0)$, the particle trajectory follows the straight line defined by $dx_{2}/dx_{1}=-1$ through the periodic boundaries and the particle does not rotate. To highlight the difference between the VA simulation and the TT simulation, the VA simulation is repeated without considering the pressure gradient, particle rotation and the external force on the particle. This simplified VA simulation is referred to as SVA simulation.

The time evolution of the particle velocity $v_{p1}$ from these different simulations are compared in Fig. 8. The results of the VA and SVA simulations are similar to that of the fully-resolved IBM simulation, which we take as the reference case. On the other hand, the particle behaviour predicted by the TT simulation exhibits large difference from the reference case. One of the most significant differences between the SVA and the TT simulations is the estimation of the drag force. As the flow disturbance is non-negligible for the finite-size particle, the force estimation according to Eq. (50) without considering the local flow disturbance underestimates the drag force and results in the smaller acceleration of the particle in the TT simulation.

To investigate the effect of the particle on the vortex, we define the induced flow disturbance as $(\alpha_{c}\left<\bm{u}\right>_{c}+\alpha_{d}\left<\bm{v}\right>_{d}-\bm{U}_{\rm ud})$. For comparison, the induced flow disturbance for the fully-resolved simulation is defined as $(\alpha_{c}\left<\bm{u}\right>_{c}+\alpha_{d}\left<\bm{v}\right>_{d}-\bm{U}_{\rm ud})$ using the local velocities only in the region where $\alpha_{c}>0$, while $(\bm{u}-\bm{U}_{\rm ud})$ is used in the other region. Figure 9(a) shows the induced flow disturbances at time $(A/L)t=33.13$ in the $x_{1}$-$x_{2}$ cross-section cutting through the particle for both the VA simulation (solid arrow) and the fully-resolved simulation (dashed arrow). This figure indicates that the disturbances at larger scales than the particle size are very close to each other. Relatively larger differences are found in the area closer to the particle due to the difference in the position of the particle. As shown in Fig. 9(b), by extracting the data from the VA simulation at time $(A/L)t=33.36$, to match the particle position to that of the fully resolved case, the difference in the flow disturbance becomes smaller. To summarise, the VA simulation shows a better agreement with the fully-resolved results for both the flow disturbance and the particle motion in comparison to the one-way coupling and the TT models.

4.3.2 $\rho_{d}/\rho_{c}=10$

The initial particle position is given as $(x_{1}/L,x_{2}/L,x_{3}/L)=(\pi/2,\pi,0)$ so that the particle trajectory bends due to the vortical flow. The particle trajectories for the different simulation methods are compared in Fig. 10. The result of the VA simulation is very similar to that of the reference fully-resolved IBM simulation. The effects of the particle rotation, pressure gradient and external force are not significant as the result of the SVA simulation is also very close to the two previous cases. As discussed in Sec. 4.3.1, the drag force estimated in the TT simulation is smaller than that in the VA simulation, which gives smaller acceleration in the $x_{1}$-direction at the early stage. In the present case where the effect of the pressure gradient is not so significant, the VA simulation effectively reproduces the curved particle trajectory with significantly less spatial resolution.

Figure 11 shows the time evolution of the angular velocity $\Omega_{p3}$. The result of the VA simulation shows good agreement with that of the fully-resolved simulation. Therefore, the contribution of the vorticity to $\Omega_{p3}$ is reasonably reproduced by the proposed model.

Finally, Fig. 12 (a) shows the induced disturbance velocity field at $(A/L)t=33.13$ in the $x_{1}$-$x_{2}$ cross-section cutting through the particle for the VA and the fully-resolved IBM simulations. The disturbances at larger scales than the particle size show good agreement with each other. As shown in Fig. 12 (b), the disturbances around the particles give an even better agreement when the particle position of the VA simulation is adjusted to that of the fully-resolved simulation by slightly changing the time ($(A/L)t=32.53$).

4.3.3 $\rho_{d}/\rho_{c}=1$

We next shortly consider particles of density equal to that of the fluid with initial particle position $(x_{1}/L,x_{2}/L,x_{3}/L)=(\pi/2,\pi,0)$. The particle trajectories are compared in Fig. 13. As the density ratio is 1, the particle velocity fluctuations are relatively large. The result of the O-NL simulation including all the forces except for the history effect and external force, Eq. (28), shows good agreement with that of the IBM simulation; the streamlines are almost closed. On the other hand, the O-NL simulation further neglecting the pressure gradient and added mass forces, Eq. (29), shows a totally different trend, suggesting that the pressure gradient gives an important contribution. The result of the VA simulation is also different from that of the fully-resolved IBM simulation. Therefore, the estimation of the fluid force needs to be improved for the case where the pressure force is dominant and the particle velocity fluctuations are large.

4.4 Effects of two-way coupling in vortical flow with gravity

The settling motion of a particle is investigated in an array of the Taylor-Green units. The flow configuration is

[TABLE]

as often used (see e.g. Maxey, 1987; Bergougnoux et al., 2014). Note that $A$ and $L$ in Eq. (62) are $\sqrt{2}$ times larger than those in Eq. (60). The length of the unit cell is $2\pi L$ in the $x_{1}$ and $x_{2}$-directions and the flow is maintained by the external force:

[TABLE]

As in the previous section, we compare particle trajectories and the induced disturbance velocity field obtained by different numerical models, with particular emphasis on investigating the history effect on the trajectory for different initial particle positions. Finally, the flow-mediated interaction between multiple particles at distances around $10D$ is studied. All the results in this section are obtained for density ratio $\rho_{d}/\rho_{c}=100$ and the Reynolds number ${\rm Re}=AL/\nu=30$ (corresponding to ${\rm Re}=15$ with the definition in Sec. 4.3). The time step is $(A/L)\Delta t=5.52\times 10^{-4}$ for the fully-resolved simulation and $(A/L)\Delta t=1.10\times 10^{-2}$ in the other cases. Gravity works in the negative $x_{2}$-direction. For all the simulations, the initial particle velocity is the same as the flow at the particle position and the angular velocity is 0.

4.4.1 Validation of the VA simulation

The particle initial position is $(x_{1}/L,x_{2}/L,x_{3}/L)=(\pi/2,\pi/2,0)$ where the flow velocity is 0. The particle diameter is $D/L=2\pi/16$ (${\rm St}=12.9$) and the Galileo number is ${\rm Ga}=8.44$. The domain size is $l_{1}/L=4\pi$ and $l_{2}/L=l_{3}/L=2\pi$. The number of grid points is $768\times 384\times 384$ for the fully-resolved IBM simulation and $64\times 32\times 32$ for the VA and the TT simulations.

Figure 14 (a) compares the particle trajectories for the different simulations. The result of the VA simulation shows good agreement with that of the fully-resolved IBM simulation. The trajectory is not as simple as in the no-gravity cases: the particle is accelerated by gravity initially and then is transported upward by the vortex. Interestingly, the TT simulation does not yield the upward particle motion due to the reduced value of the drag force. Figure 14 (b) compares the initial stage of the trajectories obtained by the different formulations including the one-way coupling regime. The result of the O-LB simulation, Eq. (43), is the closest (among the one-way coupling simulations) to that of the fully-resolved IBM simulation. Therefore, in the initial stage, the history effect is more important than the nonlinearity of the drag model. However, the difference in the trajectory of the O-LB simulation increases with respect to the reference case after the initial stage of Fig. 14 (b) (i.e., when the particle goes into the neighbouring vortex) because of the linear drag model and the error in the Basset term at longer times. The result of the VA simulation suggests that an appropriate two-way coupling model can reproduce the history effect without a complicated model when the effect of unsteady disturbance due to the finite-size particle is reflected.

The particle Reynolds number exhibits temporal variation up to around 10 with accelerated and decelerated motion (figure omitted).

4.4.2 History effect on particle trajectories

To highlight the importance to include the history effect, we investigate the particle trajectories for different initial particle positions and two particle diameters, $D/L=2\pi/16$ (${\rm Ga}=8.44$, ${\rm St}=12.8$) and $D/L=2\pi/32$ (${\rm Ga}=2.98$, ${\rm St}=3.21$). The domain size is $l_{1}/L=l_{2}/L=8\pi$ and $l_{3}/L=4\pi$ for the larger particle and $l_{1}/L=l_{2}/L=4\pi$ and $l_{3}/L=2\pi$ for the smaller particle. The number of grid points is $128\times 128\times 64$ for both cases.

The trajectories pertaining five different initial particle positions (along an enclosed streamline and at the vortex centre), $(x_{1}/L,x_{2}/L)=(\pi/2,\pi/2)$, $(\pi/4,\pi/2)$, $(3\pi/4,\pi/2)$, $(\pi/2,\pi/4)$ and $(\pi/2,3\pi/4)$, are displayed in Fig. 15. The trajectories are obtained with the VA approach and the O-NL simulation excluding the history effect. The trajectories obtained with the VA simulation have slightly larger curvature than those from the O-NL simulation at the early stage, which is consistent with the observations above about the role of the history effects. The distances between the corresponding trajectories increase with time. For the larger particle (Fig. 15 (a)), the differences are already non-negligible in the cell adjacent to that of the initial particle positions. For the smaller particle (Fig. 15 (b)), except for the particle with the initial position $(x_{1}/L,x_{2}/L)=(\pi/2,\pi/2)$, the differences between the two models are relatively small. This trend is explained by the fact that the history effect becomes smaller for smaller particles (Bergougnoux et al., 2014; Daitche, 2015). For the case with the initial position $(x_{1}/L,x_{2}/L)=(\pi/2,\pi/2)$, the long-time less-active motion around the vortex centre enhances the history effect on the trajectory.

4.4.3 Interaction between multiple particles

The importance of the two-way coupling simulation for the inter-particle interaction through the flow disturbance is demonstrated in the following. For the simulations presented here, the physical parameters are $D/L=2\pi/16$ and ${\rm Ga}=8.44$. The domain size is $l_{1}/L=l_{2}/L=8\pi$ and $l_{3}/L=4\pi$ and the number of grid points $128\times 128\times 64$. The interaction between particles at distances of around $10D$, which is a typical distance for volume fraction $\mathcal{O}(10^{-4})$, is simulated with 3 particles with initial positions $(x_{1}/L,x_{2}/L,x_{3}/L)=(\pi/2,5\pi/2,0)$, $(5\pi/2,5\pi/2,0)$ and $(3\pi/2,3\pi/2,0)$. Note that the first two particles are in the same relative position of the respective Taylor-Green vortex units. Figure 16 shows the disturbance flow field and the particle trajectories at two different instants. If the inter-particle interactions are ignored, the trajectories of the two particles initially at $x_{2}/L=5\pi/2$ should be the same. As shown in Fig. 16 (b) at time $(A/L)t=33.12$, however, the trajectory of the particle released from $(x_{1}/L,x_{2}/L,x_{3}/L)=(\pi/2,5\pi/2,0)$ turns to a different direction in comparison to that released from $(x_{1}/L,x_{2}/L,x_{3}/L)=(5\pi/2,5\pi/2,0)$. The present result suggests that the particle motion for $D\sim\eta$ is clearly influenced by other particles at distance around $10D$. Also, the flow disturbance around $(x_{1}/L,x_{2}/L,x_{3}/L)=(\pi/2,5\pi/2,0)$ is larger than that around $(x_{1}/L,x_{2}/L,x_{3}/L)=(5\pi/2,5\pi/2,0)$ owing to the inter-particle interaction. The spreading of the disturbance velocity over a wide region is caused by convection since the convective time scale $(L/A)$ is sufficiently smaller than the viscous time scale $(L^{2}/\nu)$ (i.e., $\nu/AL=1/30$) in our case. As the modeling of the convective effect is difficult, the two-way coupling simulation is necessary to investigate the inter-particle interaction.

5 Conclusion

For the simulation of flows laden with particles of size comparable to the smallest turbulent eddies, $D\sim\eta$, we have previously developed an interaction force model based on the volume-averaged continuity and momentum equations. In this paper, we proposed a new method to estimate the fluid force to enable simulation of the transport of particles within the same volume-averaged framework (VA simulation). The VA velocity at the particle centre is correlated with the particle Reynolds number. At the same time, the effects of the pressure gradient, the velocity gradients and the particle rotation are incorporated into the interaction force model. The qualitative advantages of the VA approach are the capability of representing the history effect without a complicated model and the better drag estimation compared to the traditional point-model based method.

To test the proposed model, we set up configurations of increasing complexity and compared the results with those obtained with interface-resolved simulations based on the immersed-boundary method (IBM). When considering a single settling particle in a stationary fluid, we showed that the history effect was captured in the VA simulation without any specific model. We then examined the flow around a rotating particle at ${\rm Re}_{p}\leq 20$ and showed that the direction of the lift force was represented by the model, and the magnitudes for ${\rm Re}_{p}<10$ agreed with those in other studies (Bagchi and Balachandar, 2002; Bluemink et al., 2010). Therefore, in the present cases, the proposed drag estimation method reflects the disturbance flow that contributes to the history and the lift forces.

To show the applicability of VA simulation for the study of turbulent modulation, the simulation for the Taylor-Green vortex at Reynolds numbers ${\rm Re}=15$ and $18$ was carried out with the particle diameter being $\mathcal{O}(10)$ times smaller than the vortex. For density ratio $\rho_{d}/\rho_{c}\geq 10$, the particle motion obtained by the VA simulation showed much better agreement with that of the fully-resolved simulation than the traditional two-way coupling simulation. The disturbance flow also showed good agreement with that of the fully-resolved simulation. On the other hand, for density ratio $\rho_{d}/\rho_{c}=1$, the VA simulation model needs to be improved. For a further improvement of the estimation of the fluid force, unsteadiness and non-uniformity of the flow need to be considered. However, we consider the method as promising as the computational cost of the VA simulation is $\mathcal{O}(10^{4})$ times lower than that of the fully-resolved IBM simulation in the present paper.

The importance of two-way coupling for the correct prediction of particle trajectories in vortical flows was confirmed for $\rho_{d}/\rho_{c}=100$. For particles released in a vortical array, the trajectory curvature in the initial stage increased due to the history effects, which clearly influenced the future dynamics. The history effect estimated in the VA simulation tends to be larger for larger particle as supported by Bergougnoux et al. (2014) and Daitche (2015). For particles initially placed at the vortex centre, the long residence time around the initial position increases the importance of the role of the history effects on the trajectory. It is also found that the particle interactions, assuming an average inter-particle distance of about $10D$, influence the particle motion in vortical flows. These results suggest that the history effects and inter-particle flow-mediated interactions need to be considered by two-way coupling simulations even in dilute particle-laden turbulence.

Acknowledgements

One of the authors, T.F., gratefully acknowledges the financial support of the Japan Society for the Promotion of Science (JSPS) KAKENHI Grant No. 15J00439. L.B. and W.F. acknowledge financial support from the European Research Council grant no. ERC-2013- CoG-616186, TRITOS, by the Swedish Research Council (VR) and computer time provided by SNIC (Swedish National Infrastructure for Computing). This work is partly supported by Grant-in-Aid (B) No.16H04271 and No.17H03174 of the JSPS. The authors gratefully acknowledge the financial support by Mr. Yibin Lin (Hangzhou, China).

Appendix A: Differentiability of volume-averaged quantities

The form of the viscous term in Eq. (6) is justified in this section. In the following discussion, $V$ is assumed to be larger than the particle as used in this paper. The volume integral $Q(\bm{x})$ is defined as

[TABLE]

where $q$ is a bounded function defined in both fluid and solid. We also assume that Taylor expansion of $q$ is possible except on the interface. The volume-averaged quantities correspond to $Q(\bm{x})/V$. For example, $\alpha_{c}\left<u_{1}\right>_{c}$ is constructed from

[TABLE]

and $(\alpha_{c}\left<u_{1}\right>_{c}+\alpha_{d}\left<v_{1}\right>_{d})$ corresponds to

[TABLE]

In the following, the first and second-order derivatives of $Q(\bm{x})$ are considered.

Figure A1 schematically shows the geometric relation between $V(\bm{x})$ and $V(\bm{x}+h\bm{e}_{i})$, and we focus on the integration of $q$ over the volume denoted as $Q(\bm{x})$ and $Q(\bm{x}+h\bm{e}_{i})$. The volume integrals of $q$ in the shaded regions are denoted as $Q^{+}$, $Q^{-}$ and $Q_{\rm cut}$. The surface of $V(\bm{x})$ is denoted as $S$. The outward unit normal vector on $S$ is denoted as $\bm{n}_{V}$. The surface $S$ is divided into a region denoted as $S^{+}$, where $\bm{n}_{V}\cdot\bm{e}_{i}\geq 0$, and $S^{-}$, defined by $\bm{n}_{V}\cdot\bm{e}_{i}\leq 0$. According to Fig. A1, the volume integrals are

[TABLE]

where $\bm{S}$ indicates the position on $S$. Therefore, we obtain the following equation:

[TABLE]

To deal with the interface between the phases, we define the fraction of the surface $S_{\rm jump}(\subset S)$ such that

[TABLE]

(see Fig. A1). With this decomposition of the surface $S$,

[TABLE]

Taking the limit of $h\to 0$, $S_{\rm jump}$ converges to 0 and

[TABLE]

on $S_{\rm jump}$ is bounded since $q$ is bounded. We can therefore write

[TABLE]

As Eq. (A.9) holds regardless of the sign of $h$, we obtain the derivative

[TABLE]

Note that $S_{\rm jump}\to 0$ is guaranteed by the size difference between $V$ and the particle. In this case, the continuity of $q$ on the interface is not necessary for the first-order derivative.

For simplicity, ${\partial Q(\bm{x})}/{\partial x_{i}}$ is denoted as $Q_{,i}(\bm{x})$. For higher-order derivatives, we consider that

[TABLE]

so that

[TABLE]

By denoting the interface as $(\bm{S}+k(\bm{S})\bm{e}_{i})$ with $0\leq k\leq h$, we obtain

[TABLE]

where $[\cdot]$ represents the jump of the function at the interface. Taking the limit of $h\to 0$, Eq. (A.12) yields:

[TABLE]

In general, the right-hand side of Eq. (A.14) depends on the sign of $h$ (e.g., the second-order derivative of $\alpha_{d}$ is not determined). On the other hand, when $[q(\bm{S}+k\bm{e}_{i})]=0$, we can define the second-order derivative as

[TABLE]

As the velocities $\bm{u}$ and $\bm{v}$ are continuous across the interface, the viscous term $\nabla^{2}(\alpha_{c}\left<\bm{u}\right>_{c}+\alpha_{d}\left<\bm{v}\right>_{d})$ is well-defined. Note, however, that the decomposition into $\nabla^{2}(\alpha_{c}\left<\bm{u}\right>_{c})+\nabla^{2}(\alpha_{d}\left<\bm{v}\right>_{d})$ is not allowed.

Appendix B: Calculation of interaction force

The interaction force can be written as

[TABLE]

where $\bm{s}$ is the stress on the surface:

[TABLE]

In Sec. 2.2, the stress vectors:

[TABLE]

and

[TABLE]

are used for the modeling of Eqs. (18) and (19). Therefore, the following two integrals are enough for the derivation of the interaction force models:

[TABLE]

Note that these integrals are also enough for the interaction force modelling in Fukada et al. (2016).

We consider three unit vectors $\bm{A}$, $\bm{B}$ and $\bm{C}$ with $\bm{A}={(\bm{x}-\bm{x}_{p})}/{|\bm{x}-\bm{x}_{p}|}$ and $\bm{A}\cdot\bm{B}=\bm{B}\cdot\bm{C}=\bm{C}\cdot\bm{A}=0$. The basis vector $\bm{e}_{i}$ and $\bm{n}$ can be written as

[TABLE]

so that the integrands become

[TABLE]

Given the symmetry about the direction $\bm{A}$, the integrals reduce to

[TABLE]

where $t$ is the angle between $\bm{A}$ and $\bm{n}$. In this derivation, we used the following relations:

[TABLE]

According to Eqs. (D.10) and (D.11) (shown later), we obtain

[TABLE]

Appendix C: Estimation of fluid force

The applicability of the current approximation of the fluid force is tested for unsteady flows. For the flow fields around the particle obtained by the fully-resolved simulation, the contributions of $F_{\rm drag}$ and $\nabla P_{\rm ud}$ on the particle acceleration (28) are computed by Eqs. (13), (20) and (26). Instead, the volume-averaged values are directly computed from the flow fields. In particular, we consider time $(A/L)t=33.13$ for the cases of $\rho_{d}/\rho_{c}=1,10,1000$ in Sec. 4.3 and time $(A/L)t=55.22$ for the case $\rho_{d}/\rho_{c}=100$ in Sec. 4.4.1. Table 1 shows the following three dimensionless accelerations (for comparing the contribution of each term to the right-hand side of Eq. (28)):

[TABLE]

where $d\bm{v}_{p}/dt$ is the net acceleration obtained from the fully-resolved IBM simulation. Note that $\bm{a}^{sv+pg}$ is equal to $\bm{a}^{FR}$ when the errors in both the model (28) and the estimation of each term are ignored. Especially for $\rho_{d}/\rho_{c}\leq 100$, the effect of the pressure gradient is reasonably captured by Eq. (26) as $\bm{a}^{sv+pg}$ is considerably improved from $\bm{a}^{sv}$. However, as indicated by the differences between $\bm{a}^{sv+pg}$ and $\bm{a}^{FR}$, unsteadiness and non-uniformity of the flow need to be considered to improve the total estimation method.

Appendix D: Calculations of geometrical functions

The radius $R$ of the averaging volume $V$ is larger than the particle radius $r_{d}$. The geometrical functions $\xi$, $\alpha_{d}$ and $\alpha_{d}\left<\bm{v}\right>_{d}$ are obvious for $|\bm{x}-\bm{x}_{p}|<R-r_{d}$ and $|\bm{x}-\bm{x}_{p}|>R+r_{d}$. Therefore, we only consider $R-r_{d}\leq|\bm{x}-\bm{x}_{p}|\leq R+r_{d}$. Figure D1 shows the definitions of the variables considered in the following. The origin is at the particle centre and $y$ is equal to $|\bm{x}-\bm{x}_{p}|$. The variable $a$ satisfies

[TABLE]

Note that $a$ varies in the range $-r_{d}\leq a\leq r_{d}$. The surface fraction $\xi=S_{d}/\pi D^{2}$ is

[TABLE]

The volume fraction $\alpha_{d}=V_{d}/V$ is

[TABLE]

The centre of gravity $x_{G}$ of $V_{d}$ is

[TABLE]

In general, the rigid-body velocity at $\bm{x}$ can be written as

[TABLE]

where $\bm{v}_{0}$ and $\bm{\Omega}_{0}$ are the origin velocity and the angular velocity around the origin. Introducing

[TABLE]

where $V_{r}$ is the volume of the rigid body, the average velocity in the volume becomes

[TABLE]

Therefore, the averaged velocity $\alpha_{d}\left<\bm{v}\right>_{d}$ is

[TABLE]

Eqs. (32), (33) and (34) are obtained by Eqs. (D.3), (D.4), (D.5) and (D.9) using Eq. (D.2).

Finally, the integrals used in Appendix B are calculated as follows:

[TABLE]

Appendix E: Effects of size of the averaging volume and grid width

To investigate the effect of the size of the averaging volume, the VA simulation is expanded for the case of $R=1.5D$. For this case, the fitting functions (20), (23) and (24) are replaced by

[TABLE]

[TABLE]

Based on these functions, the VA simulations under $R=1.5D$ and different $\Delta x$ are carried out for the same configurations as in Secs. 4.1 and 4.4.1.

Figure E1 shows the result for the particle settling problem in a stationary fluid. As for the case of $R=0.75D$, the result of the fine grid ($D/\Delta x=4$) is almost the same as that for $D/\Delta x=2$. On the other hand, for the coarse grid ($D/\Delta x=1$), the result is quite different from the others. Therefore, the grid resolution of $D/\Delta x=2$ is necessary to capture the averaged flow distribution around the particle. As for the case of $R=1.5D$, the result of the coarse grid ($D/\Delta x=1$) is not so different from that of $D/\Delta x=2$. Therefore, choosing a larger $R$ is better for $D/\Delta x=1$, while having a smaller $R$ is appropriate for $D/\Delta x=2$. However, to outperform the O-NL simulation (shown in Fig. 6 (a)), a fine grid with small $R$ ($=0.75D$) is necessary.

Figure E2 shows the result for the particle under gravity in an array of Taylor-Green vortices. The tendency of the effects of $R$ and $\Delta x$ are quite similar to that for the settling particle in the stationary fluid. The results for $R=1.5D$ ($D/\Delta x=1$ and 2) are still better than that of the O-NL simulation. Therefore, for the larger $R$, the history effect is qualitatively captured in the same manner as for the smaller $R$. The upper limit of $R$ is considered to be based on the length scale of the background flow. In conclusion, smaller $R$ gives better results for sufficiently fine grids. On the other hand, for coarse grid, larger $R$ is preferable.