A Dynamic Anisotropic Spatially-Averaged Two-Fluid Model for Moderately Dense Gas-Particle Flows

A Dynamic Anisotropic Spatially-Averaged Two-Fluid Model for Moderately Dense Gas-Particle Flows

Journal Pre-proof

A Dynamic Anisotropic Spatially-Averaged Two-Fluid Model for Moderately Dense Gas-Particle Flows Stefanie Rauchenzauner, Simon Schneiderbauer PII: DOI: Reference:

S0301-9322(19)30488-4 https://doi.org/10.1016/j.ijmultiphaseflow.2020.103237 IJMF 103237

To appear in:

International Journal of Multiphase Flow

Received date: Revised date: Accepted date:

4 July 2019 22 January 2020 26 January 2020

Please cite this article as: Stefanie Rauchenzauner, Simon Schneiderbauer , A Dynamic Anisotropic Spatially-Averaged Two-Fluid Model for Moderately Dense Gas-Particle Flows, International Journal of Multiphase Flow (2020), doi: https://doi.org/10.1016/j.ijmultiphaseflow.2020.103237

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2020 Published by Elsevier Ltd.

Highlights • An anisotropic turbulence model for moderately dense gas-particle flows. • Anisotropic closure models for the Reynolds-stresses are applied. • Closure models for the individual components of the drift velocity are derived. • Dynamic adjustment of the correlation coefficients in coarse-grid simulations. • An a priori validation shows good agreement with highly resolved simulation data.

1

A Dynamic Anisotropic Spatially-Averaged Two-Fluid Model for Moderately Dense Gas-Particle Flows Stefanie Rauchenzauner∗, Simon Schneiderbauer Christian-Doppler Laboratory for Multi-scale Modelling of Multiphase Processes, Department of Particulate Flow Modelling, Johannes Kepler University, Altenbergerstrasse 69, 4040 Linz, Austria

Abstract We present an anisotropic, dynamic multi-phase turbulence model for moderately dense gas-particle flows by spatially averaging the kinetic-theory based two-fluid model (Schneiderbauer, AIChE J., 2017; 63(8): 3544 – 3562). The filtered gasparticle drag force can be approximated by the resolved drag force corrected by the drift velocity (Parmentier et al., AIChE J., 2012; 58(4): 1084 – 1098). The drift velocity can be expressed as a correlation between the gas-phase velocity and the solid volume fraction. We propose to calculate the correlation coefficients locally and dynamically by application of a scale-similarity approach. Therefore, we show that test-filters can be employed to estimate the correlation coefficients in coarse-grid simulations. Furthermore, transport equations for the components of the highly anisotropic Reynolds-stress tensor are derived, and closure models known from single-phase LES modelling are applied to the unresolved terms with the exception of the interfacial work term, which is expressed by the correlations between the velocities of the phases. In addition, the cluster-induced turbulence production term (Capecelatro et al., J. Fluid Mech., 2015; 780: 578 – 635) arising in the gas-phase Reynolds-stress transport equation is closed using the drift velocity. An a priori validation of the developed closure models against filtered fine-grid simulation data of unbound fluidization for Geldart type A and B particles, as well as an a posteriori verification in wall-bounded fluidized beds of Geldart type A and B particles ∗ Corresponding author at: Christian-Doppler Laboratory for Multi-scale Modelling of Multiphase Processes, Department of Particulate Flow Modelling, Johannes Kepler University, Altenbergerstrasse 69, 4040 Linz, Austria Email address: [email protected] (Stefanie Rauchenzauner)

Preprint submitted to International Journal of Multiphase Flow

January 27, 2020

are conducted. Keywords: cluster-induced turbulence, multi-phase turbulence, turbulent kinetic energy, kinetic theory based two-fluid model, coarse-grid simulations 1. Introduction Gas-solid flows, like fluidized beds or risers, are used for a wide range of applications, from food drying and polymer production over catalyst cracking to the direct reduction of iron-ore. In addition to experiments, simulations have proven to be a useful tool in the design of such processes. While for small-scale fluidized beds, simulation is quite feasible, the occurring physical phenomena cannot be fully analyzed in full-scale reactors yet, due to the large range of involved scales. In the presence of a mean body force, such as gravity, mesoscale structures form due to the momentum coupling between the two phases. More precisely, variations in the solid volume fraction and fluctuations of the gas-phase velocity lead to clustering. These mesoscale structures, in turn, are a source for turbulence in the gas-phase, giving rise to a class of multiphase turbulence referred to as cluster induced turbulence (CIT) (Capecelatro et al., 2014, 2015, 2016a,b, 2018; Fox, 2014). Typically, in fluidized bed reactors with expansions of several meters, these smallest stable structures are only several particle diameters wide (Agrawal et al., 2001). The microscopic behaviour of individual particles, i.e. the gas-solid drag force, the intra- and interphase heat transfer between the particles and the gas, as well as collisions cannot be fully resolved due to computational limitations. The kinetic theory based two-fluid model (TFM) (Anderson and Jackson, 1967; Ishii, 1975; Lun et al., 1984), which consists of averaged equations of motion on the micro-scale has proven to predict the macro-scale hydrodynamics of gas-solid flows in a reasonable way. However, coarse numerical grids are used to further reduce the demand on computational resources. These coarse grids, in turn, do not resolve all the mesoscale structures, such as clusters and streamers. These unresolved heterogeneous mesoscale structures, however, have a significant influence on the macro-scale flow properties. An overestimation of the gas-solid drag force by not accounting for the unresolved mesoscale structures leads, for example, to an overestimation of the bed expansion (Schneiderbauer et al., 2013). In order to overcome this problem,

3

many sub-grid drag modifications have been proposed in literature. A filtered TFM approach was used in multiple previous studies (Capecelatro et al., 2014, 2015, 2016a,b; Cloete et al., 2018, 2019; Fox, 2014; Igci et al., 2008; Igci and Sundaresan, 2011; Milioli et al., 2013; Ozel et al., 2013, 2017; Parmentier et al., 2012; Schneiderbauer, 2017, 2018a,b; Schneiderbauer and Pirker, 2014; Schneiderbauer et al., 2013). Thereby, constitutive relations for the unresolved terms appearing in the filtered TFM balance equations are derived in different ways. In a large number of studies, these sub-filter modifications are deduced from either highly resolved TFM or Euler-Lagrange simulations. The simulation data is averaged using filters of different sizes and binned in terms of markers characterizing the sub-filter-scale state, such as the solid volume fraction and the slip velocity (Igci et al., 2008; Igci and Sundaresan, 2011; Milioli et al., 2013; Schneiderbauer and Pirker, 2014). Recently, anisotropic closure models were introduced by Cloete et al. (2018, 2019) in this framework. A different approach is to use concepts of turbulence modelling to derive constitutive relations for the unresolved terms. In this context, Fox (2014) and Capecelatro et al. (2014, 2015, 2016a,b, 2018) applied the concept of Reynolds-averaged NavierStokes (RANS) modelling, Schneiderbauer (2017, 2018a,b) used spatial filtering to derive a large eddy simulation (LES) type model, the spatially-averaged two-fluid model (SA-TFM). In addition to functional turbulence models, structural turbulence models have recently been employed for moderately-dense gas-particle flows. Schneiderbauer and Saeedipour (2018, 2019), for example, showed that the approximate deconvolution model (ADM) (Stolz and Adams, 1999) is capable of predicting the influence of CIT on the filtered gas-particle drag force. Furthermore, it was previously shown that the filtered drag force can be approximated by the resolved drag force reduced by a so-called drift velocity, the gas-phase velocity fluctuations seen by the particles (Capecelatro et al., 2016a,b, 2018; Fox, 2014; Ozel et al., 2013, 2017; Parmentier et al., 2012; Schneiderbauer and Saeedipour, 2018). Remarkably, this assumption is not only applicable in the case of low particle Reynolds number and low particle mass loading (Capecelatro et al., 2014, 2015, 2016a,b, 2018; Fox, 2014), but also in moderately dense gas-particle flows (Ozel et al., 2013, 2017; Parmentier et al., 2012; Schneiderbauer and Saeedipour,

4

2018). In other words, the drift velocity is a measure for the drag reduction due to the unresolved mesoscale structures. As discussed above, due to the large range of involved scales, up-scaling of simulations, which are verified against fine-grid simulations or table-top experiments, is non-trivial. Therefore, the models also need to be verified against experimental results obtained from full-size gas-particle flow reactors. One of the most prominent examples is the NETL challenge problem (Panday et al., 2014), which consist of pressure-gradient and velocity measurements inside the riser section of a circulating fluidized bed of Geldart type A and B particles under various operating conditions. Additional experimental data on the riser section of a circulating fluidized bed is, for example, provided by Issangya et al. (1999, 2000); Pärssinen and Zhu (2001a,b); Weber et al. (2018). Experimental data on dense bubbling fluidized beds can, for example be found in Zhu et al. (2008). In this work, we propose a dynamic closure model for the drift velocity as a function of the anisotropic components of the turbulent kinetic energy (TKE) of the gas-phase and the variance of the solid volume fraction. The model includes dynamic estimation of correlation coefficients and model parameters by application of testfilters in coarse-grid simulations. We discuss the derivation of a transport equation for the Reynolds-stresses, i.e. the anisotropic components of the TKE, as well as the scale-similarity of the correlation coefficients. The developed models are validated against highly resolved TFM simulation data of unbound fluidization for Geldart type A and B particles in an a priori analysis. Finally, an a posteriori verification of the developed models in two cases of 3-dimensional wall-bounded fluidized beds of Geldart type A and B particles is provided. 2. Case Description Following previous work (Milioli et al., 2013; Parmentier et al., 2012; Schneiderbauer, 2017, 2018a,b; Schneiderbauer et al., 2012), we performed highly resolved simulations of Geldart group type A and type B particles in three-dimensional periodic boxes. The main physical properties relevant for this study are summarized in Table 1. We employed a kinetic theory based TFM (Lun et al., 1984; Schneiderbauer et al., 2012). The gas-phase pressure gradient was adjusted to balance the gravi5

tational acceleration of the mixture density. In dense regions, i.e., where the solid volume fraction is close to the maximum packing conditions and the inter-particle forces are dominated by long enduring multiple frictional contacts, the solid stress was closed by using an inertial number dependent rheology (Schneiderbauer et al., 2012). The convective terms in the momentum equations were discretized by using a second-order upwind scheme. The derivatives appearing in the viscous terms were computed by a least squares method, and the pressure-velocity coupling was achieved by the SIMPLE algorithm, whereas the face pressures were computed as the average of the pressure values in the adjacent cells (linear interpolation). Snapshots of the flow field were collected at various times after the system reached a statistically steady state with persistent temporally and spatially distributed clusters (Milioli et al., 2013; Schneiderbauer, 2017, 2018a). Following this procedure, we obtained about 15 million data points for the evaluation of the presented models. The data was filtered using box filters of different sizes. We made the filter-size dimensionless by using the characteristic length scale of mesoscale structures Lch , see Table 1. We used a grid spacing of 3 · ds and 7 · ds for the Geldart type A and B particles,

respectively. This grid spacing is assumed to be fine enough to resolve all rele-

vant heterogeneous structures (Fullmer and Hrenya, 2016; Schneiderbauer, 2018a; Schneiderbauer and Saeedipour, 2019; Uddin and Coronella, 2017). 3. Spatially-Averaged TFM Approach The TFM equations as formulated by Ishii for dispersed two-phase flow and described by Enwald for dispersed gas-solid flow (Enwald et al., 1996; Ishii, 1975; van Wachem et al., 2001) are suited for application in moderately dense gas-particle flows. Therefore, they are used as a base for the derivation of the SA-TFM (Schneiderbauer, 2017). The continuity equations read, ∂φ ∂ + (φuk ) = 0, ∂t ∂xk ∂(1 − φ) ∂ + ((1 − φ)vk ) = 0, ∂t ∂xk

6

Property

Case 1

Case 2

φmax

0.6

φfr

0.4

µst i

0.28

ρg

1.224 kg m−3

µg

1.78 · 10−5 kg m−1 s−1

ds

150 · 10−6 m

75 · 10−6 m

ρs

2500 kg m−3

1500 kg m−3

ut

0.96 m s−1

0.218 m s−1

Fr

626

65

St

5.6 · 10−5

3.1 · 10−5

Lch hφid ˆd ∆

1.3 · 10−3 m

0.1

3 · 10−4 m

256 x 64 x 1024

Table 1: Material Properties and Operating Conditions for the Fine Grid Simulations (Schneiderbauer, 2018a; Schneiderbauer and Pirker, 2014; Schneiderbauer et al., 2013). Note that these properties correspond to Geldart group type B and A particles (φmax : solid volume fraction at maximum packing conditions; φfr : threshold solid volume fraction for which the frictional stress becomes active; hφid : domain averaged solid volume fraction; µst i : statistic coefficient of internal

friction; ds : particle diameter; ρs : particle density; ρg : gas density; µg : gas dynamic viscosity; ut : terminal settling velocity; Fr: particle Froude number, Fr = u2t /(ds g), where g is the magnitude of the gravitational acceleration; St: particle Stokes number, St = τp /τf , where τp is the particle relaxation time, τp = 1/hβi, where the average is taken with respect to the dispersed phase (Février et al., 2005), and τf = kg /g is the fluid Lagrangian time scale (Haworth and Pope, 1986; Moreau et al., 2010); Lch : characteristic length scale of mesoscale structures (Radl and Sundaresan, 2014), ˆ d : domain size made dimensionless using Lch ). Lch = (u2t /g)Fr−2/3 ; ∆

for the solid- and gas-phase, respectively, whereas the balance of momentum equations are given by:
∂ρs φui ∂ ∂p ∂ fr + (ρs φui uk ) = −φ + Σkc ik + Σik + β(vi − ui ) + φρs gi , ∂t ∂xk ∂xi ∂xk

7

(1)

∂ ∂p ∂ g ∂ρg (1 − φ)vi + (ρg (1 − φ)vi vk ) = −(1 − φ) + Σ ∂t ∂xk ∂xi ∂xk ik

(2)

− β(vi − ui ) + (1 − φ)ρg gi .

Here, φ denotes the solid volume faction, ρs and ρg the particle- and gas-phase densities, ∂p/∂xi the gas-phase pressure gradient, which is distributed over the phases, g fr Σkc ij the particle-phase kinetic theory stresses, Σij the frictional stresses, Σij the

components of the gas-phase stress tensor (Table 2) and β the microscopic drag coefficient (Schneiderbauer et al., 2012). We employed the Beetstra et al. (2007) drag model for the microscopic gas-particle drag, which is suitable for moderately dense flow. The spatially averaged, or filtered, part of a function ζ(x, t) is given by Z ζ(x, t) = ω(x, y, ∆f ) ζ(y, t) dVy , V

with a filter kernel ω(x, y, ∆f ), which satisfies Z ω(x, y, ∆f )dVy = 1, V

for every point in space x and every control volume V , whose elongations are larger than the filter size ∆f . Using this definition, we denote the Favre averages of the solid- and gas-phase velocity components, ui and vi , respectively, as φui , φ (1 − φ)vi hvi ig = . (1 − φ)

hui is =

The sub-filter fluctuations of the variables can, thus, be defined as: φ0 = φ − φ, u00i = ui − hui is , vi000 = vi − hvi ig . Note that, unlike for Reynolds-averaging, the spatial averages do not fulfill the identity ζ = ζ in general. Therefore, we will define the central moments of the velocities similar to Germano (1992): hvi000 vk000 ig := hvi vk ig − hvi ig hvk ig , hu00i u00k is := hui uk is − hui is huk is , φ0 vi000 := φvi − φvi . 8

(4)

Transport equation for the pseudo-thermal energy: 3 2



 ∂φρs Θ ∂ ∂qk ∂ui + (φρs uk Θ) = Σkc − + Γ s − J v − γΘ . ik ∂t ∂xk ∂xk ∂xk

(3)

Radial distribution function: g0 = min



1 1 φ2 3 φ 1 , + + 1 − φ/φmax 1 − φ 2 (1 − φ)2 2 (1 − φ)3



.

Gas-phase and solids phase stress tensors: Σgik

= 2(1 − g Sik

g φ)µg Sik ,

Σkc ik

1 = 2





µ∗ g0 ηs (2 − ηs )

= 

∂vi ∂vk + ∂xk ∂xi

s 2µkc s Sik

−

+

1 ∂vk , 3 ∂xk



∂uj λkc s ∂xj

s Sik =

1 2

− 

pkc s



fr s fr Σfr ik = 2µs Sik − ps I,  ∂ui ∂uk 1 ∂uk + − . ∂xk ∂xi 3 ∂xk

I,

Solids viscosity: 

µkc s = µ∗ =

2+α 3

   1 8 8 3 + φηs g0 1 + ηs (3ηs − 2)φg0 + ηs µb , 1 + ls /L 5 5 5 √ 2 5ρs ds πΘ 256µφ g0 8 1 µ= , µb = , α = , ηs = (1 + es ). 96 5π 5 2

µ , 2βµ 1+ (φρs )2 g0 Θ



Pseudo-thermal energy (PTE) flux vector qj , rate of dissipation of PTE γΘ , rate of dissipation of PTE by viscous damping Jv and rate of production of PTE by gas-particle slip Γs : 
12 2 64 ∂Θ ηs (4ηs − 3)φg0 + (41 − 33ηs )ηs2 φ2 g02 , 5 25π ∂xj √ ρs φ2 48 κ 75ρs ds πΘ 3/2 ∗ g0 Θ , κ = , κ= , γΘ = √ ηs (1 − ηs ) 6βκ ds 48ηs (41 − 33ηs ) π 1+ 5(φρs )2 g0 Θ p 81φµ2g (vi − ui )2 √ , Jv = 3βΘ. Γs = g0 d3s ρs πΘ

qj = −

κ∗ g0



12 1 + ηs φg0 1 + ls /L 5



1+

Solids pressure and bulk viscosity: pkc s

= φρs



 1 + 4ηs φg0 Θ, 1 + ls /L

λkc s

8 = ηs φ2 ρs ds g0 3

r

Θ . π

Frictional pressure and viscosity: 

 pfr s = 4ρs

µi (Is )pfr s

fr s µfr , s (Is , ps , S ) = q s S s /2 2 Sij ij

bds

q

s S s /2 Sij ij

φmax − φ

µi (Is ) = µst i +

µci

2

 ,

µst i

− , I0 /Is + 1

q s S s /2d Sij s ij p Is = . pfr s /ρs 2

Table 2: Kinetic theory of granular flows constitutive equations (Schneiderbauer et al. (2012) and References cited therein).

9

Consequently, we can define the sub-filter scale turbulent kinetic energies of both phases as, 1 ks = hu00i u00i is , 2 1 kg = hvi000 vi000 ig . 2

(5)

Together with the filtered kinetic energies these constitute the total kinetic energies of the phases. By applying spatial filters to the TFM balance equations and using the definitions of the central moments (4), we obtain balance equations for the filtered quantities: ∂φ ∂ + (φhuk is ) = 0, ∂t ∂xk ∂(1 − φ) ∂ + ((1 − φ)hvk ig ) = 0, ∂t ∂xk ∂ρs φhui is ∂ + (6) (ρs φhui is huk is ) ∂t ∂xk  ∂p ∂ ∂  kc fr = −φ Σik + Σik + β(vi − ui ) + φρs gi − + (ρs φhu00i u00k is ), ∂xi ∂xk ∂xk ∂ρg (1 − φ)hvi ig ∂ + (ρg (1 − φ)hvi ig hvk ig ) ∂t ∂xk ∂p ∂ g = − (1 − φ) Σ + ∂xi ∂xk ik − β(vi − ui ) + (1 − φ)ρg gi −

(7) ∂ (ρg (1 − φ)hvi000 vk000 ig ). ∂xk

The filtered balance equations have a form similar to the microscopic TFM balance equations, with the filtered variables in place of the microscopic ones. The last terms in the momentum balance equations are Reynolds-stress-like contributions. In the following, we discuss closure models for the individual unresolved terms. 3.1. Dynamic Adjustment It has to be emphasized that we propose a dynamic estimation of the arising model parameters by application of the scale-invariance theory (Pope, 2000; Schneiderbauer, 2017) in this study. Similarly, we argue that the coefficients and constants present in the models do not depend on the specific choice of the weight function in the filter operation, as they do not depend on the filter-size (Ozel et al., 2017). 10

Therefore, these coefficients can be determined dynamically by applying test-filters (Germano et al., 1991; Lilly, 1992) to the resolved variables in coarse-grid simulations. In this work, in accordance with Parmentier et al. (2012), we used a test-filter, which is defined in 3-dimensions as: 1 b ζ(x, t) = ζ(x, t) + ζ(x + ∆f ex , t) + ζ(x − ∆f ex , t) 7

+ ζ(x + ∆f ey , t) + ζ(x − ∆f ey , t)  + ζ(x + ∆f ez , t) + ζ(x − ∆f ez , t) ,

(8)

where ei is the unit vector in i-th direction and ∆f is the filter-size. The filter-size of this test filter is approximately 2 times the size of the grid filter in the case of a Gaussian filter (Saeedipour et al., 2019). Therefore, it fulfills the requirements for test filters according to Germano et al. (1991). 3.2. Closure Models 3.2.1. Filtered Stresses The Reynolds-stresses can, for example, be approximated using an isotropic Boussinesq-like hypothesis (Schneiderbauer, 2017, 2018a). They are, however, highly anisotropic (Cloete et al., 2018; Schneiderbauer and Saeedipour, 2018). Therefore, we will derive transport equations for each of the diagonal components of the Reynolds-stresses and close the shear components as correlations between the diagonal components. This will be discussed in more detail in Section 5. The filtered microscopic phase stress tensors were shown to be negligible compared to the Reynolds-stress like velocity fluctuations for sufficiently large filter sizes in various studies (Igci et al., 2008; Igci and Sundaresan, 2011; Schneiderbauer, b f  1 in this study, such that we do 2017), therefore, we restrict ourselves to ∆ not need to solve the transport equation for the granular temperature. The dom-

inance of the turbulent viscosity over the filtered kinetic-theory viscosity increases with increasing filter-size (Igci and Sundaresan, 2011). A dimensionless filter-size b f = 3 does not necessarily satisfy the constraint, however it was shown by Igci of ∆

and Sundaresan (2011) and Schneiderbauer and Pirker (2014) that even in this case

the filtered kinetic-theory stresses are small compared to the Reynolds-stresses and that solving the granular temperature equation (3) would unnecessarily slow down the coarse-grid simulations. 11

3.2.2. mesoscale Interphase Force The term involving the gas-phase pressure gradient is referred to as the mesoscale interphase force (Cloete et al., 2018; Zhang and VanderHeyden, 2002). Similar to previous publications (Ozel et al., 2017; Sarkar et al., 2016; Schneiderbauer, 2017), our data suggests that it is small compared to the filtered micro-scale interphase force, i.e. the filtered drag force, in all considered Cases. It is found to be about 10% of the total interphase force acting on the mesoscale. A closure model was proposed by De Wilde (De Wilde, 2005) by decomposing the gas-phase pressure gradient into its static and dynamic components, ∂ ∂ p = ρm gi + p˜, ∂xi ∂xi

(9)

ρm = φρs + (1 − φ)ρg .

(10)

with the mixture density,

Using this identity and the averaged solid mass fraction, ms =

φρs , ρm

he showed that the vertical component of the mesoscale interphase force has a direct contribution of the form, φ0

∂ 0 φ02 ∂ p ≈ ms p, ∂z ∂z φ

(11)

and that the term involving the gradient of p˜0 has a negligible contribution (De Wilde, 2005; Fox, 2014). In Figure 1, it can be seen that this closure model (11) fits reasonably well with the predictions obtained by filtering the fine grid simulation data in the considered Cases. The model, however, seems to perform worse in Case 1. This could be because p˜0 is not entirely negligible for Geldart type B particles which have a larger Froude number or larger particle diameter. In wall-bounded flows, the non-hydrostatic part of the pressure gradient p˜0 is not negligible, since the particles cannot fall freely in regions close to the walls. Therefore, one might need to introduce an additional model parameter, which estimates the importance of p˜0 (similar to the pressure dilatation, which is discussed in Section 4.1.5). In coarse-grid simulations, this model, together with the constitutive relation for the variance of the solid volume fraction, equation (48) (Schneiderbauer, 2017), can be 12

used to approximate the mesoscale interphase force. In addition, this term could also be implemented as a virtual added mass force (De Wilde, 2005; Zhang and VanderHeyden, 2002), which will be investigated further in a future a posteriori study. a)

b) 500

1400 1200

400

1000 800

300

600

200

400 100

200 0

0

0.2

0.4

0.6

0.8

0

1

0

0.2

0.4

0.6

0.8

1

Figure 1: Closure model for the magnitude of the mesoscale interphase force in vertical direction ∂ 0 φ0 ∂z p compared to the predictions obtained by filtering the fine grid simulation data for different

ˆ f , as a function of the normalized filtered solid volume fraction, a) for dimensionless filter sizes ∆ Case 1 and b) for Case 2. The symbols correspond to measurements from the filtered fine grid simulation data, the lines correspond to the predictions by employing equation (11).

3.2.3. Filtered Gas-Particle Drag Force Various different previous studies (Capecelatro et al., 2014, 2015, 2016a,b, 2018; Fox, 2014) have shown that the microscopic drag coefficient divided by the solid volume fraction can be approximated by its zeroth order Taylor series expansion about the filtered variables ((1 − φ), hvig , huis ), in the case of low particle Reynolds

number and low particle mass loadings. However, it was also shown that this holds true in moderately dense gas-particle flows, where i.e. the Wen and Yu (1966) drag law or the Beetstra et al. (2007) drag law is applicable (Ozel et al., 2013, 2017; Parmentier et al., 2012; Schneiderbauer and Saeedipour, 2018). Thus, the filtered drag force can be approximated as β(vi − ui ) ≈

βe φ(vi − ui ) φ

e i ig − hui is + vd,i ), = β(hv 13

(12)

with the microscopic drag coefficient evaluated at the filtered variables βe = β((1 − φ), hvig , huis ) and the drift velocity vd,i , defined as (Parmentier et al., 2012):

vd,i := hvi is − hvi ig . The drift velocity represents the sub-filter gas-phase velocity fluctuations seen by e i ig − hui is ) severely overpredicts the filthe particles. The resolved drag force β(hv tered drag force on average, since it does not account for the sub-filter, mesoscale

heterogeneous structures (Schneiderbauer et al., 2013). In order to provide comparability with relevant publications, i.e. (Igci and Sundaresan, 2011; Milioli et al., 2013; Sarkar et al., 2016), we use the fractional correction function Hi , defined as, Hi = 1 −

β(vi − ui )

e i ig − hui is ) β(hv

,

which is a measure for the sub-grid heterogeneity implying a reduction of the resolved drag force: e i ig − hui is ). β(vi − ui ) ≈ (1 − Hi )β(hv

Using the approximation for the filtered drag force (12), we can express the fractional correction function as: Hi = 1 −

e i ig − hui is + vd,i ) β(hv . e i ig − hui is ) β(hv

(13)

The drift velocity can be expressed as a function of the covariance of the solid volume fraction and the gas-phase velocity (Cloete et al., 2019): vd,i =

φ0 vi000 . φ(1 − φ)

(14)

Following previous studies (Schneiderbauer, 2017; Schneiderbauer and Saeedipour, 2018), we express the covariance of two variables by the square root of the variances ∗ of the variables scaled by a linear correlation coefficient ξφg , such that the magnitude

of the drift velocity can be written in terms of the turbulent kinetic energy of the gas phase: |vd | ≈ −

∗ ξφg

q 2kg φ02

φ(1 − φ) 14

.

(15)

The corrected linear correlation coefficient between the magnitude of the gas-phase velocity and the solid volume fraction is defined as: φ|v| − φ |v| ∗ := p ξφg . p φ0 φ0 hvi000 vi000 ig

(16)

We introduce the correlation coefficient in this form, since it is directly related to the TKE of the gas-phase. We call it corrected, however, since the correlation coefficient between the gas-phase velocity and the solid volume fraction is usually defined regarding the spatial average and not the Favre-average: φ|v| − φ |v| q . ξφg := p φ0 φ0 vi000 vi000

(17)

The difference between hvi000 vi000 ig and vi000 vi000 is not negligible, as will be discussed in

Section 4. As can be seen in Figure 2, the correlation coefficient has large negative values on average. This represents the preference of the gas-phase to pass through regions of low solid volume fraction, since the resistance due to drag is smaller in these regions. On average, it shows weaker correlation at higher solid volume

fractions. In addition, it can also take positive values (Askarishahi, M. and Salehi, M.-S. and Radl, S., 2018), for example during the initial phase of fluidization of a packed bed, where both, the gas-phase velocity and the solid volume fraction decrease in vertical direction in the upper region of the packed bed. The model for the magnitude of the drift velocity (15), tested against the filtered fine grid simulation data in Case 1 and Case 2, is shown in Figure 3. It can be observed that high variance in the solid volume fraction together with pronounced velocity fluctuations in the gas phase are the main sources for drag reduction. The velocity fluctuations of the particle phase have an indirect contribution to the filtered drag force, since the variance of the solid volume fraction depends on the TKE of the solid phase (Schneiderbauer, 2017) (equation (48)). The drift velocity, however, has a preferred orientation (Cloete et al., 2018, 2019; Schneiderbauer, 2018a; Schneiderbauer and Saeedipour, 2018). By comparing Figures 3 and 4, for example, one can conclude that the vertical component, i.e. the component in direction of the gravitational acceleration, contributes most to the magnitude of the drift velocity. Therefore, closure models for each of the individual

15

0

1

-0.2

0.8

-0.4

0.6

-0.6

0.4

-0.8

0.2

2 1.5 1

a)

-1

0

0.2

0.4

0.6

0.8

0

1

0

1

-0.2

0.8

-0.4

0.6

-0.6

0.4

-0.8

0.2

0.5

0

0.2

0.4

0.6

0.8

0

1

0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

2 1.5 1

b)

-1

0

0.2

0.4

0.6

0.8

0

1

0.5

0

0.2

0.4

0.6

0.8

1

0

Figure 2: Correlation coefficients as a function of the normalized filtered solid volume fraction, for ˆ f , a) in Case 1 and b) in Case 2. different dimensionless filter sizes ∆ The symbols correspond to the calculation of the correlation coefficients on the fine grid using the resolved variables v, u and φ, the lines correspond to the predictions by the calculation on the coarse grid using test-filters on the filtered variables hvig , huis and φ.

b)

a) 0

0 -0.05

-0.2

-0.1 -0.4 -0.15 -0.6

-0.8

-0.2

0

0.2

0.4

0.6

0.8

-0.25

1

0

0.2

0.4

0.6

0.8

1

Figure 3: Closure model for the magnitude of the drift velocity |vd | compared to the predictions

ˆ f , as a obtained by filtering the fine grid simulation data for different dimensionless filter sizes ∆ function of the normalized filtered solid volume fraction, a) for Case 1 and b) for Case 2. The symbols correspond to measurements from the filtered fine grid simulation data, the lines correspond to the predictions by employing equation (15).

drift velocity components are required, which may be written in the following form, q ∗ ξφg,ι 2kg,ι φ02 vd,ι = , (18) φ(1 − φ)

16

with the anisotropic components of the TKE of the gas-phase kg,ι defined as 1 kg,ι = hvι000 vι000 ig , 2

(19)

where we do not use Einstein’s summing convention on greek indices. These are essentially the diagonal components of the anisotropic Reynolds-stress tensor. The corrected linear correlation coefficient between the ι-th component of the gas-phase velocity and the solid volume fraction can be defined as: φvι − φ vι ∗ ξφg,ι := p . p φ0 φ0 hvι000 vι000 ig

(20)

Figure 4 demonstrates that the representation of the drift velocity in vertical direction in equation (18) is exact when the correlation coefficient (20) is calculated on the fine grid. The same holds true for the significantly smaller lateral components in wall-bounded fluidized beds, as can be seen in Figure 5. In fact, equation (20) is a relation between the covariance in equation (14) and the variances of the variables, kg,ι and φ02 . In other words, this shows that expression (14) holds true for all components of the drift velocity. In Figure 3, it can be observed that we introduce a small error if we use the same expression for the magnitude of the drift velocity ∗ as defined in (16). This is due to the fact that with the correlation coefficient ξφg

the single components of the drift velocity vector can have different signs in certain flow configurations. In coarse-grid simulations, the correlation coefficients can be calculated dynamically using test-filters, as explained in Section 3.1. The correlation coefficient ξbφg,ι , for example, can be calculated on the coarse grid as: bd \ φhv ι ig − φhvι ig ξbφg,ι = q . q 0 0 d 000 φ φ hvι\ i000 g hvι ig

(21)

In Figure 2, it can be seen that this procedure yields good approximations for the correlation coefficients and that the correlation coefficients have a very low dependence on the filter-size indeed (Schneiderbauer, 2017). The fractional correction functions calculated using the closures for the individual components of the drift velocity (equation (18)) in equation (13) are a good approximation for the filtered drag force reduction in the presence of mesoscale structures, as shown in Figure 6. It can be observed that, in addition to the drift 17

velocity, also the magnitude of the fractional correction function in lateral directions differs from the magnitude in vertical direction (Cloete et al., 2018). Meaning that the ratio between the drift- and slip-velocities is smaller in horizontal- than in vertical direction. a)

b) 0

0

-0.05

-0.2

-0.1 -0.4 -0.15 -0.6

-0.8

-0.2

0

0.2

0.4

0.6

0.8

-0.25

1

0

0.2

0.4

0.6

0.8

1

Figure 4: Closure model for the vertical component of the drift velocity vd,z compared to the predictions obtained by filtering the fine grid simulation data for different dimensionless filter sizes ˆ f , as a function of the normalized filtered solid volume fraction, a) for Case 1 and b) for Case 2. ∆ The symbols correspond to measurements from the filtered fine grid simulation data, the lines correspond to the predictions by employing equation (18).

4. Transport Equations for the Turbulent Kinetic Energies In order to obtain a closure model for the individual components of the drift velocity, we need to solve transport equations or employ constitutive relations for the anisotropic components of the TKE of the gas phase kg,ι , as well as for the variance of the solid volume fraction φ02 . Schneiderbauer (2017) derived transport equations for the variance of the solid volume fraction (constitutive equation written here in terms of the anisotropic components of the TKEs, equation (48)) and the TKE of the solid phase (28). We will follow this procedure to derive a transport equation for the TKE of the gas phase. Multiplying the filtered gas-phase momentum balance equation (7) by the phase-averaged gas-phase velocity hvj ig yields, hvj ig

∂((1 − φ)hvi ig ) ∂ hvj ig + hvj ig ((1 − φ)hvi ig hvk ig ) = (fg,i + Rig ), ∂t ∂xk ρg

18

a)

b)

0.01

0.01

0.005

0.005

0

0

-0.005

-0.005

-0.01

-0.01

-0.015

-0.015

-0.02

0

0.2

0.4

0.6

0.8

-0.02

1

0.01

0.01

0.005

0.005

0

0

-0.005

-0.005

-0.01

-0.01

-0.015

-0.015

-0.02

0

0.2

0.4

0.6

0.8

-0.02

1

0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

Figure 5: Closure models for the lateral components of the drift velocity vd,x and vd,y compared to the predictions obtained by filtering the fine grid simulation data for different dimensionless filter ˆ f , as a function of the normalized filtered solid volume fraction, a) for Case V1 and b) for sizes ∆ Case V2. Note that these cases correspond to three-dimensional wall-bounded simulations whose parameters are summarized in Table 5. The symbols correspond to measurements from the filtered fine grid simulation data, the lines correspond to the predictions by employing equation (18).

with the filtered forces acting on the gas-phase fg,i , and the Reynolds-stress like contributions: Rig

 ∂  000 000 ρg (1 − φ)vi vk . =− ∂xk

Similarly, hvi ig

∂((1 − φ)hvj ig ) ∂ hvi ig + hvi ig ((1 − φ)hvj ig hvk ig ) = (fg,j + Rjg ). ∂t ∂xk ρg

19

a)

1

1

1

0.8

0.8

0.8

0.6

0.6

0.6

0.4

0.4

0.4

0.2

0.2

0.2

0

b)

0

0.2

0.4

0.6

0.8

1

0

0

0.2

0.4

0.6

0.8

1

0

1

1

1

0.8

0.8

0.8

0.6

0.6

0.6

0.4

0.4

0.4

0.2

0.2

0.2

0

0

0.2

0.4

0.6

0.8

1

0

0

0.2

0.4

0.6

0.8

1

0

0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

Figure 6: Fractional correction function H calculated using the drift velocity closure for the filtered drag force, compared to the predictions obtained by filtering the fine grid simulation data for ˆ f , as a function of the normalized filtered solid volume fraction, different dimensionless filter sizes ∆ a) for Case 1 and b) for Case 2. The symbols correspond to measurements from the filtered fine grid simulation data, the lines correspond to the predictions by employing equations (13) and (18).

A combination of these equations yields a transport equation for the gas-phase mean kinetic energy: ∂((1 − φ)hvi ig hvj ig ) ∂ + ((1 − φ)hvi ig hvj ig hvk ig ) ∂t ∂xk hvj ig hvi ig = (fg,i + Rig ) + (fg,j + Rjg ). ρg ρg

(22)

Furthermore, by multiplying the momentum balance equation of the gas phase (2) by the gas-phase velocity vj and averaging, we obtain: ∂ 1 ∂(1 − φ)vi vj + ((1 − φ)vi vj vk ) = (vi fg,j + vj fg,i ). ∂t ∂xk ρg

(23)

Using the definition of the variances (4), and the definition of the triple correlation, as given by Germano (1992), hvi000 vj000 vk000 ig := hvi vj vk ig − hvi ig hvj000 vk000 ig − hvj ig hvi000 vk000 ig − hvk ig hvi000 vj000 ig − hvi ig hvj ig hvk ig ,

20

in (23) and subtracting equation (22) yields: ∂(1 − φ)vi000 vj000 ∂ + ((1 − φ)vi00 vj00 hvk ig ) = (24) ∂t ∂xk   hvj ig hvi ig 1 g g (vi fg,j + vj fg,i ) − (fg,i + Ri ) + (fg,j + Rj ) ρg ρg ρg  ∂  (1 − φ)vi000 vk000 hvj ig + (1 − φ)vj000 vk000 hvi ig + (1 − φ)vi000 vj000 vk000 . − ∂xk

Taking the trace of equation (24) and using the definition of the TKE of the gasphase (5) finally yields a transport equation for the TKE of the gas phase: ∂(1 − φ)kg ∂ + ((1 − φ)kg hvj ig ) = (25) ∂t ∂xj 1 ∂hvi ig 1 ∂ (vi fg,i − hvi ig fg,i ) − (1 − φ)vi000 vj000 − (1 − φ)vi000 vi000 vj000 . ρg ∂xj 2 ∂xj The second term on the right-hand side of this equation is a large-scale shear production term, similar to single-phase turbulence. The last term corresponds to the diffusion of TKE. The total contribution of the work done by the forces acting on the phase is given by: vi fg,i − hvi ig fg,i

  ∂p ∂p = − vi (1 − φ) − hvi ig (1 − φ) ∂xi ∂xi   ∂ g ∂ g Σij − hvi ig Σ + vi ∂xj ∂xj ij   − vi β(vi − ui ) − hvi ig β(vi − ui ) .

Using the zeroth order approximation for the filtered drag force (12), we can express the interfacial work term as (Fox, 2014):     β β − vi β(vi − ui ) − hvi ig β(vi − ui ) = − vi φ(vi − ui ) − hvi ig φ(vi − ui ) φ φ

βe ≈− φvi vi − φvi ui − hvi ig φvi − φui (26) φ βe = φ (hu00i vi000 is − hvi000 vi000 is − hvi000 is (hvi ig − hui is )) . φ Where we define the mixed averages as

hu00i vi000 is = hui vi is − hui is hvi is hvi000 vi000 is = hvi vi is − 2hvi is hvi ig + hvi ig hvi ig , in the framework of Germano averaging (Germano, 1992). 21

Finally, the complete transport equation for the TKE of the gas-phase is given by:

∂(1 − φ)kg ∂ ∂hvi ig + ((1 − φ)kg hvj ig ) = − (1 − φ)vi000 vj000 (27) ∂t ∂xj ∂xj | {z } | {z } convection shear production   1 g ∂ 000 βe  00 000  000 000 000 Σ v + hui vi is − hvi vi is − hvi is (hvi ig − hui is ) − {z } | {z } ρg | ρg ij ∂xj i | {z } exchange+dissipation CIT production −

∂ 0 1 p (1 − φ)vi000 ρg ∂xi | {z }

pressure dilatation+transport

∂ − ∂x | j



dissipation

 1 1 000 g 000 000 000 (1 − φ)vi vi vj − vi Σij . 2 ρg {z } diffusion+transport

The transport equation for ks reads (Schneiderbauer, 2017):

∂ ∂hui is 00 00 βe ∂φks + (φks huj is ) = − φui uj + (hu00i vi000 is − hu00i u00i is ) ∂t ∂xj ∂xj ρ {z } | {z } | {z } |s convection

1 kc ∂ 00 Σ u − ρs ij ∂xj i | {z } dissipation

exchange+dissipation

shear production

−

(28)

1 00 ∂ 0 φu p ρs i ∂xi | {z }

pressure dilatation+transport

−

∂ ∂x | j



 1 00 kc 1 00 00 00 φu u u − u Σ . 2 i i j ρs i ij {z } diffusion+transport

The second term on the left-hand side of the transport equations is the convection term, the terms on the right-hand side are the production due to mean shear, the interphase-exchange term, the viscous dissipation, the pressure dilatation and the diffusion term. Note that equivalent transport equations for the sub-filter scale kinetic energies were introduced by Simonin et al. (He and Simonin, 1993; Simonin, 1991), first in the framework of RANS modelling. All terms also arise in the TKE transport equation in single-phase compressible flows, with the exception of the interfacial work term. The interfacial work term in the transport equation for the TKE of the gas phase contains an additional turbulence production term, which is not present in the transport equation of the TKE of the solid phase. It stems from phase-averaging and, in particular, the approximation for the filtered drag force (12). Thus, turbulent kinetic energy is produced even in the absence of mean shear. This additional turbulence generation mechanism depends on the drift velocity. Therefore, it is zero if the particle-phase is homogeneous at the mesoscale, i.e. if there is no clustering. It is generally referred to as cluster-induced turbulence 22

(CIT) (Capecelatro et al., 2015). This is similar to pseudo-turbulence, where the presence of particles generates turbulence in the gas-phase at the micro-scale (Fox, 2014). Note that the viscous dissipation of the TKE appears as a source term in the filtered granular temperature equation (Capecelatro et al., 2015; Fox, 2014), or in the filtered internal energy equation, analogous to single-phase flow. As discussed before, however, in this work we will not solve a transport equation for the filtered granular temperature (3). While it is valid to assume that the TKEs are in a statistically steady state, it is generally not true that the TKEs are dissipated locally and that transport can be neglected, see Figure 7. The magnitudes of the convective and diffusive terms are on average of the same order of magnitude or of higher order than the viscous dissipation. While for the solid phase, the shear production is the most relevant term, the CIT production in the interfacial work term is on average larger than the shear production in the gas phase. Due to the large density of the solid phase, the pressure-dilatation only plays a minor role. For the gas phase, however, the pressure gradient, which consists mostly of the hydrostatic pressure of the mixture density (10), is much higher in dense gas-particle flows than it would be in single-phase gas flows, since the mixture density is considerably higher compared to the gas-phase density. Thus, similar to high Mach number flows, the pressure dilatation becomes the main TKE dissipation mechanism. In the following, we propose closure models for all terms individually. 4.1. Closure Models for the Turbulent Kinetic Energy Transport Equations 4.1.1. Interfacial Work Closure relationships for the unresolved mixed averages in the interfacial work term (26) are needed in order to evaluate the TKE. These represent the fluid phase properties seen by the particles. We want to relate these mixed averages to the Favre-averages of the relevant phase. For example, the variance of the gas-phase velocity seen by the particles hvi000 vi000 is needs to be related to the turbulent kinetic energy of the gas phase. Therefore, we introduce the spatial average of the square of the gas-phase velocity fluctuations, defined in the framework of Germano averaging

23

a)

b)

106

106

104

104

102

102

10

0

10 0

0.2

0.4

0.6

0.8

1

106

10

0

0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

106

4

10

4

102

102

100

100 0

0.2

0.4

0.6

0.8

1

Figure 7: Budget analysis of all terms in the transport equation for the TKE of both phases scaled by their respective densities as a function of the normalized filtered solid volume fraction, for a ˆ f = 7, a) in Case 1 and b) in Case 2. dimensionless filter size ∆

(Germano, 1992) as: vi000 vi000 := vi vi − 2hvi ig vi + hvi ig hvi ig . It can be expressed as a combination of the phase averaged squares of the velocity fluctuations: vi000 vi000 = φhvi000 vi000 is + (1 − φ)hvi000 vi000 ig . We use the approximation vi000 vi000 − hvi000 vi000 ig = vi000 vi000 −

((1 − φ) − φ0 )vi000 vi000 φ0 vi000 vi000 ≈ , (1 − φ) (1 − φ)

(29)

where we introduced an error in the last step, since spatial averages generally do not fulfill φ = φ. Thus, we can approximate the gas-phase velocity fluctuations seen by the particles by the TKE of the gas-phase and a triple correlation between the 24

solid volume fraction and the gas-phase velocity (Fox, 2014): hvi000 vi000 is ≈ hvi000 vi000 ig +

φ0 vi000 vi000 . φ(1 − φ)

Hence, we can express the mixed average entirely in terms of the TKE of the gasphase, hvi000 vi000 is ≈ 2ξg@s kg , with a coefficient ξg@s given by: ξg@s = 1 +

φ0 vi000 vi000 . 2kg (1 − φ)φ

(30)

Therefore, the triple correlation between the gas-phase velocity and the solid volume fraction needs to be closed. It can be approximated as, φ0 vi000 vi000 ≈ φ0 vi0 vi0 − 2(hvi ig − vi )(φvi − φvi ), where we, again, assumed φ0 = 0. Now, we can use the definition of the triple correlation φ0 vi0 vi0 as a function of the generalized central moments given by Germano (1992): φ0 vi0 vi0 = φvi vi − φ vi0 vi0 − 2vi φ0 vi0 − φvi vi = φvi vi − φvi vi + φvi vi − 2vi (φvi − φvi ) − φvi vi = φvi vi − φvi vi − 2vi (φvi − φvi ). Thus, the triple correlation φ0 vi000 vi000 can be approximated as: φ0 vi000 vi000 ≈ φvi vi − φ vi vi − 2hvi ig (φvi − φvi ).

(31)

Therefore, we can define a ’triple correlation coefficient’, ξφgg :=

φvi vi − φ vi vi − 2hvi ig (φvi − φ vi ) q q , p 000 000 000 000 0 0 φ φ vi vi vi vi

and express the triple correlation as: φ0 vi000 vi000

≈ ξφgg

q

φ02 vi000 vi000 .

(32)

The ’triple correlation coefficient’ ξφgg , relates the triple correlation to the turbulent kinetic energy of the gas-phase and the variance of the solid volume fraction. It is, 25

however, not a linear correlation coefficient by definition. It can have a magnitude of greater than one. Technically, it represents the correlation between the solid volume fraction and the square of the magnitude of the gas-phase velocity corrected by the individual correlation between the solid volume fraction and the gas-phase velocity. From equation (29), it follows that vi000 vi000 = 2kg +

φ0 vi000 vi000 , (1 − φ)

therefore, equation (32) becomes φ0 vi000 vi000

≈

2ξφgg √  ξφgg φ02 1 − (1−φ)

q φ02 kg .

Thus, the correlation coefficient (30) is given by: q ξφgg φ02  ξg@s = 1 + √ . ξφgg φ02 φ(1 − φ) 1 − (1−φ) Finally, the mixed average between the gas- and solid-phase velocity fluctuations, can be expressed by the turbulent kinetic energies of the phases, ∗ hu00i vi000 is ≈ 2ξgs

p ks kg ,

(33)

∗ is the corrected linear correlation coefficient defined as where ξgs

with

v u u u ∗ ξgs = ξgs u1 + t

q ξφgg φ02  √ , ξφgg φ02 φ(1 − φ) 1 − (1−φ)

hui vi is − hui is hvi is ξgs := p 00 00 p 000 000 . hui ui is hvi vi is

It has to be emphasized that Fox (2014) proposed the same model as (33) for the mixed fluid-particle velocity covariance with ξgs = 1, in the case of low particle Reynolds number, dilute gas-particle flow (one-way coupling). In our case of moderately dense gas-particle flow, the correlation coefficient ξgs is dependent on the solid volume fraction and takes the value of approximately 0.5 for a large number of volume fractions, since in the case of four-way coupling the velocities of the two phases 26

are not as strongly correlated as in the case of one-way coupling. This behaviour was also reported previously by Schneiderbauer (2017). As can be seen in Figure 2, in regions of low solid volume fraction the correlation is slightly stronger. In addition, the triple correlation between the gas phase velocity and the solid volume cannot be neglected in moderately-dense gas-particle flow regimes, where the TKE of the gas-phase seen by the particles differs significantly from the TKE of the gas-phase, for example in regions with clusters. The dependence of the triple correlation on the filtered solid volume fraction is shown in Figure 8. For filtered solid volume fractions larger than about φ = 0.15, the TKE of the gas-phase seen by the particles becomes smaller than the actual TKE of the gas-phase, with a minimum at about φ = 0.3. For smaller volume fractions, the TKE of the gas-phase seen by the particles appears to be larger. Therefore, the ’triple-correlation coefficient’ ξg@s is on average larger than unity for small solid volume fractions up to φ = 0.15. As for the other correlation coefficients, the dynamic test-filter approach also yields good estimates for the triple correlation coefficient, as can be seen in Figure 2. a)

b)

0.04

2

0.02

1

10 -3

0

0

-1 -0.02 -2 -0.04

-3

-0.06 -0.08

-4 0

0.2

0.4

0.6

0.8

-5

1

0

0.2

0.4

0.6

0.8

1

Figure 8: Triple correlation between the solid volume fraction and the gas-phase velocity, as a function of the normalized filtered solid volume fraction, compared to the measurements of the ˆ f , a) in Case 1 and b) in filtered fine grid simulation data, for different dimensionless filter sizes ∆ Case 2. The symbols correspond to measurements of the filtered fine grid simulation data, the lines correspond to the predictions made by employing equation (31).

Thus, in addition to the closure for the drift velocity vd,i = hvi000 is , 27

equation (15), we relate the arising mixed averages to the TKEs of the phases: hvi000 vi000 is ≈ 2ξg@s kg , p ∗ hu00i vi000 is ≈ 2ξgs kg ks .

(34)

Substituting the identities (34) and (15) into equation (26) yields an approximation for the interfacial work terms,  p  00 000 00 00 ∗ e e β (hui vi is − hui ui is ) ≈ 2β ξgs kg ks − ks ,

(35)

for the solid phase and

βe (hu00i vi000 is − hvi000 vi000 is − hvi000 is (hvi ig − hui is )) ≈    q ∗  p ξφg ∗ βe 2 ξgs kg ks − ξg@s kg − 2kg φ02 |hvig − huis | , (1 − φ)φ

(36)

for the gas phase. This is a different approach to closing the fluid-particle covariance, for which full transport equations can be solved (Simonin et al., 1993). Note that ∗ the correlation coefficient ξφg is negative on average (see Figure 2), thus, the last

term in equation (36) does indeed represent a positive contribution to the turbulent kinetic energy of the gas phase (CIT). In addition, it is important to note that during the energy transfer between the gas- and the solid phase, a large portion of TKE is dissipated. 4.1.2. Shear Production Considering the other terms in the transport equations, we re-express the shear production in terms of the TKEs of both phases, respectively, as, σg := −

∂hvi ig (1 − φ)vi000 vj000 ≈ 2(1 − φ)kg ∂xj

∂hui is 00 00 φui uj ≈ 2φks σs := − ∂xj where Seijg and Seijs are defined as:

q Seg Seg ,

q Seijs Seijs ,

∂ Seijs := hui is . ∂xj

∂ Seijg := hvi ig , ∂xj

ij

ij

(37)

(38)

(39)

Testing against the filtered fine grid simulation data shows that this approximation is valid, see Figure 9. As discussed in Section 5, in coarse-grid simulations, both, 28

the Reynolds-stresses, as well as the gradients of the phase-averaged velocities are known. Therefore, no additional closure models are needed, and one can simply use the left-hand side of equations (37) and (38). This yields the additional benefit of also being able to incorporate TKE loss due to shear, which is generally not covered by i.e. a turbulent viscosity ansatz. a)

10

b)

2

10

2

10 0

10 0

10 -2

10 -2 0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

Figure 9: Closure models for the shear production of the TKE of both phases σg and σs , respectively, compared to the predictions obtained by filtering the fine grid simulation data for different ˆ f , as a function of the normalized filtered solid volume fraction, a) for dimensionless filter sizes ∆ Case 1 and b) for Case 2. The symbols correspond to measurements from the filtered fine grid simulation data, the lines correspond to the predictions by employing equations (37) and (38).

4.1.3. Viscous Dissipation In addition, we use Prandtl’s mixing length theory (White, 1991), which states that the dissipation takes place over a characteristic turbulent length scale, the mixing length: `mg

p kg . =q g eg e Sij Sij

Thus, we approximate the dissipation as:

3/2

g :=

1 g ∂ 000 kg Σij vi ≈ Cg (1 − φ) ρg ∂xj `mg q = Cg (1 − φ)kg Seijg Seijg ,

1 ∂ 00 s := Σsij u ≈ Cs φks ρs ∂xj i 29

q

Seijs Seijs .

(40) (41)

The constants Cq , q = g, s are not filter-size dependent, but especially Cg was observed to vary by orders of magnitude depending on the considered particle properties. In Table 3, the calculated values of the model constants are summarized for the considered Cases. As can be seen in Figure 7, the ratio between the dissipation of gas-phase TKE through pressure-dilatation versus viscous damping differs by an order of magnitude between the considered Geldart B and A type particles. In Case 1, the contribution of viscous damping to the total TKE dissipation is an order of magnitude smaller than in Case 2. Similarly, it was observed that the value of the constant Cg was an order of magnitude smaller in Case 1 than in Case 2. At the same time, it can be observed that the loss of TKE during the energy transfer between the phases is an order of magnitude higher in Case 1 than in Case 2. Meaning that for the Geldart type B particles more gas-phase TKE is dissipated through drag exchange than for the Geldart type A particles. Even if the drag coefficient is larger in the latter case, the ratio between the gas-phase and the solid-phase TKE is larger in Case 1 leading to a higher loss of TKE. This could be an indication that the nature of the gas-phase TKE dissipation mechanism changes depending on the particle properties. The drift velocity and the TKE of the gas-phase were observed to be larger for larger particle sizes. The viscous dissipation of gas-phase TKE, however, seems to be nearly identical in both Cases and independent of the particle properties (see Figure 7). This would explain why, with larger values of vi000 , the constant in the single-phase LES-like closure model is smaller. On the other hand, the dissipation due to drag exchange is larger for the Geldart type B particles, indicating that there is a change in the main dissipation scale of the flow. Most of the gas-phase TKE is dissipated through drag exchange and not at the viscous dissipation scale of the gas-phase. Note that the viscous dissipation takes place in the boundary layers and wakes of the particles at scales which are not resolved in a fine-grid two-fluid model simulation, where a drag coefficient is employed instead. As can be seen in Figure 10, the models show promising results in a comparison against filtered fine grid simulation data. It has to be emphasized that the values of the constants Cq (q = g, s), which are summarized in Table 3, differ significantly from those given by Schneiderbauer (2017), Cq ≈ 1. This is because Schneiderbauer

(2017) did not consider the pressure-dilatation as an additional dissipation term, but

30

modelled it together with the viscous dissipation. a)

b)

10 1

10 1

10 0

10 0

10 -1

10 -1

10 -2

10 -2

10 -3

10 -3

10 -4

0

0.2

0.4

0.6

0.8

10 -4

1

0

0.2

0.4

0.6

0.8

1

Figure 10: Closure models for the viscous dissipation of the TKE of both phases g and s , respectively, compared to the predictions obtained by filtering the fine grid simulation data for different ˆ f , as a function of the normalized filtered solid volume fraction, a) for dimensionless filter sizes ∆ Case 1 and b) for Case 2. The symbols correspond to measurements from the filtered fine grid simulation data, the lines correspond to the predictions by employing equations (40) and (41).

4.1.4. Diffusion Using the gradient assumption, the diffusion of TKE can be approximated as,     1 1 000 g νtg ∂kg ∂ ∂ 000 000 000 (1 − φ)vi vi vj − vi Σij ≈ (1 − φ) , (42) δg := − ∂xj 2 ρg ∂xj Sctg ∂xj

where the turbulent viscosity is defined as (mixing length model (Igci et al., 2008; White, 1991)): νtg ≈ kg1/2 `mg . In the case of the solid phase: δs := −

∂ ∂xj



1 00 00 00 1 φui ui uj − u00i Σkc ij 2 ρs



≈





∂ks  ∂  φ k q s . ∂xj Scts 2 Ses Ses ∂xj ij ij

(43)

The turbulent Schmidt numbers Sctq are of order 1, see Table 3. As can be seen in Figure 11, these single-phase LES-type closure models are also valid in the considered multiphase framework. Furthermore, it can be seen that the diffusion of the TKE of the gas-phase is important especially in regions with low solid volume fraction, while the diffusion of the TKE of the solid phase becomes important for larger solid volume fractions. 31

a)

b)

10 2

10 2

10 0

10 0

10 -2

10 -2

10 -4

10 -4 0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

Figure 11: Closure models for the diffusion of the TKE of both phases δg and δs , respectively, compared to the predictions obtained by filtering the fine grid simulation data for different dimenˆ f , as a function of the normalized filtered solid volume fraction, a) for Case 1 sionless filter sizes ∆ and b) for Case 2. The symbols correspond to measurements from the filtered fine grid simulation data, the lines correspond to the predictions by employing equations (42) and (43).

4.1.5. Pressure Dilatation Finally, we will discuss a closure model for the term involving the gas-phase pressure gradient, which consists mainly of the pressure-dilatation. As discussed above, the gradient of the pressure fluctuation can be decomposed into its static and dynamic components (Fox, 2014), ∂ 0 ∂ 0 p = ρ0m gi + p˜ , ∂xi ∂xi with the mixture density ρm as defined in (10). Thus, the pressure-dilatation term can be decomposed as: (1 − φ)vi000

∂ 0 ∂ 0 p = (1 − φ)hvi000 ρ0m gi ig + (1 − φ)hvi000 p˜ ig . ∂xi ∂xi

The static term can be re-expressed as: (1 − φ)hvi000 ρ0m gi ig = (1 − φ)(ρs − ρg )gi hφ0 vi000 ig . The dynamic contribution is assumed to be negligible (Fox, 2014). It could, however, have a significant contribution, especially in wall-bounded flows where particles cannot fall freely close to the walls. Therefore, we introduce a constant Cpg and propose the following closure model for the pressure dilatation: −

∂ 0 Cpg 1 (1 − φ)vi000 p ≈− (1 − φ)(ρs − ρg )gi hφ0 vi000 ig . ρg ∂xi ρg 32

Our data suggests that the triple correlation φ0 φ0 vi000 ≈ 0 is negligible, therefore, hφ0 vi000 ig ≈ φ0 vi000 . ∗ Thus, we can use the correlation coefficient ξφg (as defined in Section 3.2.3) in the

closure model: πg := −

∂ 0 Cpg 1 (1 − φ)vi000 p ≈− (1 − φ)(ρs − ρg )gi φ0 vi000 ρg ∂xi ρg q Cpg ∗ ≈− (1 − φ)(ρs − ρg )gi ξφg,i 2kg,i φ02 . ρg

(44)

As can be seen in Figure 12, this approximated model fits well with the predictions obtained by filtering the fine-grid simulation data. Similar to the constant in the dissipation closure, Cpg can be calculated on the coarse grid for each case individually, in the herein considered cases, however, Cpg ≈ 1. Note that the vector of

gravitational acceleration g is negative in vertical direction, since the correlation ∗ coefficient ξφg,i is also negative on average, the pressure-dilatation term is indeed a

dissipative term. Given its large density, the term involving the gas-phase pressure gradient is negligible for the solid phase. The complete transport equations with all closures are summarized in Table 4. The model constants are summarized in Table 3 for both Cases. 5. Transport Equations for the Reynolds-Stresses In fact, equation (24) is a transport equation for the Reynolds-stresses, which can be solved individually for each component. Thus, if we do not take the trace of equation (24), but rather derive the transport equation for the ιι-th component of the variance of the velocity, we obtain the same form of transport equation as (27) for kg,ι , as defined in (19). Using similar closure arguments as above, we can solve a transport equation for the anisotropic components of the TKEs, here written for

33

a)

b) 250

1400 1200

200

1000 800

150

600

100

400 50

200 0

0

0.2

0.4

0.6

0.8

0

1

0

0.2

0.4

0.6

0.8

1

Figure 12: Closure model for the pressure dilatation πg compared to the predictions obtained by ˆ f , as a function of the filtering the fine grid simulation data for different dimensionless filter sizes ∆ normalized filtered solid volume fraction, a) for Case 1 and b) for Case 2. The symbols correspond to measurements from the filtered fine grid simulation data, the lines correspond to the predictions by employing equation (44).

the gas phase (similar for the solid phase): q ∂ ∂(1 − φ)kg,ι g eg + ((1 − φ)kg,ι hvj ig ) = 2(1 − φ)kg,ι Seιj Sιj ∂t ∂xj    q ∗  p ξφg,ι βeι ∗ kg,ι ks,ι − ξg@s,ι kg,ι − 2kg,ι φ02 (hvι ig − huι is ) + 2 ξgs,ι ρg (1 − φ)φ q q C pg g eg ∗ − Cg (1 − φ)kg,ι Seιj Sιj − 2kg,ι φ02 (1 − φ)(ρs − ρg )gι ξφg,ι ρg   ∂  (1 − φ) ∂kg,ι  kg,ι q + . ∂xj Sctg 2 Seg Seg ∂xj ιj

(45)

ιj

The filtered drag coefficient βeι , is the microscopic drag coefficient evaluated from

the filtered velocities in ι-th direction. Instead of solving transport equations for all components of the Reynolds-stress tensor, we suggest calculating the non-diagonal components from the diagonal stresses using correlation coefficients: hvi000 vj000 ig ≈ 2ξij

p kg,i kg,j .

Therefore, no additional closure for the Reynold-stresses, such as for example an (isotropic) Boussinesq hypothesis, is required. The correlation coefficients between the components of the gas-phase velocities show very low negative correlations. They can be determined in coarse-grid simulation by using test-filters as described in Section 3.1. As an example, the correlation coefficient between the two lateral velocities vx and vy is shown in Figure 13. It has to be noted that the low values of 34

Constant

Case 1

Case 2

Cg

0.04

0.6

Cs

0.2

0.4 0.2

† Cφg · Cg

0.3

Cφs · Cs†

1

Cpg Sctg

2

2

Scts

1

1

Table 3: Model constants for the SA-TFM. †

Taken from Schneiderbauer (2017).

the correlation coefficients, which might be fluctuating around zero, could introduce errors in the coarse-grid simulations. Although we did not observe this phenomenon in our short verification study in this work (Section 6), we will further investigate this in a future a posteriori study. a)

b)

1

1

0.5

0.5

0

0

-0.5

-0.5

-1

0

0.2

0.4

0.6

0.8

-1

1

0

0.2

0.4

0.6

0.8

1

Figure 13: Correlation coefficient between the x and y components of the gas-phase velocity, as a ˆ f, function of the normalized filtered solid volume fraction, for different dimensionless filter sizes ∆ a) in Case 1 and b) in Case 2. The symbols correspond to the calculation of the correlation coefficient on the fine grid using the resolved variables v, u and φ, the lines correspond to the predictions by the calculation on the coarse grid using test-filters on the filtered variables hvig , huis and φ.

35

Transport equation for the TKE of the gas phase: q ∂(1 − φ)kg ∂ g eg + ((1 − φ)kg hvj ig ) = 2(1 − φ)kg Seij Sij ∂t ∂xj    q q ∗  p ξφg βe ∗ g eg Sij 2 ξgs kg ks − ξg@s kg − 2kg φ02 |hvig − huis | − Cg (1 − φ)kg Seij + ρg (1 − φ)φ   q Cpg ∂  (1 − φ) kg ∂kg  ∗ q − 2kg,i φ02 + (1 − φ)(ρs − ρg )gi ξφg,i . (46) ρg ∂xj Sctg 2 Seg Seg ∂xj ij

ij

Transport equation for the TKE of the solid phase:

q  ∂φks ∂ 2βe  ∗ p s es + (φks huj is ) = 2φks Seij Sij + ξgs kg ks − ks ∂t ∂xj ρs   q ∂ φ k ∂k s s s es  . q Sij + − Cs φks Seij ∂xj Scts 2 Ses Ses ∂xj ij

(47)

ij

Constitutive relation for the variance of the solid volume fraction (Schneiderbauer, 2017) in terms of the ∗ anisotropic parts of the TKE of the solid phase and the correlation coefficients ξφs,i :

φ0 2 =

 2 p ∂φ ∗ 8 ξφs,i ks,i ∂xi

1/2

∂huk is ks + Cφs Cs ∂xk lms

!2 .

Table 4: Transport equations for the TKE of both phases and constitutive relation for the variance of the solid volume fraction.

36

(48)

6. Verification of the Dynamic Anisotropic SA-TFM In the following we provide a verification of the herein developed dynamic, anisotropic spatially-averaged two-fluid model in the case of 3-dimensional wallbounded fluidized beds of Geldart type A and B particles, respectively. The relevant physical properties of both verification cases, hereafter called Case V1 and Case V2, are summarized in Table 5. Moreover, for the reference fine-grid simulation, we used a grid spacing of 4 · ds and 8 · ds for the Geldart type A and B particles, respectively, in order to resolve all relevant heterogeneous structures (Fullmer and Hrenya, 2016;

Schneiderbauer, 2018a; Schneiderbauer and Saeedipour, 2019; Uddin and Coronella, ˆ g /∆ ˆ f = 6, 8, 12, where 2017). Three different coarse-graining ratios are considered, ∆ ˆ g is the grid spacing of the coarse-grid simulations. As shown by Kong et al. ∆ (2017), the nature of the flow, in particular, the size and structure of clusters and thus, the flow statistics depend strongly on the ratio between the domain size and the characteristic length scale of the mesoscale structures. Therefore, we chose to keep the domain size constant in the considered test cases and to use less grid points in the simulations in order to achieve the different coarse graining ratios. All cases were simulated using the open-source CFD software OpenFOAM. In particular, we added additional turbulence models to the OpenFOAM solver twoPhaseEulerFoam, which solve the transport equations for the Reynolds-stresses and the algebraic equation for the variance of the solid volume fraction, equations (45) - (48). The components of the drift velocity are calculated according to equation (18), the mesoscale interphase force is calculated according to equation (11). The former was implemented as an additional term in the momentum equations, the latter as an additional term in the pressure equation. Thus, the original implementation of the drag terms was not modified. We employed the Wen and Yu (1966) drag law in this a posteriori analysis. Again, following our previous studies, we used a kinetictheory based TFM (Lun et al., 1984; Schneiderbauer et al., 2012) with additional treatment of friction (Schneiderbauer et al., 2012). In particular, Case V1 is the same as in Schneiderbauer and Saeedipour (2019), where the authors discussed the application of an approximate deconvolution model. For the gas-phase we applied a no-slip boundary condition, for the solid phase we used a partial slip boundary condition accounting for the wall friction. For the Reynolds-stresses we employed 37

Property

Case V1

Case V2

φmax

0.6

φfr

0.4

hφid /φmax

0.22

µst i ds ρs

0.28 150 · 10−6 m 2500 kg m−3

75 · 10−6 m

1500 kg m−3

ρg

1.224 kg m−3

µg

1.78 · 10−5 kg m−1 s−1

ut

0.96 m s−1

0.218 m s−1

Fr

626

65

Lch

1.3 · 10−3 m

3 · 10−4 m

Wg

0.63 m s−1

0.07 m s−1

ˆf ∆

1

1

ˆd ∆

125 x 33 x 1000

125 x 33 x 1040

Table 5: Material Properties and Operating Conditions for the Verification Cases (Schneiderbauer, 2018a; Schneiderbauer and Pirker, 2014; Schneiderbauer et al., 2013). Note that these properties correspond to Geldart group type B and A particles (φmax : solids volume fraction at maximum packing conditions; φfr : threshold solids volume fraction for which the frictional stress becomes active; hφid : domain averaged solids volume fraction; µst i : statistic coefficient of internal friction; ds : particle diameter; ρs : particle density; ρg : gas density; µg : gas dynamic viscosity; ut : ter-

minal settling velocity; Fr: particle Froude number; Lch : characteristic length scale of mesoscale structures (Radl and Sundaresan, 2014), Lch = (u2t /g)Fr−2/3 , where g is the magnitude of the ˆ f : dimensionless grid spacing of fine-grid gravitational acceleration; Wg : gas superficial velocity; ∆ ˆ d : dimensionless domain size). simulation; ∆

slip boundary conditions at the walls. The convective terms in the balance equations, as well as, the spatial gradients were discretized using a second-order scheme. The pressure-velocity coupling was solved using the PIMPLE algorithm, where we used linear interpolation to compute the face pressure values. Time advancement was achieved by a variable time-step procedure, where the time step is limited by a maximum Courant number of 0.25. The correlation coefficients were calculated dynamically using the test-filter defined 38

Constant

Case 1

Case 2

Cg

0.1

0.6

Cs

0.2

0.3 0.2

† Cφg · Cg

Cφs · Cs†

0.3

Cpg

0.6

0.5

Sctg

2

2

Scts

1

1

Table 6: Model constants used in the Verification Cases.

in equation (8). We based the implementation of this filter on the OpenFOAM simpleFilter, which interpolates and summates the values of the neighbouring cells. In cells adjacent to the boundary, the values specified at the boundary are taken into account in the summation, alongside with the interpolations of the other neighbouring cells. In this study, we did not use any other adjustment of the test-filters, which still yields reasonable values for the correlation coefficients in the vicinity of the walls. As suggested by Lilly (1992) and others, after the local, dynamic calculation of the correlation coefficients, the values were smoothed using an additional spatial filter. The model parameters were chosen to be constant as listed in Table 6. It has to be emphasized that we had to adjust the value of the model constant for the pressure dilatation, Cpg , in the presence of walls. As discussed before, particles cannot fall freely close to the walls and the pressure gradient will be different from the hydrostatic pressure gradient. In our narrow test cases, the influence of the walls was observed to be particularly large. The simulations were started from a slightly non-uniform settled state. After the initial transient period of bed expansion, the system reached a statistically steady state with persistent temporally and spatially distributed clusters. Snapshots of the flow field were collected at various times after the initialization phase (Milioli et al., 2013; Schneiderbauer, 2017, 2018a). In Figure 14, snapshots of the flow field at time t = 5 s are shown for the different coarse-grid simulations of Case V1. It can be observed that the models correctly 39

predict the bed expansion, as well as the nature of the particle distribution inside the fluidized bed. For coarser grids, however, the transition to the free-board becomes more blurred. This could be due to the dynamic adjustment on the very coarse grids, which leads to smoothed-out correlation coefficients. In addition, the applied regularization yields even more uniform coefficients, thus, local flow variations may not be adequately captured. In particular, this can also be observed in Figure 15 a), which shows the time-averaged axial profile of the solid volume fraction, where h i

depicts the time average over a range of 2 s after the initial transient bed expansion phase, as well as the average over the lateral components. In Figure 15 b), the variance of the solid volume fraction evaluated from the fine-grid simulation is compared to the one evaluated from the coarse-grid simulations by solving equation (48). It has to be noted that the algebraic equation for the solid volume fraction variance seems to predict almost identical values for the different filter-sizes, while in reality the variance will increase with the filter-size. In Figure 14, it can also be observed that the particle clusters become less pronounced for coarser grids. We assume that there is a loss of information through using the simplified algebraic equation in place of a full transport equation for the variance of the solid volume fraction. In fact, the assumption that the volume fraction fluctuations would decay simultaneously with the turbulence in the absence of mean gradients in the flow was used to close the unresoloved terms in the transport equation (Schneiderbauer, 2017), which does not seem to be applicable in the considered cases. This will be investigated in future studies. The anisotropy of the Reynolds-stresses, which was discussed above, as well as in various recent studies (Cloete et al., 2018; Fullmer and Hrenya, 2016; Schneiderbauer and Saeedipour, 2019) is supposed to have a great influence on the macroscopic flow properties. Cloete et al. (2018) suggested that an isotropic treatment of the mesoscale stresses, which leads to an overestimation of the horizontal stresses, yields an overestimation of the bed expansion due to the particles being pushed away from the walls. In Figure 16, it can be observed that the Reynolds-stresses in vertical direction kq,z , q = g, s, are indeed larger than in horizontal direction kq,x . ˆ 6/7 (Milioli In addition, our data shows that the Reynolds-stresses scale with ∆ filter et al., 2013; Schneiderbauer and Saeedipour, 2019). The models for the anisotropic

40

a)

c)

b)

d)

e)

2.5

2.5

2.5

2

2

2

1.5

1.5

1.5

1

1

1

0.5

0.5

0.5

0 -0.07

0

0.07

0 -0.07

0

0.07

0 -0.07

0

0.07

Figure 14: Snapshots of the solid volume fraction of Case V1 at t = 5 s in a slice of the domain: ˆ g = 1∆ ˆ f ; b) fine-grid filtered with a filter of size ∆ ˆ filter = 8∆ ˆ f ; c) coarse-grid with a a) fine-grid, ∆ ˆ g = 6∆ ˆ f ; c) coarse-grid with a grid-size of ∆ ˆ g = 8∆ ˆ f ; c) coarse-grid with a grid-size grid-size of ∆ ˆ g = 12∆ ˆ f. of ∆

components of the Reynolds-stresses (45), seem to correctly reproduce the fine-grid statistics. Note that we used a grid refinement near the walls only in the case of ˆ g = 12, since otherwise the number of cells in horizontal y-direction would have ∆ been less than 5 (the dimensionless domain size is given in Table 5). We suppose that there should be a minimum amount of grid cells in order to reproduce the ˆ g = 8, we used a uniformly spaced grid resulting macroscopic flow correctly. For ∆ in 5 cells in y-direction, which still lead to reasonable macroscopic flow behaviour. The Reynolds-stresses in horizontal x-direction, however, are underpredicted by the simulations on the coarser grids, supporting the hypothesis of Schneiderbauer and Saeedipour (2019) that the grid resolution (number of cells in lateral directions) is insufficient to correctly reproduce the fine-grid statistics. In Figure 16 d), it can especially be observed that the horizontal components of the solid-phase mesoscale stress tensor are reproduced correctly on the coarser grid with wall refinement, and ˆ g = 12, but underesthus, with a larger number of grid cells in lateral direction, ∆

41

b)

a) 2.5

2.5

2

2

1.5

1.5

1

1

0.5

0.5

0

0

0.2

0.4

0.6

0.8

0

1

0

0.02

0.04

0.06

0.08

0.1

ˆg Figure 15: a) Time-averaged axial profile of the solid volume fraction for different grid spacings ∆ and b) time averaged profile of the variance of the solid volume fraction in Case V1. The symbols correspond to the measurements from the fine-grid simulation data, which was filtered using filters ˆ filter , the lines correspond to the measurements from the coarse-grid simulations of different sizes ∆ ˆ g. with different grid spacings ∆ The averaging symbol h i depicts the time average over a range of 2 s after the initial transient

bed expansion phase, as well as the average over the lateral components.

ˆ g = 8, where less grid points were used in timated in the overall finer simulations ∆ horizontal direction. In addition, it can be observed that the vertical components of the Reynolds-stresses are slightly underestimated in the upper regions of the fluidized bed on the coarser grids. This could be a consequence of the gradually more smeared-out mesoscale structures in the higher parts of the fluidized bed, as shown in Figure 14. Similarly, there is a slightly less pronounced transition between the fluidized bed and the free-board for larger filter sizes in Case V2, as can be observed in Figure 17. The solid volume fraction distribution inside this homogeneously bubbling fluidized bed, however, is correctly reproduced, as can be observed in Figure 18 a) as well. Figure 18 b) shows that the variance of the solid volume fraction is overestimated for the smaller coarse-graining ratios. Similar to the discussion of Figure 15 b), the predicted variance of the solid volume fraction is observed to have the same value for all considered filter-sizes. The vertical components of the Reynolds-stress tensors seem to be overestimated in Case V2, as can be seen in Figure 19. The overestimation, however, seems to become less pronounced for larger grid spacings.

42

a)

b)

2.5

2.5

2

2

1.5

1.5

1

1

0.5

0.5

0 10-2

10-1

100

101

0 10-2

102

10-1

c) 2.5

2

2

1.5

1.5

1

1

0.5

0.5

10-1

100

101

102

101

102

d)

2.5

0 10-2

100

101

0 10-2

102

10-1

100

Figure 16: Time-averaged axial profiles of the vertical and horizontal components of the Reynoldsstress tensors in Case V1. The symbols correspond to the measurements from the fine-grid simulation data, which was filtered ˆ filter , the lines correspond to the measurements from the coarse-grid using filters of different sizes ∆ ˆ g . The averaging symbol h i depicts the time average simulations with different grid spacings ∆

over a range of 2 s after the initial transient bed expansion phase, as well as the average over the lateral components.

The lateral components of the Reynolds-stress tensor seem to be underpredicted for all considered grid sizes. Overall, the macroscopic flow properties, as well as the fluctuating components of the solid volume fraction and the velocities could be reproduced to a satisfying extent in both Verification Cases. 7. Conclusion and Outlook In this paper, we validated closure models for the individual components of the drift velocity and, thus, for the filtered drag force, against filtered fine-grid 43

a)

b)

c)

d)

e)

2.5

2.5

2.5

2

2

2

1.5

1.5

1.5

1

1

1

0.5

0.5

0.5

0 -0.07

0

0.07

0 -0.07

0

0.07

0 -0.07

0

0.07

Figure 17: Snapshots of the solid volume fraction of Case V2 at t = 5 s in a slice of the domain: ˆ g = 1∆ ˆ f ; b) fine-grid filtered with a filter of size ∆ ˆ filter = 8∆ ˆ f ; c) coarse-grid with a a) fine-grid, ∆ ˆ g = 6∆ ˆ f ; c) coarse-grid with a grid-size of ∆ ˆ g = 8∆ ˆ f ; c) coarse-grid with a grid-size grid-size of ∆ ˆ g = 12∆ ˆ f. of ∆

simulation data. Highly resolved TFM simulations of Geldart type A and B particles in moderately-dense unbound fluidization were conducted. The drift velocity was expressed as a correlation between the variance of the solid volume fraction and the turbulent kinetic energy of the gas phase. The TKE of the gas phase was observed to be highly anisotropic. Therefore, transport equations for the anisotropic parts of the TKE, in other words, for the Reynolds-stresses, were derived and closure models known from single-phase turbulence modelling were applied to the arising unresolved terms in the transport equation. Thus, it can be concluded that multiphase turbulence models can be viewed analogous to single-phase compressible turbulence models, with the exception of the contribution of the interfacial forces. In fact, gas-phase TKE is generated and sustained via interfacial work if the drift velocity, as well as the slip velocity are non-zero. This additional turbulence generation mechanism is generally referred to as cluster induced turbulence (Capecelatro et al., 2015). The mixed averages in the interfacial work term were related to the TKEs of 44

b)

a) 2.5

2.5

2

2

1.5

1.5

1

1

0.5

0.5

0

0

0.2

0.4

0.6

0.8

0

1

0

0.02

0.04

0.06

0.08

ˆg Figure 18: a) Time-averaged axial profile of the solid volume fraction for different grid spacings ∆ and b) time averaged profile of the variance of the solid volume fraction in Case V2. The symbols correspond to the measurements from the fine-grid simulation data, which was filtered using filters ˆ filter , the lines correspond to the measurements from the coarse-grid simulations of different sizes ∆ ˆ g. with different grid spacings ∆ The averaging symbol h i depicts the time average over a range of 2 s after the initial transient

bed expansion phase, as well as the average over the lateral components.

the phases using correlation coefficients. In contrast to low particle Reynolds number and low mass loading flows, the gas-phase properties seen by the particles were observed to differ significantly from the actual gas-phase properties in moderately dense regimes. Therefore, the triple correlation between the gas-phase velocity and the solid volume fraction, for example, was found to be non-negligible. It was modelled as a function of the turbulent kinetic energy of the gas-phase and the variance of the solid volume fraction scaled by a ’triple correlation coefficient’. Following the scale-similarity theory, we proposed to calculate the correlation coefficients present in the model equations, as well as the model constants dynamically using test-filters. An a posteriori verification of the presented models in two cases of 3-dimensional wall-bounded fluidized beds showed that the dynamic adjustment is applicable even in the presence of walls and that the fine-grid simulation statistics can be reproduced by the presented models to a satisfying extent. No additional wall treatment was applied in the considered cases. We suspect, however, that in riser flows where i.e. the gas-phase velocity profile is crucial, some wall treatment needs to be considered. Future work will focus on more in-depth validations against fine-grid simulations 45

0.1

a)

b)

2.5

2.5

2

2

1.5

1.5

1

1

0.5

0.5

0 10-3

10-2

10-1

100

0 10-3

101

10-2

c) 2.5

2

2

1.5

1.5

1

1

0.5

0.5

10-2

10-1

100

101

100

101

d)

2.5

0 10-3

10-1

100

0 10-3

101

10-2

10-1

Figure 19: Time-averaged axial profiles of the vertical and horizontal components of the Reynoldsstress tensors in Case V2. The symbols correspond to the measurements from the fine-grid simulation data, which was filtered ˆ filter , the lines correspond to the measurements from the coarse-grid using filters of different sizes ∆ ˆ g . The averaging symbol h i depicts the time average simulations with different grid spacings ∆

over a range of 2 s after the initial transient bed expansion phase, as well as the average over the lateral components.

and experimental data. In addition, scalar transport, specifically heat transfer, will be considered in this framework. Acknowledgement The financial support by the Austrian Federal Ministry for Digital and Economic Affairs and the National Foundation for Research, Technology and Development is gratefully acknowledged. The author also wants to acknowledge the financial support from the K1MET center for metallurgical research in Austria (www.k1-met.com). 46

References References Agrawal, K., Loezos, P. N., Syamlal, M., Sundaresan, S., 2001. The role of meso-scale structures in rapid gas-solid flows. J. Fluid Mech. 445, 151–185. Anderson, T. B., Jackson, R., 1967. A Fluid Mechanical Description of Fluidized Beds. Ind. Eng. Chem. Fund. 6 (4), 527–539. Askarishahi, M. and Salehi, M.-S. and Radl, S., 2018. Voidage correction algorithm for unresolved Euler–Lagrange simulations. Comput. Part. Mech. 5 (4), 607–625. Beetstra, R., van der Hoef, M. A., Kuipers, J. A. M., 2007. Drag force of intermediate Reynolds number flow past mono- and bidisperse arrays of spheres. AIChE J. 53 (2), 489–501. Capecelatro, J., Desjardins, O., Fox, R. O., 2014. Numerical study of collisional particle dynamics in cluster-induced turbulence. J. Fluid Mech. 747, R2. Capecelatro, J., Desjardins, O., Fox, R. O., 2015. On fluid–particle dynamics in fully developed cluster-induced turbulence. J. Fluid Mech. 780, 578–635. Capecelatro, J., Desjardins, O., Fox, R. O., 2016a. Strongly coupled fluid-particle flows in vertical channels. I. Reynolds-averaged two-phase turbulence statistics. Phys. Fluids 28, 033306. Capecelatro, J., Desjardins, O., Fox, R. O., 2016b. Strongly coupled fluid-particle flows in vertical channels. II. Turbulence modeling. Phys. Fluids 28, 033307. Capecelatro, J., Desjardins, O., Fox, R. O., 2018. On the transition between turbulence regimes in particle-laden channel flows. J. Fluid Mech. 845, 499–519. Cloete, J. H., Cloete, S., Municchi, F., Radl, S., Amini, S., 2018. Development and verification of anisotropic drag closures for filtered Two Fluid Models. Chem. Eng. Sci. 192, 930–954. Cloete, J. H., Cloete, S., Radl, S., Amini, S., 2019. On the choice of closure complexity in anisotropic drag closures for filtered Two Fluid Models. Chem. Eng. Sci. 207, 379–396. 47

De Wilde, J., 2005. Reformulating and quantifying the generalized added mass in filtered gas-solid flow models. Phys. Fluid 17, 113304. Enwald, H., Peirano, E., Almstedt, A.-E., 1996. Eulerian two-phase flow theory applied to fluidization. Int. J. Multiph. Flow 22, 21–66. Fox, R. O., 2014. On multiphase turbulence models for collisional fluid-particle flows. J. Fluid Mech. 742, 368–424. Fullmer, W. D., Hrenya, C. M., 2016. Quantitative assessment of fine-grid kinetictheory-based predictions of mean-slip in unbounded fluidization. AIChE J. 62 (1), 11–17. Février, P., Simonin, O., Squires, K. D., 2005. Partitioning of particle velocities in gas–solid turbulent flows into a continuous field and a spatially uncorrelated random distribution: theoretical formalism and numerical study. J. Fluid Mech. 533, 1–46. Germano, M., 1992. Turbulence: the filtering approach. J. Fluid Mech. 238, 325–336. Germano, M., Piomelli, U., Moin, P., Cabot, W. H., 1991. A dynamic subgrid-scale eddy viscosity model. Phys. Fluids A 3 (7), 1760–1765. Haworth, D. C., Pope, S. B., 1986. A generalized Langevin model for turbulent flows. The Physics of Fluids 29 (2), 387–405. He, J., Simonin, O., 1993. Non-equilibrium prediction of the particle-phase stress tensor in vertical pneumatic conveying. Gas-Solid Flows, ASME FED 166, 253– 263. Igci, Y., Andrews, A. T., Sundaresan, S., Pannala, S., O’Brien, T., 2008. Filtered two-fluid models for fluidized gas-particle suspensions. AIChE J. 54 (6), 1431– 1448. Igci, Y., Sundaresan, S., 2011. Constitutive Models for Filtered Two-Fluid Models of Fluidized Gas–Particle Flows. Ind. Eng. Chem. Res. 50 (23), 13190–13201. Ishii, M., 1975. Thermo-Fluid Dynamic Theory of Two-Phase Flow. Direction des Etudes et Recherches d’Electricité de France, Eyrolles, Paris. 48

Issangya, A., Bai, D., Bi, H., Lim, K., Zhu, J., Grace, J., 1999. Suspension densities in a high-density circulating fluidized bed riser. Chem. Eng. Sci. 54 (22), 5451– 5460. Issangya, A. S., Grace, J. R., Bai, D., Zhu, J., 2000. Further measurements of flow dynamics in a high-density circulating fluidized bed riser. Powder Technol. 111 (1), 104–113. Kong, B., Fox, R. O., Feng, H., Capecelatro, J., Patel, R., Desjardins, O., Fox, R. O., 2017. Euler–euler anisotropic gaussian mesoscale simulation of homogeneous cluster-induced gas–particle turbulence. AIChE J. 63 (7), 2630–2643. Lilly, D. K., 1992. A proposed modification of the Germano subgrid-scale closure method. Phys. Fluids A 4 (3), 633–635. Lun, C. K. K., Savage, S. B., Jeffrey, D. J., Chepurniy, N., 1984. Kinetic theories for granular flow: Inelastic particles in Couette flow and slightly inelastic particles in a general flowfield. J. Fluid Mech. 140, 223–256. Milioli, C. C., Milioli, F. E., Holloway, W., Agrawal, K., Sundaresan, S., 2013. Filtered Two-Fluid Models of Fluidized Gas-Particle Flows: New Constitutive Relations. AIChE J. 59 (9), 3265–3275. Moreau, M., Simonin, O., Bédat, B., 2010. Development of Gas-Particle Euler-Euler LES Approach: A Priori Analysis of Particle Sub-Grid Models in Homogeneous Isotropic Turbulence. Flow Turbul. Combust. 84, 295–324. Ozel, A., Fede, P., Simonin, O., 2013. Development of filtered Euler-Euler two-phase model for circulating fluidised bed: High resolution simulation, formulation and a priori analyses. Int. J. Multiph. Flow 55, 43–63. Ozel, A., Gu, Y., Milioli, C. C., Kolehmainen, J., Sundaresan, S., 2017. Towards filtered drag force model for non-cohesive and cohesive particle-gas flows. Phys. Fluids 29, 103308. Panday, R., Shadle, L., Shahnam, M., Cocco, R., Issangya, A., Spenik, J., Ludlow, J., Gopalan, B., Shaffer, F., Syamlal, M., Guenther, C., Karri, S., Knowlton, 49

T., 2014. Challenge problem: 1. Model validation of circulating fluidized beds. Powder Technol. 258, 370–391. Parmentier, J.-F., Simonin, O., Delsart, O., 2012. A functional subgrid drift velocity model for filtered drag prediction in dense fluidized bed. AIChE J. 58 (4), 1084– 1098. Pope, S. B., 2000. Turbulent Flows. Cambridge University Press. Pärssinen, J., Zhu, J.-X., 2001a. Particle velocity and flow development in a long and high-flux circulating fluidized bed riser. Chem. Eng. Sci. 56 (18), 5295–5303. Pärssinen, J. H., Zhu, J.-X., 2001b. Axial and radial solids distribution in a long and high-flux CFB riser. AIChE J. 47 (10), 2197–2205. Radl, S., Sundaresan, S., 2014. A drag model for filtered Euler–Lagrange simulations of clustered gas–particle suspensions. Chem. Eng. Sci. 117, 416–425. Saeedipour, M., Vincent, S., Pirker, S., 2019. Large eddy simulation of turbulent interfacial flows using Approximate Deconvolution Model. Int. J. Multiph. Flow, 286–299. Sarkar, A., Milioli, F. E., Ozarkar, S., Li, T., Sun, X., Sundaresan, S., 2016. Filtered sub-grid constitutive models for fluidized gas-particle flows constructed from 3-D simulations. Chem. Eng. Sci. 152, 443–456. Schneiderbauer, S., 2017. A Spatially-Averaged Two-Fluid Model for Dense LargeScale Gas-Solid Flows. AIChE J. 63 (8), 3544–3562. Schneiderbauer, S., 2018a. Validation Study on Spatially Averaged Two-Fluid Model for Gas-Solid Flows: I. A Priori Analsis of Wall Bounded Flows. AIChE J. 64 (5), 1591–1605. Schneiderbauer, S., 2018b. Validation study on spatially averaged two-fluid model for gas-solid flows: II. Application to risers and fluidized beds. AIChE J. 64 (5), 1606–1617.

50

Schneiderbauer, S., Aigner, A., Pirker, S., 2012. A comprehensive frictional-kinetic model for gas-particle flows: analysis of fluidized and moving bed regimes. Chem. Eng. Sci. 80, 279–292. Schneiderbauer, S., Pirker, S., 2014. Filtered and heterogeneity based subgrid modifications for gas-solid drag and solids stresses in bubbling fluidized beds. AIChE J. 60(3), 839–854. Schneiderbauer, S., Puttinger, S., Pirker, S., 2013. Comparative analysis of subgrid drag modifications for dense gas-particle flows in bubbling fluidized beds. AIChE J. 59, 4077–4099. Schneiderbauer, S., Saeedipour, M., 2018. Approximate deconvolution model for the simulation of turbulent gas-solid flows: An a priori analysis. Phys. Fluids 30, 023301. Schneiderbauer, S., Saeedipour, M., 2019. Numerical simulation of turbulent gas–solid flow using an approximate deconvolution model. Int. J. Multiph. Flow 114, 287–302. Simonin, O., 1991. Second-moment prediction of dispersed phase turbulence in particle-laden flows. In: Proc. 8th Int. Symp. on Turbulent Shear Flows (University of Munich). Vol. 1. pp. 741–6. Simonin, O., Deutsch, E., Minier, J.-P., 1993. Eulerian Prediction of the Fluid/Particle Correlated Motion in Turbulent Two-Phase Flows. Appl. Sci. Res. 51, 275–283. Stolz, S., Adams, N. A., 1999. An approximate deconvolution procedure for largeeddy simulation. Phys. Fluids 11 (7), 1699–1701. Uddin, M. H., Coronella, C. J., 2017. Effects of grid size on predictions of bed expansion in bubbling fluidized beds of Geldart B particles: A generalized rule for a grid-independent solution of TFM simulations. Particuology 34, 61–69. van Wachem, B. G. M., Schouten, J. C., van den Bleek, C. M., Krishna, R., Sinclair, J. L., 2001. Comparative analysis of CFD models of dense gas–solid systems. AIChE J. 47 (5), 1035–1051. 51

Weber, J. M., Bobek, M. M., Breault, R. W., Mei, J. S., Shadle, L. J., 2018. Investigation of core-annular flow in an industrial scale circulating fluidized bed riser with electrical capacitance volume tomography (ECVT). Powder Technol. 327, 524–535. Wen, C. Y., Yu, Y. H., 1966. Mechanics of Fluidization. Chem. Eng. Prog., Symp. Ser. No. 62, 100. White, F. M., 1991. Viscous fluid flow. The McGraw-Hill, Inc., New York. Zhang, D. Z., VanderHeyden, W., 2002. The effects of mesoscale structures on the macroscopic momentum equations for two-phase flows. Int. J. Multiph. Flow 28 (5), 805–822. Zhu, H., Zhu, J., Li, G., Li, F., 2008. Detailed measurements of flow structure inside a dense gas–solids fluidized bed. Powder Technol. 180 (3), 339–349.

52

Author Statement Stefanie Rauchenzauner: Methododlogy, Software, Validation, Formal analysis, Writing – Original draft Simon Schneiderbauer: Conceptuation, Methodology, Software, Writing – Review & Editing, Supervision, Project administration, Funding acquisition Declaration of interests The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

53