# Machine Learning Approaches to Close the Filtered Two-Fluid Model for Gas-Solid Flows: Models for Subgrid Drag Force and Solid Phase Stress

Baptiste Hardy, Stefanie Rauchenzauner, Pascal Fede, Simon Schneiderbauer, Olivier Simonin, Sankaran Sundaresan, Ali Ozel

► To cite this version:

Baptiste Hardy, Stefanie Rauchenzauner, Pascal Fede, Simon Schneiderbauer, Olivier Simonin, et al. Machine Learning Approaches to Close the Filtered Two-Fluid Model for Gas-Solid Flows: Models for Subgrid Drag Force and Solid Phase Stress. Industrial and engineering chemistry research, 2024, 63 (18), pp.8383-8400. 10.1021/acs.iecr.3c04652. hal-04660088

HAL Id: hal-04660088

https://ut3-toulouseinp.hal.science/hal-04660088v1

Submitted on 23 Jul 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Machine Learning Approaches to Close the Filtered Two-Fluid Model for Gas–Solid Flows: Models for Subgrid Drag Force and Solid Phase Stress

Published as part of Industrial & Engineering Chemistry Research virtual special issue "Multiscale Modeling and Artificial Intelligence for Multiphase Flow Science".

Baptiste Hardy,* Stefanie Rauchenzauner, Pascal Fede, Simon Schneiderbauer, Olivier Simonin, Sankaran Sundaresan, and Ali Ozel

<figure><img src="image_3.png" /></figure>

Cite This: *Ind. Eng. Chem. Res.* 2024, 63, 8383–8400

<figure><img src="image_4.png" /></figure>

Read Online

## ACCESS |

<figure><img src="image_5.png" /></figure>

Metrics & More

<figure><img src="image_6.png" /></figure>

Article Recommendations

**ABSTRACT:** Gas-particle flows are commonly simulated through a two-fluid model at the industrial scale. However, these simulations need a very fine grid to have accurate flow predictions, which is prohibitively demanding in terms of computational resources. To circumvent this problem, the filtered two-fluid model has been developed, where the large-scale flow field is numerically resolved and small-scale fluctuations are accounted for through subgrid-scale modeling. In this study, we have performed fine-grid two-fluid simulations of dilute gas-particle flows in periodic domains and applied explicit filtering to generate data sets. Then, these data sets have been used to develop artificial neural network (ANN) models for closures such as the filtered drag force and solid phase stress for the filtered two-fluid model. The set of input variables for the subgrid drag force ANN model that has been found previously to work well for dense flow regimes is found to work as well for the dilute regime. In

addition, we present a Galilean invariant tensor basis neural network (TBNN) model for the filtered solid phase stress, which can nicely capture the anisotropic nature of the solid phase stress arising from subgrid-scale velocity fluctuations. Finally, the predictions provided by this new TBNN model are compared to those obtained from a simple eddy-viscosity ANN model.

<figure><img src="image_7.png" /></figure>

## INTRODUCTION

Gas-particle flows arise in technological applications, e.g., fluidized and circulating fluidized bed reactors, and in nature, e.g., dust storms. There is much interest in studying the characteristics of these flows via mathematical modeling with complementary computer-aided simulations. Reliable modeling and simulations can aid in the design, retrofit, and troubleshooting of industrial processes. As industrial-scale fluidized beds contain trillions of particles, it is impractical to analyze the flow behavior by following the motion of individual particles. In contrast, two-fluid models (TFMs),¹,² which treat the fluid and particle phases as interpenetrating continua, are more viable to analyze and simulate the flows in industrial-scale applications. They have been useful in the analysis of the onset of instabilities in fluidization, and the emergence of inhomogeneous flow structures (e.g., see refs ³–⁵). The TFM can readily be solved numerically using commercial codes (e.g., ANSYS Fluent) and open-source simulation platforms (e.g., MFIX, OpenFOAM).

It is now well-known that fluidized and circulating fluidized beds manifest structures that span a wide range of length and time scales. The scale of the spatial structures can range from a few particle diameters to the size of the vessel, which can be as large as 10³–10⁵ particle diameters. The macroscale flow structures, termed coherent flow structures, can have a large effect on the overall flow hydrodynamics in the device. At the same time, mesoscale structures (such as streamers, clusters, and small bubble-like voids) are also important as they affect the emergence of the macroscale structures. As a result, accurate simulations of TFM equations often require fine spatial resolution down to the scale of a few particle

Received: December 29, 2023

Revised: March 18, 2024

Accepted: March 20, 2024

Published: April 24, 2024

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>
diameters⁶⁻⁸ which, in turn, require very small time steps as well. Such highly resolved simulations of industrial-scale processes are not feasible due to the high computational demands.⁹ This consideration led to the development of filtered two-fluid models (fTFMs)⁷,⁸,¹⁰⁻¹² where the TFM model equations are filtered by a convolution kernel as in the development of Large Eddy Simulation (LES) equations for turbulent flows by averaging the Navier–Stokes equations.¹³

The fTFM equations contain several terms representing the consequences of subfilter-scale fluctuations on the spatiotemporal evolution of the filtered variables. Similar to the LES modeling, the importance of subgrid terms in the fTFM equations could be studied by a priori tests. It has been shown through the budget analysis of the filtered solid momentum balance generated by filtering fine-grid simulations⁷,⁸,¹¹ that the correction to the fluid-particle drag force is of principal importance for gas-particle flows with high mass loading of particles because of the subfilter-scale inhomogeneous distribution of particles. The solid phase stress associated with the subfilter-scale particle velocity fluctuations is of secondary importance,⁸ while all the other corrections are essentially negligible. As a result, the literature on fTFM model development has focused primarily on correction to the fluid-particle drag force and, to a lesser extent, on the (filtered) solid phase stress.

The drag correction models accounting for the effects of unresolved drag due to particle clustering at mesoscale in the literature can be classified as follows. The Energy Minimization Multi-Scale (EMMS) model¹⁴⁻¹⁶ describes subgrid structures through a heterogeneous index, which is used to estimate the effective drag force. In the framework of fTFM, the explicit correlations were proposed by Igci et al.,¹⁰ Milioli et al.,¹⁷ Sarkar et al.¹⁸ for the filtered drag in terms of the filtered variables and the filter size. Parmentier et al.⁷ argued that in the presence of subfilter-scale (aka subgrid-scale, as the filter size is usually the same as the grid size in coarse simulations of the fTFM equations) inhomogeneities, the average gas velocity seen by the particles is not the same as the filtered gas velocity and expressed the drag force correction in terms of a subgrid quantity known as the drift velocity. Algebraic models for the drift velocity have been proposed in several studies.⁷,⁸,¹⁹ Rauchenzauner and Schneiderbauer²⁰ expressed the drift velocity in terms of the subgrid turbulent kinetic energy of the gas and the scalar variance of the particle volume fraction, which were determined by solving corresponding dynamic transport equations. In a recent study, Hardy et al.²¹ found that the drift velocity could be expressed in terms of the scalar variance of the particle volume fraction, with the same model applying to all filter sizes.

Several research groups have applied machine learning (ML) techniques to arrive at models for the filtered drag force. Jiang et al.²² developed the transport equation for the drift velocity and performed a budget analysis of the terms in the developed equation to analyze their importance. They concluded that an algebraic model for the filtered drift velocity would be sufficient for dense fluidized beds. The algebraic model relates the filtered drift velocity to the filtered gas phase pressure gradient and solid volume fraction and the difference between the filtered gas and solid phase velocities (referred to as the slip velocity). These variables are taken as physics-inspired inputs to an artificial neural network (ANN) model (specifically, a multilayer perceptron, (MLP)) for the drift flux (which is the product of filtered drift velocity and the filtered solid volume

fraction), for a given filter size and physical properties of the gas and particles. The drift flux was then used to predict the filtered drag force, as illustrated by Parmentier et al.⁷ Jiang et al.²³ extended the analysis to include the filter sizes and the Reynolds number based on the terminal settling velocity of the particle as additional inputs so that the ANN model can be used for a wide range of fluidized bed applications. Interestingly, Jiang et al.²² concluded that a good correlation could be obtained with their ANN model only when the output variable was chosen to be the drift flux; in contrast, their ANN model performed poorly when the drag correction or drift velocity was used as the output variable. In a very similar context, Zhang et al.²⁴ found that an ANN model where the hidden layers included a combination of convolutional layers and fully connected layers revealed better predictions for a priori tests of the filtered drag. Zhu et al.²⁵ recently developed data-driven models for the filtered drag force, filtered heat transfer coefficient, and filtered reaction rate in gas-solid flows and performed a posteriori validation analyses for turbulent bubbling fluidized beds and risers. All these studies considered systems without interparticle forces; Tausendschön et al.²⁶ report that the strength of the drag correction is affected by interparticle forces and proposed an ANN model including the Bond number as a measure of the attractive interparticle forces. In the Eulerian–Lagrangian approach, Lu et al.²⁷ developed a filtered drag force ANN model from fine-grid CFD-DEM simulations of dense fluidized beds and they coupled this new drag correction model with the MFiX software for large-scale simulations. For a more extensive review of ML-based modeling efforts in multiphase flow reactors, the reader is referred to the recent study by Zhu et al.²⁸

In spite of the progress in the application of ML methods for formulating subgrid drag models for gas-solid flows without the interparticle forces, unanswered questions remain. The first objective of the present study is to address the following questions:

* Why have prior studies found better predictions with some output variables (namely, the drift flux) but not others?

* Is this observation related to the underlying physics of the problem or the use of suboptimal neural networks?

* Will the input variables for the ANN model identified by Jiang et al.²³ using data sets generated through simulations of dense fluidized beds be sufficient to model the filtered drift flux in dilute flows such as those encountered in risers?

The last question has been addressed by the development of neural network modeling of fTFM closures using two-dimensional fine-grid TFM simulations of turbulent riser flows and performing a posteriori fTFM simulations with developed models.²⁵,²⁹ Different from these studies, we focus on data-driven modeling of the filtered drag force in the dilute regime using fine-grid, three-dimensional, unbounded fluidization simulations and further study the effects of physical parameters on the model development. Besides, it must be stressed that, in this study, we characterize as dilute gas-solid mixtures found in the lower range of the moderately dense regime³⁰ where particle-particle collisions start to contribute significantly to the solid phase microscopic stress model.

As noted above, the filtered solid phase stress associated with the subfilter particle velocity variations is the second most
important correction. Although the budget analysis found this stress term to be secondary, it could play a role in correctly capturing the smaller-scale structures resolved in the fTFM simulations, which in turn could affect the emergence of the macroscale structures in some flow problems. Ozel et al.,⁸ Milioli et al.,¹⁷ Schneiderbauer³¹ proposed an isotropic model for the stress which resembled the Smagorinsky model for the stress in single-phase turbulent flows, and advanced explicit functional models for the pressure and effective viscosity associated with subfilter-scale fluctuations in the velocities of both phases. The importance of anisotropy has been discussed by several researchers.¹⁸,³²⁻³⁴ Rauchenzauner and Schneider-bauer³⁴ have proposed a dynamic multiphase turbulence model for coarse-grid simulations, which includes transport equations for the scalar variance of the solid volume fraction and the individual components of the turbulent kinetic energies of both phases, requiring additional closure models. These extra transport equations improve predictions but add to the computational cost. Thus, as a second objective of the present study, we explore the use of tensor-based neural network models, which have found use in single-phase turbulent flows,³⁵ to constitute the filtered solid phase stress. The predictions offered by this new Galilean invariant model are also compared with simple eddy-viscosity ANN models, similar to earlier proposals in the literature.³⁶ The proposed Galilean-invariant tensor-based neural network model is able to find the components of the filtered solid phase stress tensor, and we demonstrate that it can capture the anisotropic particle phase stress stemming from the subgrid velocity fluctuations. Finally, a sample data set and Python ML model source codes are placed in the GitHub repository (https://github.com/bahardy/fTFM_ANN_modeling.git) for broader use.

The rest of the article is organized as follows. First, the filtered two-fluid model equations and the subgrid terms to be modeled are briefly recalled. Then, the flow configuration and the filtering procedure used to generate data sets are detailed. Subsequently, the Artificial Neural Network architecture chosen to predict the filtered drag force is presented, and the results obtained by this model are discussed. Finally, we introduce the Tensor Basis Neural Network architecture for the particle phase stress; we present the predictions yielded by this more advanced model and compare them with a simple eddy-viscosity approach. The review ends with a summary of the achievements of this work and suggestions for future research.

**Filtered Two-Fluid Model.** As detailed in the work of Igci et al.¹⁰ and others,⁷,⁸,³⁷ the filtering of the mass and momentum balance equations of the “micro-scale” TFM leads to the following set of equations for the gas and solid phases:

$$ \frac{\partial}{\partial t}(\rho_g \bar{\phi}_g) + \nabla \cdot (\rho_g \bar{\phi}_g \tilde{\mathbf{u}}_g) = 0 \qquad (1) $$

$$ \frac{\partial}{\partial t}(\rho_g \bar{\phi}_g) + \nabla \cdot (\rho_g \bar{\phi}_g \tilde{\mathbf{u}}_g) = 0 \qquad (2) $$

$$ \frac{\partial}{\partial t}(\rho_g \bar{\phi}_g \tilde{\mathbf{u}}_g) + \nabla \cdot (\rho_g \bar{\phi}_g \tilde{\mathbf{u}}_g \tilde{\mathbf{u}}_g) = -\nabla \cdot \sigma_g - \bar{\phi}_g \nabla p_g - \nabla \cdot \bar{\Sigma}_g^d - \bar{\mathbf{I}}_{gs} + \rho_g \bar{\phi}_g \mathbf{g} \qquad (3) $$

$$ \begin{aligned} \frac{\partial}{\partial t}(\rho_s \bar{\phi}_s \tilde{\mathbf{u}}_s) + \nabla \cdot (\rho_s \bar{\phi}_s \tilde{\mathbf{u}}_s \tilde{\mathbf{u}}_s) &= -\nabla \cdot \sigma_s - \bar{\phi}_s \nabla p_g - \nabla \cdot \bar{\Sigma}_s + \bar{\mathbf{I}}_{gs} + \rho_s \bar{\phi}_s \mathbf{g} \end{aligned} \qquad (4) $$

Here, $\rho_g$ and $\rho_s$ are the gas and solid phase densities, respectively, $p_g$ is the gas phase pressure, and $\mathbf{g}$ is the gravitational acceleration. The filtered volume fractions for gas ($k = g$) and solid ($k = s$) phases are defined as

$$ \bar{\phi}_k(\mathbf{x}, t) = \int_V \phi_k(\mathbf{y}, t) G(\mathbf{y} - \mathbf{x}) \, d\mathbf{y}, \quad k = g, s \qquad (5) $$

where $\phi_k(\mathbf{x}, t)$ is the volume fraction for each phase given by the “micro-scale” TFM, $G(\mathbf{r})$ is the filter convolution kernel satisfying $\int_V G(\mathbf{r}) \, d\mathbf{r} = 1$. Similar to LES of compressible flows,³⁸ the phase velocities are filtered through the Favre-averaging as

$$ \tilde{\mathbf{u}}_k(\mathbf{x}, t) = \frac{1}{\bar{\phi}_k} \int_V \phi_k(\mathbf{y}, t) \mathbf{u}_k(\mathbf{y}, t) G(\mathbf{y} - \mathbf{x}) \, d\mathbf{y} = \frac{\overline{\phi_k \mathbf{u}_k}}{\bar{\phi}_k}, \quad k = g, s \qquad (6) $$

The filtered quantity for each phase, denoted with a bar, in eq 1 is defined as

$$ \bar{Q}_k(\mathbf{x}, t) = \int_V Q_k(\mathbf{y}, t) G(\mathbf{y} - \mathbf{x}) \, d\mathbf{y}, \quad k = g, s \qquad (7) $$

where $Q_k(\mathbf{y}, t)$ is a quantity for each phase in the “micro-scale” TFM. The filtered and its fluctuating quantities are described as

$$ Q_k(\mathbf{x}, t) = Q_k(\mathbf{x}, t) + Q'_k(\mathbf{x}, t) \qquad (8) $$

$\bar{\Sigma}_g^d$ and $\bar{\Sigma}_s$ correspond to the deviatoric part of the filtered gas phase microscale stresses and to the filtered solid phase microscale stresses. $\bar{\mathbf{I}}_{gs}$ is the filtered gas-solid interphase momentum exchange term; in gas-particle systems, it is principally only the drag force. The explicit expressions for the filtered microscale stress tensors can be found in earlier works.⁷,⁸,¹⁰

The filtered drag force term is defined as

$$ \bar{\mathbf{I}}_{gs} = \overline{\beta(\mathbf{u}_g - \mathbf{u}_s)} \qquad (9) $$

where $\beta$ is the microscopic drag coefficient (also defined as $\beta = \frac{\rho_g \bar{\phi}_g}{\tau_p}$, with $\tau_p$ the particle relaxation time⁷,⁸), here computed from the drag law of Wen and Yu.³⁹

In the literature, it is a very common practice¹¹,¹⁴,¹⁵,²⁶,³⁷,⁴⁰⁻⁴² to account for the drag correction required in the fTFM by introducing an effective drag coefficient $\beta_e$, namely

$$ \bar{\mathbf{I}}_{gs} \approx \beta_e (\tilde{\mathbf{u}}_g - \tilde{\mathbf{u}}_s) \qquad (10) $$

where $\beta_e$ has to be determined for fTFM simulations. The relation between the effective drag coefficient and the microscopic drag coefficient computed from the filtered quantities, here noted $\tilde{\beta}$, is usually expressed in terms of a drag correction factor $H_d$¹⁵ defined as

$$ H_d = \frac{\beta_e}{\tilde{\beta}} \qquad (11) $$

Numerous studies have sought to improve the functional description of $H_d$.⁷,⁸,¹⁰,¹⁷,¹⁸,⁴³ It must be stressed that eqs 10 and 11 implicitly assume that the required drag correction (i.e.,
the subgrid drag force term) is aligned with the filtered slip velocity and that the drag correction factor is isotropic. (Note: the terms subgrid and subfilter are used interchangeably as in fTFM simulations the filter size is usually taken to be the grid size.) As noted by Tausendschön et al.,²⁶ a more general model would be

$$ \bar{\mathbf{I}}_{\text{gs}} = \mathbf{H}_{\text{d}}\tilde{\beta}(\tilde{\mathbf{u}}_{\text{g}} - \tilde{\mathbf{u}}_{\text{s}}) \quad (12) $$

where $\mathbf{H}_{\text{d}}$ is a second-order tensor. Nevertheless, it has been verified⁷,⁸,⁴⁴ that eq 11 is adequate and that, to a very good approximation, the filtered drag force can be written as

$$ \bar{\mathbf{I}}_{\text{gs}} = \tilde{\beta}(\tilde{\mathbf{u}}_{\text{g}} + \mathbf{v}_{\text{d}} - \tilde{\mathbf{u}}_{\text{s}}) \quad (13) $$

where $\mathbf{v}_{\text{d}}$ is the so-called drift velocity, defined as

$$ \bar{\phi}_{\text{s}} \mathbf{v}_{\text{d}} = \bar{\phi}_{\text{s}} \bar{\mathbf{u}}_{\text{g}} - \bar{\phi}_{\text{s}} \tilde{\mathbf{u}}_{\text{g}} \quad (14) $$

In this paper, the product $\bar{\phi}_{\text{s}}\mathbf{v}_{\text{d}}$ will be referred to as the “drift flux”. The filtered drag and the drift flux can be explicitly linked as follows

$$ \bar{\mathbf{I}}_{\text{gs}} = \frac{\rho_{\text{s}} \bar{\phi}_{\text{s}}}{\tilde{\tau}_{\text{p}}} (\tilde{\mathbf{u}}_{\text{g}} - \tilde{\mathbf{u}}_{\text{s}}) + \frac{\rho_{\text{s}}}{\tilde{\tau}_{\text{p}}} \bar{\phi}_{\text{s}} \mathbf{v}_{\text{d}} \quad (15) $$

where $\tilde{\tau}_{\text{p}}$ is the particle relaxation time computed from filtered quantities. Most functional models proposed in the literature considered $\mathbf{v}_{\text{d}}$ to be aligned with the resolved slip velocity⁷,⁸ (which eventually boils down to an expression of the form given by eq 10) though, Ozel et al.⁴⁵ suggested that the drift velocity has an additional, albeit small, explicit dependence on the subgrid variance of the solid volume fraction. The subgrid variance of the solid volume fraction also enters the dynamic drift velocity model of Rauchenzauner and Schneiderbauer⁴⁶ while Hardy et al.²¹ recently proposed an explicit model to deduce the drift velocity from the subgrid variance, independently of the filter size. The strong alignment between the drift velocity and the resolved slip velocity was recently confirmed by Dabbagh and Schneiderbauer⁴⁷ based on the PDF of the relative angle between the two quantities.

The second most important subgrid contribution in the fTFM comes from the filtered solid phase (or mesoscale) stresses, defined as

$$ \sigma_{\text{s}} = \bar{\rho}_{\text{s}} \bar{\phi}_{\text{s}} \bar{\mathbf{u}}_{\text{s}} - \bar{\rho}_{\text{s}} \bar{\phi}_{\text{s}} \tilde{\mathbf{u}}_{\text{s}} \tilde{\mathbf{u}}_{\text{s}} \quad (16) $$

While the drag correction has been found to saturate as the filter size increases, the mesoscale stress term grows monotonically with the filter size.⁸ Its modeling is therefore becoming increasingly important, as the mesh resolution is lowered for very large-scale simulations. Prior studies⁶,¹⁰,⁴⁸ have concluded that the filtered solid phase microscale stresses ($\bar{\Sigma}_{\text{s}}$) described by the kinetic theory of granular flows are much weaker than the mesoscale stresses, even for moderately large filter sizes (typically filter sizes larger than 15–20 particle diameters).

Finally, the filtered pressure gradient term can be decomposed as the sum of a resolved and a subgrid term as

$$ \bar{\phi}_{\text{k}} \bar{\nabla}_{\bar{p}_{\text{g}}} = \bar{\phi}_{\text{k}} \nabla_{\bar{p}_{\text{g}}} + \Phi_{\text{k}}^{\text{sgs}}, \quad k = g, s \quad (17) $$

Some authors³⁷,⁴⁹ proposed to model $\Phi_{\text{k}}^{\text{sgs}}$ as an added mass term, but most studies concluded that this term was small with respect to the filtered drag and mesoscale stress terms and that it is sufficient to retain only the resolved part of the pressure gradient, i.e., $\bar{\phi}_{\text{k}} \nabla_{\bar{p}_{\text{g}}} \approx \bar{\phi}_{\text{k}} \nabla_{\bar{p}_{\text{g}}}$.

**Flow Configurations and Generating Data Set through TFM Simulations.** A prerequisite for the training and a priori validation of ANN models for the filtered drag and filtered stresses is the generation of a database of fine-grid TFM simulation results covering a wide range of physical parameters. As noted earlier, Jiang et al.²³ performed fine-grid simulations of dense fluidized beds. Their data sets had only a sparse representation of dilute flow conditions. Therefore, their findings apply to flow conditions with particle volume fractions exceeding 0.1 and become less accurate under more dilute conditions. One of the goals of this study is to examine if the input variables required for the ANN model for drag correction are any different for dilute flow conditions. For that reason, we have performed a number of dilute flow simulations using the two(multi)-fluid solver neptune_cfd.⁵⁰

The computational domain is triperiodic with the gravitational acceleration acting along the z-direction. A source term mimicking an external pressure gradient is added to the gas and solid phase momentum transport equations to compensate for the weight of the mixture against gravitational acceleration. The list of simulated cases and their physical parameters (particle diameter $d_p$ and density $\rho_s$) are summarized in Table 1. The following quantities are fixed through all simulations:

<table><caption>**Table 1.** Physical and Flow Parameters of the Fine-Grid TFM Simulations of a Gas–Solid Unbounded Fluidized Bed in a Tri-Periodic Domain</caption><thead><tr><th></th><th>d<sub>p</sub> (μm)</th><th>ρ<sub>s</sub> (kg m<sup>-3</sup>)</th><th>Re<sub>p</sub></th><th>Fr<sub>p</sub></th><th>Fr<sub>Δ</sub></th></tr></thead><tbody><tr><td>case 1</td><td>75</td><td>1500</td><td>1.10</td><td>65.33</td><td>20.41</td></tr><tr><td>case 2</td><td>90</td><td>1500</td><td>1.79</td><td>101.09</td><td>20.41</td></tr><tr><td>case 3</td><td>100</td><td>1500</td><td>2.37</td><td>128.37</td><td>20.41</td></tr><tr><td>case 4</td><td>75</td><td>3000</td><td>2.05</td><td>228.50</td><td>20.41</td></tr><tr><td>case 5</td><td>150</td><td>2500</td><td>9.88</td><td>663.01</td><td>76.03</td></tr><tr><td>case 6</td><td>150</td><td>1800</td><td>7.63</td><td>395.71</td><td>53.90</td></tr><tr><td>case 7</td><td>180</td><td>1600</td><td>10.68</td><td>448.16</td><td>58.56</td></tr><tr><td>case 8</td><td>130</td><td>1800</td><td>5.40</td><td>304.23</td><td>45.23</td></tr><tr><td>case 9</td><td>180</td><td>2500</td><td>14.98</td><td>882.90</td><td>92.03</td></tr><tr><td>case 10</td><td>120</td><td>2000</td><td>4.83</td><td>309.96</td><td>45.80</td></tr></tbody></table>

the gas density $\rho_g = 1.2 \text{ kg m}^{-3}$, the gas dynamic viscosity $\mu_g = 1.8 \times 10^{-5} \text{ Pa s}$, the particle restitution coefficient $e_c = 0.9$, and the mean solid volume fraction $\langle \phi_s \rangle = 0.05$. We also report the particle Reynolds number based on terminal settling velocity $Re_p = \frac{\rho_g U_d d_p}{\mu_g}$, the particle Froude number $Fr_p = \frac{U_t^2}{g d_p}$, and the Froude number based on the grid size $Fr_Δ = \frac{U_t^2}{g Δ}$. The terminal settling velocity $U_t$ was estimated from Schiller and Naumann drag law.⁵¹ The mesh size $Δ$ is uniform (and can be derived from $Fr_Δ$ in Table 1), with 640 cells along the vertical direction and 160 cells along each horizontal direction.

In the remainder of this study, Case 1 will be considered as the reference case, as it corresponds to typical conditions for the fluidization of Geldart A particles with air. It must be emphasized that this case has already been studied extensively in the literature. A snapshot of the solid volume fraction field from Case 1 is shown in Figure 1, highlighting the formation of typical elongated clusters.

**Filtering of Fine-Grid TFM data.** Filtered and subgrid terms in the fTFM equations can be computed by applying an explicit filter on the fine-grid simulation data.¹⁰ In the present work, we use a box (or top-hat) filter $\bar{G}$ in the physical space:
<figure><img src="image_10.png" /><figcaption>Figure 1. Instantaneous solid volume fraction $\phi_s$ in a TFM simulation of a triperiodic fluidized bed in Case 1 (see Table 1).</figcaption></figure>

$$ \bar{G}(\mathbf{r}) = \begin{cases} \frac{1}{\bar{\Delta}^3}, & |r_i| < 0.5\bar{\Delta}, \quad i = x, y, z \\ 0, & \text{otherwise} \end{cases} \qquad (18) $$

where $\bar{\Delta}$ is the filter width. By doing so, the filtered drag, the drift velocity (and drift flux), and the filtered solid phase stresses have been computed from eqs 13, 14, and 16, respectively, for a range of filter-to-grid size ratios: $\bar{\Delta}/\bar{\Delta} = [2, 4, 6, 8, 10, 12, 16]$.

To address the first objective, we build new ANN-based models for the filtered drag force and solid stresses from a set of well-chosen markers. The selection of those markers should be inspired by the underlying physics and by the many studies that have sought explicit coarser-grained models for fTFM. The case of the filtered drag force and solid stresses will be addressed separately later in more detail, but a few key quantities can already be identified. The filtered slip velocity

$$ \tilde{\mathbf{u}}_{\text{slip}} = \tilde{\mathbf{u}}_{\text{g}} - \tilde{\mathbf{u}}_{\text{s}} \qquad (19) $$

appears in numerous functional models next to the filtered solid volume fraction $\bar{\phi}_s$ and the filter size $\bar{\Delta}$ (see Cloete et al.⁵² for a comparative study of different existing models). It has been more recently observed by Jiang et al.²² and confirmed by later studies²³,²⁴,⁵³,⁵⁴ that the addition of the filtered pressure gradient acting against the gravitational acceleration $\frac{\partial \bar{\phi}_s}{\partial z}$ as an input to their filtered drag ANN model dramatically improved their results. Based on that, Jiang et al.⁵⁵ formulated an explicit drift velocity model with an additional dependence on the filtered pressure gradient. Besides, recent studies⁸,²¹,⁴⁵,⁴⁶ identified the subgrid variance of the solid volume fraction as a potentially good marker for the filtered drag. The underlying idea is that the drift velocity (and therefore the drag correction) originates from a heterogeneous distribution of the particles at the subgrid-scale and that the subgrid variance of the particle volume fraction (henceforth

simply termed SV) is a good measure of this level of heterogeneity. In this study, the SV will be defined as

$$ \bar{\phi}_s^2 = \bar{\phi}_s^3 - \bar{\phi}_s^2 \qquad (20) $$

The SV was previously introduced for filtered drag force modeling by Schneiderbauer¹¹ where it appears in the expression of their drag correction factor. Their model additionally depends on the subgrid correlated kinetic energy of the solid phase, $k_s^{\text{sgs}}$, defined as

$$ k_s^{\text{sgs}} = \frac{1}{2} \text{tr}(\sigma_s) \qquad (21) $$

As for the filtered solid stresses modeling, single-phase turbulence models⁵⁶,⁵⁷ and previous functional modeling efforts for gas-solid flows⁸,¹⁷,³¹ suggest that the deviatoric part of the filtered rate-of-deformation tensor, defined as

$$ \tilde{\mathbf{S}}_s = \frac{1}{2}(\nabla \tilde{\mathbf{u}}_s + \nabla \tilde{\mathbf{u}}_s^\text{T}) - \frac{1}{3} \nabla \cdot \tilde{\mathbf{u}}_s \mathbf{I} \qquad (22) $$

should definitely play a role. The rotation-rate tensor

$$ \bar{\mathbf{R}}_s = \frac{1}{2}(\nabla \tilde{\mathbf{u}}_s - \nabla \tilde{\mathbf{u}}_s^\text{T}) \qquad (23) $$

will also be involved in the ANN modeling of the filtered stresses later on. Overall, the gradient of the phase-filtered velocities and of the filtered solid volume fraction $\nabla \bar{\phi}_s$ have been extracted from our fine-grid TFM simulation results as they contain nonlocal information that can potentially improve the description of the filtered drag and filtered solid stresses.

## NEURAL NETWORK MODELING OF FILTERED DRAG FORCE

The first study to exploit the neural network approach for modeling the filtered drag (FD) force was proposed by Jiang et al.²² These authors first developed the transport equation of the drift velocity. They used the fine-grid TFM simulation results of a bubbling fluidized bed with Geldart-A type particles to assess the importance of unclosed terms in the transport equation. The transport equation was then simplified to the algebraic model, which was closed by using a 3-marker model using the filtered solid volume fraction $\bar{\phi}_s$, the filtered slip velocity $\tilde{\mathbf{u}}_{\text{slip},z}$ and the gas phase gradient $\frac{\partial \bar{\phi}_s}{\partial z}$ in the gravitational acceleration direction. Instead of proposing an explicit functional form for the 3-marker model, they used an MLP made of three hidden layers of 128, 32, and 8 nodes to close the component of the drift flux aligned with gravitational acceleration, i.e. $\bar{\phi}_s v_{d,z}$. The filtered drag force was then deduced from an explicit relation derived from eq 13. These authors observed that the inclusion of the filtered gas phase pressure gradient significantly improved the degree of correlation of their model with the exact filtered drag. Different neural networks were trained for each filter size so that the filter size was not considered as a distinct marker at that point. They were unable to achieve a similar level of correlation when the filtered drag force or the drift velocity was set as an output variable of their ANN, instead of the drift flux. In a later study,²³ the same group verified that a filter-size dependent ANN model could be applied to large-scale simulations. In addition, they extended the range of application of their filtered drag force model by training their ANN on a wide range of physical parameters and by using appropriate scaling for the
input markers. They also concluded that the particle Reynolds number or the Archimedes number should come as an additional marker to account for variations in the physical properties. The prediction improvement offered by the addition of the filtered gas phase pressure gradient was also confirmed by Ouyang et al.⁵⁴ through interpretable ML metrics. Zhang et al.²⁴ proposed a convolutional neural network (CNN) architecture to predict the filtered drag in periodic unbounded gas-solid flows. Information of neighboring grid points was therefore inherently included in their filtered drag force model by the structure of the network, and the authors report increased performance with respect to MLP models or explicit dynamic models. They also concluded that the filtered gas phase pressure gradient improves the predictions of their model but to a lesser extent than with the MLP architecture since information from neighboring cells is already embedded in the model. It is however expected that such a CNN model will be more computationally intensive and therefore less tractable for practical fTFM simulations at an industrial scale. Zhu et al.⁵³ compared a classical MLP and a gradient boosting framework to predict the filtered drag and integrated these ML models into fTFM simulations of bubbling and turbulent fluidized beds. They subsequently validated their ANN-assisted fTFM simulations against experimental results of a bubbling fluidized bed and found satisfactory agreement between the two approaches. Very recently, Tausendschön et al.²⁶ used Eulerian-Lagrangian simulation results to propose distinct ANN models for the different components of the drift flux. These authors found that these anisotropic drag models lead to better results than the isotropic counterpart, as previously noticed by Cloete et al.³²,⁵⁸ for explicit models. These authors also added the Bond number as an additional marker to account for cohesive effects in gas-solid flows.

**Artificial Neural Network (ANN) Architecture for Filtered Drag Force.** To address whether the many different markers used in the previous NN models represent the underlying physics or the use of suboptimal neural networks, we developed our ANN architecture, which is a feedforward MLP similar to Jiang et al.²² This neural network architecture is sketched in Figure 2. It consists of one input layer, some hidden layers, and one output layer. The input layer

<figure><img src="image_11.png" /><figcaption>**Figure 2.** Multilayer perceptron (MLP) architecture for filtered drag force modeling.</figcaption></figure>

corresponds to the physical markers selected to predict the target quantity. Hidden layers are made of a number of nodes or neurons. Because the network is densely connected, each node takes as input all of the nodes in the previous layer. The outcome of a single neuron i within layer n is described by

$$z_i^{(n)} = f(\mathbf{w}^T \mathbf{z}^{(n-1)} + w_0) \quad (24)$$

where $\mathbf{z}^{(n-1)}$ is the output vector of layer $n-1$, $\mathbf{w}$ is the weight vector of the node, $w_0$ is the bias and $f$ is the activation function. The number of hidden layers and the number of nodes per layer are two hyperparameters of the problem that are discussed further below. The output layer is made of a single node whose value (the target quantity $y$) should allow us to compute the filtered drag in a straightforward manner, i.e., using an explicit model.

Using the data sets generated for dilute systems in the present study, to that end, we investigated the following strategies:

1. develop an ANN model for the filtered drag force (FD-ANN): the target quantity $y$ is the filtered drag $\bar{\mathbf{I}}_{gs}$

2. develop an ANN model for the drift velocity (DV-ANN): the target quantity is the drift velocity $v_d$ and one uses an explicit expression similar to eq 13 to compute the filtered drag force, and

3. develop an ANN model for the drift flux (DF-ANN): the target quantity is the drift flux $\bar{\phi}_s v_d$ and one uses an explicit expression similar to eq 15 to compute the filtered drag force.

In what follows, the architecture of the ANN for filtered drag force prediction is kept unchanged in order to compare the capabilities of the aforementioned strategies for the same level of complexity. The current ANN is made of 3 hidden layers with 128, 32, and 8 nodes, respectively, using a ReLU (Rectified Linear Unit) activation function. We tested the ANN structures of 4 hidden layers with 64, 32, 16, and 8 nodes and even a deeper structure of 5 layers with 32, 32, 32, 16 and 8 nodes. These structures provided a very similar root-mean-squared error (RMSE) for the filtered drag force predictions. Therefore, we preferred to use the ANN structure of 3 hidden layers with 128, 32, and 8 nodes as a reference.

The output layer has a linear activation function to return the regression result. The loss function used to train the network is the mean absolute error (MAE), defined as

$$\text{MAE} = \frac{1}{N} \sum_{i}^{N} |\hat{y}_i - y_i| \quad (25)$$

where $N$ is the number of data points used in a training batch. The data set for the reference case is made of $1.792 \times 10^6$ entries (accounting for the 7 filter widths). Training and testing procedures use subsets corresponding to 80 and 20% of the data set, respectively. Among the training data set, 20% of the data are further preserved for validation and to prevent overfitting during the learning process.

The prediction accuracy of the developed ANN models will be assessed using the coefficient of determination $R^2$ defined as

$$R^2 = 1 - \frac{\sum_i (y_i - \hat{y}_i)^2}{\sum_i (y_i - \langle y \rangle)^2} \quad (26)$$

where $y_i$ is the exact value of the target quantity for the ith observation (specifically, the target quantity obtained by
filtering the fine-grid TFM simulation results), $\hat{y}_i$ is the corresponding ANN model output value and $\langle y \rangle$ is the mean value of $y$ over the data set. The quality of the model can also be quantified by the Probability Density Function (PDF) of the relative error $e$, defined as

$$e_i = \frac{(y_i - \hat{y}_i)}{y_i} \qquad (27)$$

## ■ RESULTS FROM THE ANN MODEL ANALYSIS OF THE FILTERED DRAG FORCE

We begin this section by comparing the filtered drag force results of Jiang et al.²² and Rauchenzauner and Schneiderbauer⁴⁶ for dense fluidized beds. As these authors employed different sets of input variables to model the filtered drag force, the comparison described below illustrates the nonuniqueness of the choices of input variables to model the filtered drag force. We then turn our attention to ANN models for dilute flows based on the simulations performed in our study, where we consider the quality of predictions afforded by three different combinations of input and output variables in the ANN models.

**Dense Bubbling Fluidized Bed: Different Choices of Input Variables To Estimate the Drag Correction.** Jiang et al. developed an ANN model for the filtered drag force in dense, bubbling fluidized beds with average solid volume fractions in the range of 0.25–0.4.²³ The authors trained the ANN on about 20 different combinations of particle properties in the Geldart-A and A/B groups. In order to obtain a universally applicable model, the ANN model is based on dimensionless input and output variables. For a specified gas-particle system, the ANN model requires the following four input quantities (and hence their model is referred to as the 4-marker model): the filter size, the filtered solid volume fraction, and the filtered slip velocity, and the component of the filtered gas phase pressure gradient in the gravitational acceleration direction, all of which are available in an fTFM model simulation. They found that this model can be made applicable to different gas-particle systems by including either the particle Reynolds number based on the terminal velocity

$Re_p$ or the Archimedes number $Ar = \frac{(\rho_s - \rho_g)\rho_g d_p^3 g}{\mu_g^2}$ as an additional input variable. They reported a model that employed the Reynolds number as the extra input, which we consider here. Their drift flux ANN model takes the form:

$$\frac{\bar{\phi}_i v_{d,z}}{\phi_{s,max} U_t} = f\left(\frac{\bar{\phi}_i}{\phi_{s,max}}, \frac{\tilde{u}_{slip,z}}{U_t}, \frac{\bar{\Delta}}{L_c}, \frac{1}{\rho_g} \frac{\partial \bar{p}_g}{\partial z}, Re_p\right) \qquad (28)$$

where the characteristic length scale $L_c$ is set equal to $d_p Fr_p^{1/3}$, as suggested by Radl and Sundaresan.⁵⁹

Rauchenzauner and Schneiderbauer⁴⁶ performed fine-grid simulations of a bubbling fluidized bed of Geldart-A type particles, filtered their simulation results, and used them to validate a functional model for drift velocity, which used the filtered solid volume fraction, the SV, and the kinetic energy associated with the subgrid particle phase velocity fluctuations (which they referred to as the turbulent kinetic energy, TKE). As SV and TKE are not available in an fTFM simulation, additional transport equations must be solved to track their spatiotemporal evolution. Nevertheless, it is interesting to see that Rauchenzauner and Schneiderbauer⁴⁶ and Jiang et al.²³

employed different sets of input variables to model the drift velocity. To test whether the data set generated by Rauchenzauner and Schneiderbauer⁴⁶ could have been captured by the ANN model supplied by Jiang et al.,²³ we compared the predictions of the Jiang et al.²³ model with the data generated by Rauchenzauner and Schneiderbauer.⁴⁶ (Specifically, we tested the results from case 2 of Rauchenzauner and Schneiderbauer⁴⁶ which closely corresponded to case 2 of Jiang et al.²³) Figure 3 reveals a fairly

<figure><img src="image_12.png" /><figcaption>Figure 3. Drift velocity predicted by the ANN of Jiang et al.,²³ compared to the filtered fine-grid simulation data produced by a different research group⁴⁶ for a dimensionless filter size $\bar{\Delta}^* = 9$. The filter-size was made dimensionless using the characteristic length-scale $L_c = d_p Fr_p^{1/3}$.</figcaption></figure>

good correlation, suggesting that the drift velocity could be modeled by either combination of input variables. The combination of inputs suggested by Jiang et al.²³ is perhaps advantageous as it does not require the simulation of additional transport equations. It also suggests that SV and TKE can be estimated in terms of the local quantities employed as input variables by Jiang et al.;²³ i.e., the transport equations for SV and TKE can be approximated by a local-equilibrium approximation. The filtered slip velocity and the filtered gas phase pressure gradient in the Jiang et al.²³ model appear in the transport equations for SV and TKE (even when they are simplified with a local-equilibrium approximation).

**New ANN Filtered Drag Force Models for Dilute-to-Moderately Dense Flows.** Model for the Reference Case. The MLP architecture described in Figure 2 was used to predict the filtered drag force for the Reference Case (Case 1 in Table 1), involving Geldart A-type particles in the dilute regime. We explored different choices of input quantities (markers) to the network to see if conclusions drawn in the dense regime apply to the dilute case as well. The present analysis focuses on the vertical component of the filtered drag force given its primary importance in fluidization.

We start by taking the filtered drag force itself as the output quantity of the network (FD-ANN). Figure 4 compares the predictions of the 3- and 4-marker FD-ANN models, respectively, defined by

$$\bar{I}_{gs,z} = f(\bar{\phi}_s, \bar{\Delta}, \tilde{u}_{slip,z}) \qquad (29)$$

and

$$\bar{I}_{gs,z} = f\left(\bar{\phi}_s, \bar{\Delta}, \tilde{u}_{slip,z}, \frac{\partial \bar{p}_g}{\partial z}\right) \qquad (30)$$
<figure><img src="image_13.png" /><figcaption>Figure 4. Assessment of the FD-ANN model predictions: (a) 3-marker model ($\bar{\phi}_s$, $\bar{\Delta}$, $\tilde{u}_{lip,z}$); (b) 4-marker model ($\bar{\phi}_s$, $\bar{\Delta}$, $\tilde{u}_{slip,z}$, $\frac{\partial \bar{\phi}_g}{\partial z}$).</figcaption></figure>

In line with Jiang et al.²² findings in the dense regime, we observe that the inclusion of the filtered gas phase pressure gradient in the markers of the ANN dramatically improves the accuracy of the model prediction in the dilute flow regime as well. Zhang et al.²⁴ came to the same conclusion with an MLP architecture, though the authors report lower correlation coefficients, even when the filtered gas phase pressure gradient is added. Jiang et al.²² found that the FD-ANN model performed poorly for dense flows, which is clearly not the case for dilute flows.

We show in Figures 5 and 6 the predictive capability of the ANN model when the target quantity is the drift flux or the drift velocity, respectively. It can be observed that the vertical component of the drift flux $\bar{\phi}_s v_{d,z}$ is more accurately captured by the ANN than the drift velocity itself, either with the 3- or 4-marker model. However, when the result of the ANN is inserted into the explicit relations given by eqs 13 and 15 to estimate the filtered drag force, the 3-marker DV-ANN and DF-ANN models prove to be considerably inferior to the 4- marker models. Both 4-marker DV-ANN and DF-ANN models perform equally well in dilute flows; in contrast, Jiang et al.²² found that DF-ANN was much better than DV-ANN in the dense flow regime. Taken together, the 4-marker DF-ANN model appears to be better suited than DV-ANN and FD-ANN models for both dilute and dense flows.

It can be inferred from Figure 7 that the 4-marker ANN model described by

$$ \bar{\phi}_s v_{d,z} = f(\bar{\phi}_s, \bar{\Delta}, \tilde{u}_{slip,z}, \bar{\phi}_s^2) \quad (31) $$

is capable of predicting the drift flux in the vertical direction with quite high accuracy ($R^2 = 0.94$). However, the quality of the model for the filtered drag force decreases significantly ($R^2 = 0.78$) as previously observed with the 3-marker DF-ANN and DV-ANN models. We can therefore conclude that the vertical component of the filtered gas phase pressure gradient performs better as the fourth marker than the SV for this reference case.

Zhang et al.²⁴ argued that the flow information on the neighboring grid cells was crucial in predicting the local filtered drag, which is inherently provided by their CNN. Yet, if one

aims to build explicit models inspired by the machine learning approach, we should identify which differential quantities add the most information. To that end, we also tested 4-marker DF-ANN models where the filtered gas phase pressure gradient in the vertical direction is replaced by the vertical component of the filtered solid volume fraction gradient or the filtered solid phase velocity gradient. These attempts yielded slightly better results than those of the 3-marker model, without achieving the same predictive capacity as the gas pressure gradient-based 4-marker model. Therefore, the filtered gas phase pressure gradient appears to be the most promising fourth marker in both the dense and dilute regimes.

Figure 8 shows the PDF of the relative error in the filtered drag force predictions. The different ANN approaches (FD, DV, DF) lead to very similar distributions of the modeling error in the 3-marker case, although the DF-ANN model displays a slightly narrower distribution. In the 4-marker case, the PDF curves of the different ANN models are not distinguishable and are symmetric. The scalar variance-based 4-marker model leads to a narrower distribution of the error than the 3-marker model, without reaching the level of accuracy of the filtered gas phase pressure gradient-based 4- marker model, as discussed above. Based on all these considerations, we consider only the DF-ANN model for the rest of this article.

**Generalized ANN Model.** For practical use in large-scale simulations, it is desirable that the filtered drag force ANN model can be generalized to a wide range of physical parameters. To that end, the input markers of the ANN must be made dimensionless using proper characteristic scales, and the ANN should be trained with different data sets to cover a large range of these dimensionless input quantities. This approach has been investigated by Jiang et al.²³ for bubbling to turbulent fluidized beds, but the corresponding range of solid volume fractions was confined to the fairly dense regime. To take the analysis further, simulation results from Cases 1 to 9 in Table 1 have been used to train a generalized drift flux ANN model. Assuming that the 4-marker DF-ANN model examined in the previous section is sufficient to capture
<figure><img src="image_14.png" /><figcaption>Figure 5. Assessment of the DV-ANN model predictions on (a, b) the drift velocity and (c, d) the filtered drag force using explicit eq 15. (a, c) 3-marker model ($\bar{\phi}_s$, $\bar{\Delta}$, $\tilde{u}_{slip,z}$). (b, d) 4-marker model ($\bar{\phi}_s$, $\bar{\Delta}$, $\tilde{u}_{slip,z}$, $\frac{\partial \bar{p}_g}{\partial z}$).</figcaption></figure>

the filtered drag force for constant physical parameters, a more general DF-ANN model can be expressed as

$$ \bar{\phi}_s v_{d,z} = f \left( \bar{\phi}_s, \tilde{u}_{slip,z}, \bar{\Delta}, \frac{\partial \bar{p}_g}{\partial z}, d_p, \rho_s, \epsilon_c, \rho_g, \mu_g \right) \quad (32) $$

The restitution coefficient is kept constant in our TFM simulations so that its influence could not be assessed. The scaling proposed by Jiang et al.²³ in the dense regime to reduce the number of independent variables is given by eq 28. They could not discriminate between the three different definitions of $L_c$ usually found in the literature: $L_{c,I} = d_p Fr_p^{1/3}$ and $L_{c,II} = \frac{U_t^2}{g} = d_p Fr_p$ and $L_{c,III} = d_p$, and the authors set $L_c = L_{c,I}$.

We started by adopting the same characteristic length scale in our analysis. Figure 9 shows the predictions of the generalized DF-ANN model described in eq 28 with $L_c = L_{c,I}$. The training has been performed using 80% of the data points of Cases 1–9 while the remaining 20% have been preserved for testing purposes (as shown in Figure 9). It can be concluded that the generalized ANN model is able to predict the filtered drag

force with very good accuracy ($R^2 = 0.94$). The model is also tested against data from Case 10 (Figure 9), which was not used for the training process. It can be observed that the performance level of the model is preserved, which means that the scaling proposed in eq 28 is also appropriate in the dilute regime and that the data set used for the training was sufficiently populated in order to build a robust model.

Then, considering that the characteristic length scale could be more generally written as $L_c = d_p \mathcal{F}(Fr_p, Re_p)$ (where $\mathcal{F}$ is some function of the Froude and Reynolds numbers), we studied the case where $L_c$ is taken equal to the particle diameter ($L_c = L_{c,III}$) and the Froude number is added as an additional marker in eq 28, i.e.,

$$ \frac{\bar{\phi}_s v_{d,z}}{\phi_{s,max} U_t} = f \left( \frac{\bar{\phi}_s}{\phi_{s,max}}, \frac{\tilde{u}_{slip,z}}{U_t}, \frac{\bar{\Delta}}{d_p}, \frac{1}{\rho_s g} \frac{\partial \bar{p}_g}{\partial z}, \mathrm{Re}_p, Fr_p \right) \quad (33) $$

It is shown in Figure 10 that the predictions of this model are not superior to the one presented above with the presumed definition of the characteristic length scale, and the choice $L_c = L_{c,I}$ appears as the best option for now, also in the dilute
<figure><img src="image_15.png" /><figcaption>Figure 6. Assessment of the DF-ANN model predictions on (a, b) the drift flux and (c, d) the filtered drag force using explicit eq 15. (a, c) 3-marker model ($\bar{\phi}_s$, $\bar{\Delta}$, $\tilde{u}_{slip,z}$). (b, d) 4-marker model ($\bar{\phi}_s$, $\bar{\Delta}$, $\tilde{u}_{slip,z}$, $\frac{\partial \bar{\phi}_g}{\partial z}$).</figcaption></figure>

regime. Nevertheless, we note that we used 9 different cases to train our model while Jiang et al.²³ trained their model with more than 20 different cases. Thus, enlarging the dilute regime data sets in future studies to include more cases spanning a wider range of parameters and solid volume fractions appears to be the best approach to further improve the model's predictive capabilities. It is reassuring to know (based on the present study) that the DF-ANN model for drag correction trained with a comprehensive set of data can bridge both dilute and dense flow conditions.

$$ \frac{\bar{\phi}_s v_{d,z}}{\phi_{s,max} U_t} = f\left(\frac{\bar{\phi}_s}{\phi_{s,max}}, \frac{\tilde{u}_{slip,z}}{U_t}, \frac{\bar{\Delta}}{d_p}, \frac{1}{\rho_s g} \frac{\partial \bar{\rho}_g}{\partial z}, Re_p, Fr_p\right) $$

## NEURAL NETWORK MODELING OF FILTERED SOLID STRESSES

Machine learning-based modeling of the filtered solid phase stresses has also been addressed in the literature, although more sparsely. By relying on an eddy-viscosity concept, Ouyang et al.³⁶ used fine-grid 2-D simulation results to train

a 7-layer ANN and predict the filtered solid phase pressure and the effective viscosity of the solid phase. These authors claim that an anisotropic model using the filtered solid volume fraction, the filter size, the filtered solid phase velocity, and its gradient improves the accuracy of the predictions with respect to isotropic markers, e.g., using only the norm of the rate-of-deformation tensor. Recently, Ouyang et al.²⁹ implemented their ANN solid stress model into an fTFM solver and compared their results with explicit closures. They conclude that anisotropic models are needed for low fluidization velocities in a laboratory-scale dense gas-solid fluidized bed. In this study, we adapt the advanced ANN architecture proposed by Ling et al.³⁵ to model the Reynolds stresses in single-phase flow turbulence with embedded Galilean invariance, and formulate a model for the subgrid solid phase stress in fTFM.

**Artificial Neural Network Architecture for Solid Phase Subgrid Stresses.** To the best of the authors' knowledge, only eddy viscosity-type ANN models inspired by the Smagorinsky⁶⁶ model in single-phase flow turbulence have been developed to close the deviatoric part of the filtered solid
<figure><img src="image_16.png" /><figcaption>Figure 7. Assessment of the DF-ANN model predictions on (a) the drift flux and (b) the filtered drag force using explicit eq 15 with the scalar variance-based 4-marker model ($\bar{\phi}_s$, $\bar{\Delta}$, $\tilde{u}_{slip,z}$, $\bar{\phi}_s'^2$).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 8. Relative error on the vertical component of the filtered drag force predictions for the various ANN models investigated.</figcaption></figure>

stresses. The filtered solid stress tensor is split into its isotropic and deviatoric parts as

$$ \sigma_s = \tau_s + \frac{1}{3} \text{tr}(\sigma_s)\mathbf{I} \qquad (34) $$

The filtered solid pressure, or so-called mesoscale pressure, is given by

$$ P_{s,meso} = \frac{1}{3} \text{tr}(\sigma_s) \qquad (35) $$

Smagorinsky-type models assume the alignment between the deviatoric solid stress tensor and the deviatoric part of the filtered solid rate-of-deformation tensor, namely,

$$ \tau_s = -2\mu_{s,meso}\tilde{\mathbf{S}}_s \qquad (36) $$

where $\mu_{s,meso}$ is the so-called mesoscale viscosity, which is usually estimated from filtered fine-grid data as

$$ \mu_{s,meso} = \frac{\sqrt{\tau_s: \tau_s}}{2\sqrt{\tilde{\mathbf{S}}_s: \tilde{\mathbf{S}}_s}} \qquad (37) $$

So far, ML-based models for the filtered solid stresses sought to describe the mesoscale pressure and mesoscale viscosity through functional models or with distinct ANNs. In single-phase turbulence, the effective viscosity models have known drawbacks: they are completely dissipative, which means that they do not unveil energy backscattering,⁶⁰ and they do not capture accurately anisotropic stresses, even in simple shear flows.⁶¹ In the Reynolds-averaged Navier-Stokes approach for single-phase flows, various ML methods have been recently developed to model the individual components of the turbulent stresses.⁶¹⁻⁶³ One of the most promising ideas has been proposed by Ling et al.³⁵ with a special neural network architecture referred to as a Tensor Basis Neural Network (TBNN). This architecture, sketched in Figure 11, satisfies the
<figure><img src="image_18.png" /><figcaption>**Figure 9.** Assessment of the generalized DF-ANN model prediction with characteristic length scale $L_c = L_{c,I}$: testing on (a) Cases 1–9 (used for the training phase) and (b) Case 10.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 10.** Assessment of the generalized DF-ANN model prediction with characteristic length scale $L_c = d_p$ and $Fr_p$ added as a distinct marker: testing on (a) Cases 1–9 (used for the training phase) and (b) Case 10.</figcaption></figure>

Galilean invariance by relying on the decomposition of the deviatoric stress tensor as a function of the basis of tensors. In this way, the output of the neural network is not modified by an arbitrary rotation of the reference frame, which is a key principle in turbulence model development.

The TBNN approach developed by Ling et al.³⁵ was inspired by the work of Pope,⁶⁴ who showed that, in the single-phase incompressible case, a general eddy viscosity model that is a function of the rate-of-deformation and rate-of-rotation tensors only could be expressed as a linear combination of 10 basis tensors:

$$ \mathbf{a} = \sum_{n=1}^{10} g^{(n)}(\lambda_1, \dots, \lambda_5) \mathbf{T}^{(n)} \quad (38) $$

where $\mathbf{a}$ is the deviatoric stress tensor normalized by the turbulent kinetic energy, while $\lambda_1, \dots, \lambda_5$ (the scalar basis) and $\mathbf{T}^{(1)}...\mathbf{T}^{(10)}$ (the tensor basis) are functions of the filtered strain- and rotation-rate tensors. As mentioned by Ling et al.,³⁵ any

tensor that can be expressed as in eq 38 will satisfy Galilean invariance by construction. The simple eddy-viscosity model is recovered by limiting the summation to $n = 1$. We adopt the same approach for the modeling of the filtered solid phase stresses in gas-solid flows, where, by analogy with Pope,⁶⁴ the tensor and scalar bases can be expressed as

$$ \begin{aligned} \mathbf{T}^{(1)} &= \tilde{\mathbf{S}}_s & \mathbf{T}^{(6)} &= \tilde{\mathbf{R}}_s^2 \tilde{\mathbf{S}}_s + \tilde{\mathbf{S}}_s \tilde{\mathbf{R}}_s^2 - \frac{2}{3} \mathrm{tr}(\tilde{\mathbf{S}}_s \tilde{\mathbf{R}}_s^2) \mathbf{I} \\ \mathbf{T}^{(2)} &= \tilde{\mathbf{S}}_s \tilde{\mathbf{R}}_s - \tilde{\mathbf{R}}_s \tilde{\mathbf{S}}_s & \mathbf{T}^{(7)} &= \tilde{\mathbf{R}}_s \tilde{\mathbf{S}}_s \tilde{\mathbf{R}}_s^2 - \tilde{\mathbf{R}}_s^2 \tilde{\mathbf{S}}_s \tilde{\mathbf{R}}_s \\ \mathbf{T}^{(3)} &= \tilde{\mathbf{S}}_s^2 - \frac{1}{3} \mathrm{tr}(\tilde{\mathbf{S}}_s^2) \mathbf{I} & \mathbf{T}^{(8)} &= \tilde{\mathbf{S}}_s \tilde{\mathbf{R}}_s \tilde{\mathbf{S}}_s^2 - \tilde{\mathbf{S}}_s^2 \tilde{\mathbf{R}}_s \tilde{\mathbf{S}}_s \\ \mathbf{T}^{(4)} &= \tilde{\mathbf{R}}_s^2 - \frac{1}{3} \mathrm{tr}(\tilde{\mathbf{R}}_s^2) \mathbf{I} & \mathbf{T}^{(9)} &= \tilde{\mathbf{R}}_s^2 \tilde{\mathbf{S}}_s^2 - \tilde{\mathbf{S}}_s^2 \tilde{\mathbf{R}}_s^2 - \frac{2}{3} \mathrm{tr}(\tilde{\mathbf{S}}_s^2 \tilde{\mathbf{R}}_s^2) \mathbf{I} \\ \mathbf{T}^{(5)} &= \tilde{\mathbf{R}}_s \tilde{\mathbf{S}}_s^2 - \tilde{\mathbf{S}}_s^2 \tilde{\mathbf{R}}_s & \mathbf{T}^{(10)} &= \tilde{\mathbf{R}}_s \tilde{\mathbf{S}}_s^2 \tilde{\mathbf{R}}_s^2 - \tilde{\mathbf{R}}_s^2 \tilde{\mathbf{S}}_s^2 \tilde{\mathbf{R}}_s \end{aligned} \quad (39) $$
<figure><img src="image_20.png" /><figcaption>Figure 11. Tensor Basis Neural Network (TBNN) architecture for modeling of the filtered solid phase stresses. $\lambda_1, \dots, \lambda_5$ refer to the scalar bases, given by eq 40, and $T^{(1)}...T^{(10)}$ refer to the tensor bases, which are functions of the strain- and rotation-rate tensors given by eq 39.</figcaption></figure>

and

$$ 
\begin{aligned}
\lambda_1 &= \text{tr}(\tilde{\mathbf{S}}_s^2) \\
\lambda_2 &= \text{tr}(\tilde{\mathbf{R}}_s^2) \\
\lambda_3 &= \text{tr}(\tilde{\mathbf{S}}_s^3) \\
\lambda_4 &= \text{tr}(\tilde{\mathbf{R}}_s^2 \tilde{\mathbf{S}}_s) \\
\lambda_5 &= \text{tr}(\tilde{\mathbf{R}}_s^2 \tilde{\mathbf{S}}_s^2)
\end{aligned}
\qquad (40) $$

where $\tilde{\mathbf{S}}_s$ and $\tilde{\mathbf{R}}_s$ are the suitably scaled (see below) rate-of-deformation and rate-of-rotation tensors, respectively, defined by eqs 22 and 23.

Every tensor in the basis defined by eq 39 is traceless and symmetric, consistent with the tensor to be modeled. The goal of the TBNN is therefore to capture the scalar functions $g^{(n)}$ in eq 38 and the construction of the final tensor is performed by the output merge layer as shown in Figure 11. Similar to Ling et al.³⁵ and Pope,⁶⁴ the tensor $\mathbf{a}$ is identified as the deviatoric part of the filtered solid stress scaled by the mesoscale pressure (or, equivalently, the subgrid turbulent kinetic energy):

$$ \mathbf{a} = \frac{\tau_s}{3P_{s,meso}} \triangleq \tau_s^* \qquad (41) $$

The input tensors $\tilde{\mathbf{S}}_s$ and $\tilde{\mathbf{R}}_s$ are scaled using the time-scale $\frac{U_1}{g}$.

The scaled deviatoric part of the solid stress tensor and its eigenvalues ($\xi_1 \ge \xi_2 \ge \xi_3$) must satisfy the following realizability conditions:⁶⁵

$$ 
\begin{aligned}
-\frac{1}{3} \le a_{\mu\mu} \le \frac{2}{3}, & \qquad \xi_1 \ge (3|\xi_2| - \xi_2)/2 \\
-\frac{1}{2} \le a_{\mu\nu} \le \frac{1}{2} & \quad \text{for } \mu \ne \nu, \quad \xi_1 \le \frac{1}{3} - \xi_2
\end{aligned}
\qquad (42) $$

where $\mu, \nu = \{1, 2, 3\}$. Ling et al.³⁵ suggest adding a postprocessing step after the TBNN model to iteratively enforce conditions given by eq 42. From our experience, the number of data points that do not satisfy the realizability condition is fairly limited. Moreover, this additional procedure is computationally costly and does not increase the level of correlation between the predicted values and the data. Therefore, we chose not to keep it for the results presented below.

However, in addition to the TBNN model for the deviatoric part of the stress, a distinct ANN model is needed to predict the mesoscale pressure. Once $P_{s,meso}$ is known, the anisotropic stress $\tau_s$ can be retrieved from eq 41 and the full stress tensor $\sigma_s$ (given by eq 34) can be closed. In addition to the scalar basis defined by eq 40, other scalar inputs specific to the modeling of filtered solid stresses might enter the network, namely, the filtered solid volume fraction and/or the filtered slip velocity. The filter width $\bar{\Delta}$ is also added as an extra marker to account for the variation in the mesh size. In the present study, we only consider the reference case (Case 1) to examine

<figure><img src="image_21.png" /><figcaption>Figure 12. Prediction of the scaled mesoscale solid phase viscosity in Case 1 by (a) the 3-marker and (b) the 14-marker ANN models.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>Figure 13. Prediction of the scaled mesoscale solid phase pressure in Case 1 by (a) the 3-marker and (b) the 14-marker ANN models.</figcaption></figure>

the merits of TBNN-based stress modeling for fTFM analysis. Additional markers accounting for variation in physical properties (such as Re<sub>p</sub>, Fr<sub>p</sub>) will enter the network (as was done for the drift flux model development) when the model is generalized to cover a wide range of gas-particle systems.

## A PRIORI BENCHMARK RESULTS ON ANN SOLID SUBGRID STRESS MODELING

**MLP Models for Mesoscale Solid Pressure and Effective Viscosity.** As a preliminary step, and for a point of comparison, simple MLP models similar to the one used for the filtered drag force are built to predict the mesoscale viscosity and mesoscale pressure. In both cases, the three-layer network architecture detailed earlier is employed with the same number of nodes (124, 32, 8). From our experience, deeper networks did not yield superior results. The mesoscale pressure and viscosity are scaled as follows:

$$ \mu_{s,meso}^* = \frac{\mu_{s,meso} g}{\rho_s U_t^3}, \quad P_{s,meso}^* = \frac{P_{s,meso}}{\rho_s U_t^2} \qquad (43) $$

In line with the work of Ouyang et al.,³⁶ we compare two ANN models for $\mu_{s,meso}^*$ and $P_{s,meso}^*$ with different input markers:

* A 3-marker model with $\lambda_1 = \text{tr}(\tilde{S}_s^2)$, $\bar{\phi}_s$, and $\bar{\Delta}$ as inputs,

* A 14-marker model with $\bar{\phi}_s$, $\bar{\Delta}$ and the components of $\tilde{u}_s$ and $\nabla \tilde{u}_s$ as inputs.

Ouyang et al.³⁶ concluded from their 2-D analysis that the second “anisotropic” model showed improved predictive capacity with respect to the first isotropic version. It must be stressed that this model is anisotropic in the sense that the input markers contain directional information, but the final eddy-viscosity model still relies on the alignment between $\tau_s$ and $\tilde{S}_s$ (see eq 36). Nevertheless, Figures 12 and 13 confirm that the more complete 14-marker ANN model fed with the individual components of $\tilde{u}_s$ and $\nabla \tilde{u}_s$ shows reduced scatter with respect to the simpler 3-marker one, both for the mesoscale pressure and viscosity.

**TBNN Model for the Subgrid Stresses.** It was shown in the previous section that the mesoscale pressure and viscosity can be quite successfully captured with simple ANN models, which can be enough if one is only interested in low-order modeling of the stresses. We now address the TBNN model described earlier, when a more complete description of the

stresses accounting for anisotropic effects is sought. The TBNN network sketched in Figure 11 was built with 8 hidden layers of 30 nodes. The number of layers was initially determined by following the work of Ling et al.³⁵ It is worth noting that we tested the TBNN model with 4 layers (30 nodes each) and found out that the RMSE of the model predictions increases by less than 1%. Although the simpler ANN structure of the TBNN model gave very similar results, in this manuscript, we report the results obtained with eight layers in order to be consistent with Ling et al.³⁵

The activation function of the hidden layers is the ReLU and the loss function is the minimum absolute error (MAE) as for the classical MLP studied above. The MAE is computed on all components of the final deviatoric stress tensor. To assess the predictive capacity of different variants of the TBNN, we also introduce the root mean squared error, computed on the six independent components of $\tau_s^*$:

$$ \text{RMSE} = \sqrt{\frac{1}{6N_{\text{data}}} \sum_{m=1}^{N_{\text{data}}} \sum_{i=1}^{3} \sum_{j=i}^{3} (\hat{\tau}_{s,ij}^* - \tau_{s,ij}^*)^2} \qquad (44) $$

Although the expression given by eq 38 is very general, simpler variants of the decomposition into basis tensors might yield results of similar accuracy. To that end, Table 2 compares

<table><caption>**Table 2.** Comparison between Linear Eddy-Viscosity and Various TBNN Models To Predict the Solid Phase Deviatoric Subgrid Stress Tensor</caption><thead><tr><th>model</th><th>markers</th><th>RMSE</th></tr></thead><tbody><tr><td>eddy-viscosity ANN model</td><td>$\lambda_1, \bar{\phi}_s, \bar{\Delta}$</td><td>0.4170</td></tr><tr><td>eddy-viscosity ANN model</td><td>$\tilde{u}_s, \nabla \tilde{u}_s, \bar{\phi}_s, \bar{\Delta}$</td><td>0.4048</td></tr><tr><td>TBNN model: $T^{(1)}, T^{(3)}$</td><td>$\lambda_1, \lambda_3, \bar{\phi}_s, \bar{\Delta}$</td><td>0.2076</td></tr><tr><td>TBNN model: $T^{(1)}-T^{(4)}$</td><td>$\lambda_1-\lambda_3, \bar{\phi}_s, \bar{\Delta}$</td><td>0.0697</td></tr><tr><td>TBNN model: $T^{(1)}-T^{(4)}$</td><td>$\lambda_1-\lambda_3, \bar{\phi}_s, \tilde{u}_{slip}, \bar{\Delta}$</td><td>0.0695</td></tr><tr><td>TBNN model: $T^{(1)}-T^{(4)}$</td><td>$\lambda_1-\lambda_5, \bar{\phi}_s, \bar{\Delta}$</td><td>0.0697</td></tr><tr><td>TBNN model: $T^{(1)}-T^{(10)}$</td><td>$\lambda_1-\lambda_5, \bar{\phi}_s, \bar{\Delta}$</td><td>0.0692</td></tr></tbody></table>

the RMSE of different versions of the TBNN model, where we vary both the number of tensors in the basis and the scalars fed to the network. We also introduce in Table 2 the RMSE of the eddy-viscosity models studied in the previous section for comparison. First, it comes out that, even though the 14-marker eddy-viscosity model was shown to be more accurate
<figure><img src="image_23.png" /><figcaption>**Figure 14.** Prediction of the TBNN model built with four basis tensors ($T^{(1)}-T^{(4)}$) and 5 input markers ($\lambda_1, \lambda_2, \lambda_3, \bar{\phi}_s, \bar{\Delta}$) in Case 1 for (a) the 6 independent components of $\tau_s^*$ and (b) the square of the norm of $\tau_s^*$.</figcaption></figure>

than the 3-marker model to predict the mesoscale viscosity, i.e. the norm of the stress tensor, both models have similar RMSE computed on the individual components of the deviatoric stress tensor. Thus, the improvement achieved by increasing the number of markers from four to fourteen is only marginal. In contrast, all variants of the TBNN model show better accuracy. This highlights the limitation of the eddy-viscosity concept for subgrid stress modeling in gas-solid flows (just as in turbulent single-phase flows). A first TBNN model can be built only with tensors $T^{(1)}$ and $T^{(3)}$, and the scalars $\lambda_1$ and $\lambda_3$, i.e. discarding quantities that involve the rotation-rate tensor $\tilde{R}_s$. This model shows a reduced RMSE with respect to the two eddy-viscosity models. However, models that account for the first four tensors in the basis display a much better performance. Among these models, we observe that using only the first three scalars $\lambda_1-\lambda_3$ defined by eq 40 appear sufficient. The last two scalars $\lambda_4$ and $\lambda_5$ do not seem to provide extra information in order to capture the solid phase subgrid stress tensor. Likewise, the addition of the slip velocity vector as an input marker only marginally reduces the RMSE value. Finally, using the complete tensor basis ($T^{(1)}-T^{(10)}$) reduces only slightly the error made on the predicted stress tensor, while increasing the complexity of the network and the computational cost of the convergence algorithm. Therefore, for practical use in fTFM simulations, we suggest using the TBNN model with the first four basis tensors ($T^{(1)}-T^{(4)}$), and, as an input to the network, the first three scalars of Pope's basis ($\lambda_1-\lambda_3$), the filter size $\bar{\Delta}$ and the filtered solid volume fraction $\bar{\phi}_s$.

The relative contribution of the different basis tensors in the construction of the final stress tensor can be estimated by the mean absolute value of the scalar functions $g^{(1)}$ to $g^{(4)}$. This analysis shows that the eddy-viscosity term, i.e., $g^{(1)}$, only contributes to 0.6% of the final model, while the second, third, and fourth tensors contribute respectively to 57.9, 24.6, and 16.9%. Future works should therefore investigate why the eddy-viscosity term virtually vanishes when a more complete tensor basis is used to build the deviatoric stress tensor and whether this conclusion is valid for a wide range of regimes.

Figure 14 shows the parity plots between the predicted values and the filtered fine-grid data for the 6 individual

components of tensor $\tau_s^*$ and for the square of its norm $\tau_s^*:\tau_s^*$. It appears that the scatter increases substantially in the second case. If the norm of the stress tensor is of higher importance for the prospective fTFM simulation, the loss function of the network could be customized for a more balanced error between the individual components and the norm of the tensor.

## ■ PUBLICLY AVAILABLE ANN PYTHON LIBRARY FOR FTFM CLOSURES

As we discussed before, the Python ML model source codes are available in the GitHub repository: https://github.com/bahardy/fTFM_ANN_modeling.git. Instructions for interested end-users who need to develop ANN models for the filtered drag force and the filtered solid phase stresses through TBNN using their own filtered data sets are as follows: The code reads the filtered data set from a `.txt` file, which should be generated by the end-user by filtering fine-grid simulation results. As an illustrative example, and to avoid sharing a very large data set of filtered simulation results generated by `neptune_cfd`,⁵⁰ we have uploaded about 6% of the filtered simulation results for a benchmark case. After training the networks, the output files from the training process are saved using Keras API:⁶⁶ (1) a JSON file that contains the neural network structure; (2) an HDF5 file that contains weight information required for the evaluation of the neural network model. These files could then be read by an open-source interface `Fortran/C++` code and are suitable for `MFIX` and `OpenFOAM` simulation platforms. This implementation allows the user to read the flow quantities during simulation runtime and evaluate the prediction with the ANN models. With ANN models and implementation approach, it is possible to assess its accuracy in a posteriori simulations.

## ■ CONCLUSIONS

Jiang et al.²²,²³ proposed an artificial neural network (ANN) model for the dense flow regime in which the filter size, the filtered particle volume fraction, the filtered slip velocity, and the filtered gas phase pressure gradient, which are available in an fTFM model simulation, are used to estimate the drift flux.
The drift flux is then used to estimate the correction of the drag force for fTFM simulations, which are feasible for industrial-scale gas-solid flows. Rauchenzauner and Schneiderbauer²⁰ used the filtered particle volume fraction, the scalar variance of the subgrid particle volume fraction variation, and the kinetic energy associated with the subgrid velocity fluctuations of the particles to find the drift velocity in the dense flow regime. Both approaches work well when applied to a common data set generated by filtering the results from a dense fluidized bed simulation.

We then extended these studies through gas-particle fluidization simulations in periodic domains in the dilute regime and examined several different approaches to finding the correction for the drag force needed for fTFM models. It was found that the approach adopted by Jiang et al.²³ for the dense flow regime works well for the dilute flow regime as well. This implies that a single ANN model that covers both regimes can be found by pooling together the dense and dilute flow regime data sets.

Furthermore, we introduce a Galilean-invariant tensor-based neural network (TBNN) model to capture the anisotropic particle phase stress stemming from the subgrid velocity fluctuations, which need closure for the fTFM approach. The proposed approach first utilizes distinct ANNs to find the filtered solid phase pressure and effective viscosity, which is a classical way of turbulence modeling in single-phase flows for subgrid velocity fluctuations. It then employs a TBNN model to find the components of the filtered solid phase stress tensor. It is demonstrated that the TBNN approach captures the anisotropy quite nicely.

Future work should strive to generate a comprehensive drift flux model that combines the data sets generated through dense and dilute flow simulations. It should also examine how the level of sophistication of the stress model, a simple Smagorinsky-like eddy viscosity model vs the TBNN model allowing for anisotropy, influences the predictions of fTFM simulations. Besides, the present data-based filtered drag and subgrid solid stress models should be compared with earlier functional anisotropic closures proposed in the literature,³²,³³ and their respective performance assessed. A further step will be to perform fTFM simulations, namely, a posteriori tests, with the developed models, and compare the predictions with the fine-grid TFM simulations and the experimental studies to assess their accuracy.

## AUTHOR INFORMATION

### Corresponding Author

Baptiste Hardy - Process and Energy Department, Faculty of Mechanical Engineering, TU Delft, Delft 2628 CB, Netherlands; ● orcid.org/0000-0001-8076-9546; Email: b.hardy@tudelft.nl

### Authors

Stefanie Rauchenzauner - Process Systems Engineering, TUM School of Life Sciences, Technical University of Munich, Freising 85354, Germany

Pascal Fede - Institut de Mécanique des Fluides de Toulouse, Université de Toulouse, CNRS, Toulouse 31400, France

Simon Schneiderbauer - Department of Particulate Flow Modelling, Johannes Kepler University, Linz 4040, Austria; ● orcid.org/0000-0001-7690-9884

Olivier Simonin - Institut de Mécanique des Fluides de Toulouse, Université de Toulouse, CNRS, Toulouse 31400, France

Sankaran Sundaresan - Department of Chemical and Biological Engineering, Princeton University, Princeton, New Jersey 08544, United States; ● orcid.org/0000-0002-5868-1975

Ali Ozel - School of Engineering and Physical Sciences, Heriot-Watt University, Edinburgh EH14 4AS, U.K.; ● orcid.org/0000-0002-0133-2002

Complete contact information is available at: https://pubs.acs.org/10.1021/acs.iecr.3c04652

### Notes

The authors declare no competing financial interest.

## ACKNOWLEDGMENTS

The netpune_cfd software is a multiphase CFD solver developed in the framework of the NEPTUNE project financially supported by EDF (Electricité de France), CEA (Commissariat à l'Energie Atomique), IRSN (Institut de Radioprotection et de Sureté Nucléaire) and Framatome. This work was financially supported by the Carnot ISIFOR Institute under the REBIS project (2022/2023) and by an Excellence Research Grant from Wallonie-Bruxelles International 2023/2024 (reference SOR/2023/550158). It also benefited from HPC resources provided by the supercomputing centers CALMIP (project P0111) and GENCI (allocation A0142B06012). A preprint version of this work prior to any peer-reviewing is publicly available on the arXiv preprint server at 10.48550/arXiv.2401.00179.

## REFERENCES

(1) Anderson, T. B.; Jackson, R. Fluid Mechanical Description of Fluidized Beds. Equations of Motion. Ind. Eng. Chem. Fundam. 1967, 6, 527-539.

(2) Gidaspow, D. Multiphase Flow and Fluidization; Academic Press: San Diego, 1994.

(3) Glasser, B. J.; Kevrekidis, I.; Sundaresan, S. One-and two-dimensional travelling wave solutions in gas-fluidized beds. J. Fluid Mech. 1996, 306, 183-221.

(4) Glasser, B. J.; Kevrekidis, I.; Sundaresan, S. Fully developed travelling wave solutions and bubble formation in fluidized beds. J. Fluid Mech. 1997, 334, 157-188.

(5) Glasser, B. J.; Sundaresan, S.; Kevrekidis, I. From bubbles to clusters in fluidized beds. Phys. Rev. Lett. 1998, 81, 1849.

(6) Agrawal, K.; Loezos, P. N.; Syamlal, M.; Sundaresan, S. The role of meso-scale structures in rapid gas-solid flows. J. Fluid Mech. 2001, 445, 151-185.

(7) Parmentier, J.-F.; Simonin, O.; Delsart, O. A functional subgrid drift velocity model for filtered drag prediction in dense fluidized bed. AIChE J. 2012, 58, 1084-1098.

(8) Ozel, A.; Fede, P.; Simonin, O. Development of filtered Euler-Euler two-phase model for circulating fluidised bed: High resolution simulation, formulation and a priori analyses. Int. J. Multiphase Flow 2013, 55, 43-63.

(9) Sundaresan, S. Modeling the hydrodynamics of multiphase flow reactors: current status and challenges. AIChE J. 2000, 46, 1102-1105.

(10) Igci, Y.; Andrews, A. T.; Sundaresan, S.; Pannala, S.; O'Brien, T. J. Filtered Two-Fluid Models for Fluidized Gas-Particle Suspensions. AIChE J. 2008, 54, 1431-1448.

(11) Schneiderbauer, S. A spatially-averaged two-fluid model for dense large-scale gas-solid flows. AIChE J. 2017, 63, 3544-3562.
(12) Sundaresan, S.; Ozel, A.; Kolehmainen, J. Toward Constitutive Models for Momentum, Species, and Energy Transport in Gas-Particle Flows. *Annu. Rev. Chem. Biomol. Eng.* **2018**, *9*, 61–81.

(13) Pope, S. B. Turbulent flows. *Meas. Sci. Technol.* **2001**, *12*, 2020–2021.

(14) Li, J.; Kwauk, M. *Particle-Fluid Two-Phase Flow: The Energy-Minimization Multi-Scale Method*; Metallurgical Industry Press: Beijing, China, 1994.

(15) Lu, B.; Wang, W.; Li, J. Searching for a mesh-independent sub-grid model for CFD simulation of gas-solid riser flows. *Chem. Eng. Sci.* **2009**, *64*, 3437–3447.

(16) Pakseresht, P.; Yao, Y.; Fan, Y.; Theuerkauf, J.; Capecelatro, J. A critical assessment of the Energy Minimization Multi-Scale (EMMS) model. *Powder Technol.* **2023**, *425*, No. 118569.

(17) Milioli, C. C.; Milioli, F. E.; Holloway, W.; Agrawal, K.; Sundaresan, S. Filtered two-fluid models of fluidized gas-particle flows: new constitutive relations. *AIChE J.* **2013**, *59*, 3265–3275.

(18) Sarkar, A.; Milioli, F. E.; Ozarkar, S. S.; Li, T.; Sun, X.; Sundaresan, S. Filtered sub-grid constitutive models for fluidized gas-particle flows constructed from 3-D simulations. *Chem. Eng. Sci.* **2016**, *152*, 443–456.

(19) Ozel, A.; Kolehmainen, J.; Radl, S.; Sundaresan, S. Fluid and particle coarsening of drag force for discrete-parcel approach. *Chem. Eng. Sci.* **2016**, *155*, 258–267.

(20) Rauchenzauner, S.; Schneiderbauer, S. Validation study of a Spatially-Averaged Two-Fluid Model for heat transport in gas-particle flows. *Int. J. Heat Mass Transfer* **2022**, *198*, No. 123382.

(21) Hardy, B.; Fede, P.; Simonin, O. Theoretical derivation and a priori validation of a new scalar variance-based sub-grid drag force model for simulation of gas-solid fluidized beds. *Powder Technol.* **2024**, *436*, No. 119454.

(22) Jiang, Y.; Kolehmainen, J.; Gu, Y.; Kevrekidis, Y. G.; Ozel, A.; Sundaresan, S. Neural-network-based filtered drag model for gas-particle flows. *Powder Technol.* **2019**, *346*, 403–413.

(23) Jiang, Y.; Chen, X.; Kolehmainen, J.; Kevrekidis, I. G.; Ozel, A.; Sundaresan, S. Development of data-driven filtered drag model for industrial-scale fluidized beds. *Chem. Eng. Sci.* **2021**, *230*, No. 116235.

(24) Zhang, Y.; Jiang, M.; Chen, X.; Yu, Y.; Zhou, Q. Modeling of the filtered drag force in gas-solid flows via a deep learning approach. *Chem. Eng. Sci.* **2020**, *225*, No. 115835.

(25) Zhu, L.-T.; Ouyang, B.; Lei, H.; Luo, Z.-H. Conventional and data-driven modeling of filtered drag, heat transfer, and reaction rate in gas-particle flows. *AIChE J.* **2021**, *67*, No. e17299.

(26) Tausendschön, J.; Sundaresan, S.; Salehi, M.; Radl, S. Machine Learning-Based Filtered Drag Model for Cohesive Gas-Particle Flows. *Chem. Eng. Technol.* **2023**, *46*, 1373–1386.

(27) Lu, L.; Gao, X.; Dietiker, J.-F.; Shahnam, M.; Rogers, W. A. Development of a Filtered CFD-DEM Drag Model with Multiscale Markers Using an Artificial Neural Network and Nonlinear Regression. *Ind. Eng. Chem. Res.* **2022**, *61*, 882–893.

(28) Zhu, L.-T.; Chen, X.-Z.; Ouyang, B.; Yan, W.-C.; Lei, H.; Chen, Z.; Luo, Z.-H. Review of Machine Learning for Hydrodynamics, Transport, and Reactions in Multiphase Flows and Reactors. *Ind. Eng. Chem. Res.* **2022**, *61*, 9901–9949.

(29) Ouyang, B.; Zhu, L.-T.; Su, Y.-H.; Luo, Z.-H. A hybrid mesoscale closure combining CFD and deep learning for coarse-grid prediction of gas-particle flow dynamics. *Chem. Eng. Sci.* **2022**, *248*, No. 117268.

(30) Fox, R. O. On multiphase turbulence models for collisional fluid-particle flows. *J. Fluid Mech.* **2014**, *742*, 368–424.

(31) Schneiderbauer, S. Validation study on spatially averaged two-fluid model for gas-solid flows: I. A priori analysis of wall bounded flows. *AIChE J.* **2018**, *64*, 1591–1605.

(32) Cloete, J. H.; Cloete, S.; Municchi, F.; Radl, S.; Amini, S. Development and verification of anisotropic drag closures for filtered Two Fluid Models. *Chem. Eng. Sci.* **2018**, *192*, 930–954.

(33) Cloete, J. H.; Cloete, S.; Radl, S.; Amini, S. Development and verification of anisotropic solids stress closures for filtered Two Fluid Models. *Chem. Eng. Sci.* **2018**, *192*, 906–929.

(34) Rauchenzauner, S.; Schneiderbauer, S. A Dynamic Spatially-Averaged Two-Fluid Model for Heat Transport in Moderately Dense Gas-Particle Flows. *Phys. Fluids* **2020**, *32*, No. 063307.

(35) Ling, J.; Kurzawski, A.; Templeton, J. Reynolds averaged turbulence modelling using deep neural networks with embedded invariance. *J. Fluid Mech.* **2016**, *807*, 155–166.

(36) Ouyang, B.; Zhu, L.-T.; Luo, Z.-H. Data-driven modeling of mesoscale solids stress closures for filtered two-fluid model in gas-particle flows. *AIChE J.* **2021**, *67*, No. e17290.

(37) Zhang, D. Z.; Vander Heyden, W. B. The effects of mesoscale structures on the macroscopic momentum equations for two-phase flows. *Int. J. Multiphase Flow* **2002**, *28*, 805–822.

(38) Garnier, E.; Adams, N.; Sagaut, P. *Large eddy simulation for compressible flows*; Springer Science & Business Media, 2009.

(39) Wen, C.; Yu, Y. Mechanics of fluidization. In *Chemical Engineering Symposium Series*; 1966, pp 100-111.

(40) Wang, J.; van der Hoef, M. A.; Kuipers, J. A. M. Coarse grid simulation of bed expansion characteristics of industrial-scale gas-solid bubbling fluidized beds. *Chem. Eng. Sci.* **2010**, *65*, 2125–2131.

(41) Schneiderbauer, S.; Pirker, S. Filtered and heterogeneity based sub-grid modifications for gas-solid drag and solids stresses in bubbling fluidized beds. *AIChE J.* **2014**, *60*, 839–854.

(42) Schneiderbauer, S.; Puttinger, S.; Pirker, S.; Aguayo, P.; Kanellopoulos, V. CFD modeling and simulation of industrial scale olefin polymerization fluidized bed reactors. *Chem. Eng. J.* **2015**, *264*, 99–112.

(43) Igci, Y.; Sundaresan, S. Constitutive Models for Filtered Two-Fluid Models of Fluidized Gas-Particle Flows. *Ind. Eng. Chem. Res.* **2011**, *50*, 13190–13201.

(44) Schneiderbauer, S.; Saeedipour, M. Approximate deconvolution model for the simulation of turbulent gas-solid flows: An a-priori analysis. *Phys. Fluids* **2018**, *30*, No. 023301.

(45) Ozel, A.; Gu, Y.; Milioli, C. C.; Kolehmainen, J.; Sundaresan, S. Towards filtered drag force model for non-cohesive and cohesive particle-gas flows. *Phys. Fluids* **2017**, *29*, 103308.

(46) Rauchenzauner, S.; Schneiderbauer, S. A Dynamic Anisotropic Spatially-Averaged Two-Fluid Model for Moderately Dense Gas-Particle Flows. *Int. J. Multiphase Flow* **2020**, *126*, No. 103237.

(47) Dabbagh, F.; Schneiderbauer, S. Small-scale flow topologies, pseudo-turbulence, and impact on filtered drag models in turbulent fluidization. *Phys. Rev. Fluids* **2021**, *6*, No. 074301.

(48) Andrews, A. T.; Loezos, P. N.; Sundaresan, S. Coarse-Grid Simulation of Gas-Particle Flows in Vertical Risers. *Ind. Eng. Chem. Res.* **2005**, *44*, 6022–6037.

(49) De Wilde, J. Reformulating and quantifying the generalized added mass in filtered gas-solid flow models. *Phys. Fluids* **2005**, *17*, 113304.

(50) Neau, H.; Pigou, M.; Fede, P.; Ansart, R.; Baudry, C.; Mérigoux, N.; Laviéville, J.; Fournier, Y.; Renon, N.; Simonin, O. Massively parallel numerical simulation using up to 36,000 CPU cores of an industrial-scale polydispersed reactive pressurized fluidized bed with a mesh of one billion cells. *Powder Technol.* **2020**, *366*, 906–924.

(51) Schiller, L. A drag coefficient correlation. *Zeit. Ver. Deutsch. Ing. 1933*, *77*, 318–320.

(52) Cloete, S.; Cloete, J. H.; Amini, S. Hydrodynamic validation study of filtered Two Fluid Models. *Chem. Eng. Sci.* **2018**, *182*, 93–107.

(53) Zhu, L.-T.; Tang, J.-X.; Luo, Z.-H. Machine learning to assist filtered two-fluid model development for dense gas-particle flows. *AIChE J.* **2020**, *66*, No. e16973.

(54) Ouyang, B.; Zhu, L.; Luo, Z. Interpretable machine learning analysis and automated modeling to simulate fluid-particle flows. *Particuology* **2023**, *80*, 42–52.

(55) Jiang, M.; Chen, X.; Zhou, Q. A gas pressure gradient-dependent subgrid drift velocity model for drag prediction in fluidized gas-particle flows. *AIChE J.* **2020**, *66*, No. e16884.

(56) Smagorinsky, J. General circulation experiments with the primitive equations: I. basic experiment. *Mon. Weather Rev.* **1963**, *91*, 99–164.
(57) Yoshizawa, A. Statistical theory for compressible turbulent shear flows, with the application to subgrid modeling. *Phys. Fluids* **1986**, *29*, 2152–2164.

(58) Cloete, J. H.; Cloete, S.; Radl, S.; Amini, S. On the choice of closure complexity in anisotropic drag closures for filtered Two Fluid Models. *Chem. Eng. Sci.* **2019**, *207*, 379–396.

(59) Radl, S.; Sundaresan, S. A drag model for filtered Euler-Lagrange simulations of clustered gas-particle suspensions. *Chem. Eng. Sci.* **2014**, *117*, 416–425.

(60) Dabbagh, F.; Schneiderbauer, S. Anisotropy characterization of turbulent fluidization. *Phys. Rev. Fluids* **2022**, *7*, No. 094301.

(61) Tracey, B.; Duraisamy, K.; Alonso, J.; Application of supervised learning to quantify uncertainties in turbulence and combustion modeling. In *51st AIAA Aerospace Sciences Meeting*, 2013; p 259.

(62) Tracey, B. D.; Duraisamy, K.; Alonso, J. J.; A machine learning strategy to assist turbulence model development. In *53rd AIAA Aerospace Sciences Meeting*, 2015; p 1287.

(63) Zhang, Z. J.; Duraisamy, K.; Machine learning methods for data-driven turbulence modeling. In *22nd AIAA Computational Fluid Dynamics Conference*, 2015; p 2460.

(64) Pope, S. B. A more general effective-viscosity hypothesis. *J. Fluid Mech.* **1975**, *72*, 331–340.

(65) Banerjee, S.; Krahl, R.; Durst, F.; Zenger, C. Presentation of anisotropy properties of turbulence, invariants versus eigenvalue approaches. *J. Turbul.* **2007**, *8*, 1–27.

(66) Chollet, F. *Deep learning with Python*; Simon and Schuster, 2021.