Machine learning approaches to close the filtered two-fluid model for gas-solid flows: Models for subgrid drag force and solid phase stress

Gas-particle flows are commonly simulated through two-fluid model at industrial-scale. However, these simulations need 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 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 datasets. Then, these datasets 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 capture nicely 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 with those obtained from a simple eddy-viscosity ANN model.


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 computeraided 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), 1,2 which treat the fluid and particle phases as inter-penetrating continua, are more viable to analyze and simulate the flows in industrialscale 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 [3][4][5] ).The TFMs can readily be solved numerically using commercial codes (e.g., ANSYS Fluent), and opensource 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 3 − 10 5 particle diameters.The macroscale flow structures, referred to as coherent flow structures, can have a large effect on the overall flow hydrodynamics in the device.At the same time, meso-scale 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 diameters [6][7][8][9] 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. 10This consideration led to the development of filtered two-fluid models (fTFMs) 7,8,[11][12][13] 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. 14e fTFM equations contain several terms representing the consequences of subfilterscale 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 7,8,12 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, 8 while all the other corrections are essentially negligible.As a result, the literature on fTFM model development has focused primarily on the 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 meso-scale in the literature can be classified as follows.The Energy Minimization Multi Scale (EMMS) model [15][16][17] 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. 11 , Milioli et al. 18 , Sarkar et al. 19 for the filtered drag in terms of the filtered variables and the filter size.Parmentier et al. 7 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. 7,8,20Rauchenzauner and Schneiderbauer 21 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. 22 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. 23 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 multi-layer perceptron, (MLP)) for the drift flux (which is the product of filtered drift velocity and the filtered solid volume fraction), for 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. 7 .Jiang et al. 24 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. 23 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. 25 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.All these studies considered systems without inter-particle forces; Tausendschön et al. 26 report that the strength of the drag correction is affected by inter-particle forces and proposed an ANN model including the Bond number as a measure of the attractive inter-particle forces.
In the Eulerian-Lagrangian approach, Lu et al. 27 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. 28 .
In spite of the progress in the application of ML methods for formulating subgrid drag models for gas-solid flows without the inter-particle 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. 24 using datasets generated through simulations of dense fluidized beds be sufficient to model the filtered drift flux in dilute flows such as those encountered in risers?
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 capturing correctly 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. 8 , Milioli et al. 18 , Schneiderbauer 29 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.1][32] Rauchenzauner and Schneiderbauer 32 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, 33 to constitute the 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. 34Finally, a sample dataset 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 paper 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 datasets are detailed.Subsequently, the Artifical 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 paper 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. 11 and others, 7,8,35 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: Here, ρ g , ρ s are the gas and solid phase densities, respectively, p g is the gas phase pressure, and g is the gravitational acceleration.The filtered volume fractions for gas (k = g) and solid (k = s) phases are defined as where ϕ k (x, t) is the volume fraction for each phase given by the "micro-scale" TFM, G(r) is the filter convolution kernel satisfying V G(r)dr = 1.Similar to LES of compressible flows, 36 the phase velocities are filtered through the Favre-averaging as The filtered quantity for each phase, denoted with a bar, in Eq. ( 1) is defined as: where Q k (y, t) is a quantity for each phase in the "micro-scale" TFM.The filtered and its fluctuating quantities are described as Σd g and Σs are respectively the deviatoric filtered gas phase and total solid phase micro-scale stresses.Ī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. 7,8,11e filtered drag force term is defined as where β is the microscopic drag coefficient (also defined as β = ρ s ϕ s τ p , with τ p the particle relaxation time 7,8 ).In the literature, it is very common practice 12,15,16,26,35,[37][38][39] to account for the drag correction required in the fTFM by introducing an effective drag coefficient β e , namely, where β 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 β, is usually expressed in terms of a drag correction factor H d 16 defined as Numerous studies have sought to improve the functional description of H d . 7,8,11,18,19,40It 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. 26 , a more general model would be where H d is a second-order tensor.Nevertheless, it has been verified 7,8,41 that Eq. ( 11) is adequate and that, to a very good approximation, the filtered drag force can be written as where v d is the so-called drift velocity, defined as In this paper, the product φs v d will be referred to as the "drift flux".The filtered drag and the drift flux can be explicitly linked as follows where τp is the particle relaxation time computed from filtered quantities.Most functional models proposed in the literature considered v d to be aligned with the resolved slip velocity 7,8 (which eventually boils down to an expression of the form given by Eq. ( 10)) though, Ozel et al. 42 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 43 while Hardy et al. 22 recently proposed an explicit model to deduce the drift velocity from the subgrid variance, independently of the filter size.
The second most important sub-grid contribution in the fTFM comes from the filtered solid phase (or meso-scale) stresses, defined as While the drag correction has been found to saturate as the filter size increases, the meso-scale stress term grows monotonically with the filter size. 8Its modeling is therefore becoming increasingly important as the mesh resolution is lowered for very large scale simulations.Prior studies 6,11,44 have concluded that the filtered solid phase micro-scale stresses (Σ s ) described by the kinetic theory of granular flows is much weaker than the meso-scale 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 Some authors 35,45 proposed to model Φ sgs k as an added mass term, but most studies concluded that this term was small with respect to the filtered drag and meso-scale stress terms and that it is sufficient to retain only the resolved part of the pressure gradient, i.e. ϕ k ∇p g ≃ φk ∇p g .

Flow Configurations and Generating Dataset 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. 24 performed fine-grid simulations of dense fluidized beds.Their datasets 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 at more dilute conditions.One of the goals of the present 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 . 46e computational domain is tri-periodic 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 the gravitational acceleration.The different investigated cases and their physical parameters (particle diameter d p , density ρ s and domain-averaged solid fraction ⟨ϕ s ⟩) are summarized in Table 1.The following variables are fixed through all simulations: the gas density ρ g = 1.2 kg m −3 , the gas dynamic viscosity µ g = 1.8 × 10 −5 Pa s and the particle restitution coefficient e c = 0.9.We also report the particle Froude number Fr p = U  1. The mesh size ∆ is uniform, with 640 cells along the vertical direction, and the aspect ratio between the vertical and horizontal dimensions of the domain is 4.
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 has to 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. 11In the present work, we use a box (or top-hat) filter   1).
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 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 details, but a few key quantities can already be identified.
The filtered slip velocity appears in numerous functional models next to the filtered solid volume fraction φs and the filter size ∆ (see Cloete et al. 48for a comparative study of different existing models).It has been more recently observed by Jiang et al. 23 and confirmed by later studies 24,25,49,50 that the addition of the filtered pressure gradient acting against the gravitational acceleration ∂ pg ∂z as an input to their filtered drag ANN model dramatically improved their results.
Based on that, Jiang et al. 51 formulated an explicit drift velocity model with an additional dependence on the filtered pressure gradient.Besides, recent studies 8,22,42,43 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 referred to as SV) is a good measure of this level of heterogeneity.In this study, the SV will be defined as The SV was previously introduced for filtered drag force modeling by Schneiderbauer 12 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 sgs s , defined as As for the filtered solid stresses modeling, single phase turbulence models 52,53 and previous functional modeling efforts for gas-solid flows 8,18,29 suggest that the deviatoric part of the filtered rate-of-deformation tensor, defined as should definitely play a role.The rotation-rate tensor 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 ∇ φs have been extracted from our fine-grid TFM simulation results as they contain non-local 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. 23 .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 φs , the filtered slip velocity ũslip,z and the gas phase gradient ∂ pg ∂z in the gravitational acceleration direction.
Instead of proposing an explicit functional form for the 3-marker model, they used an MLP made of 3 hidden layers of 128, 64 and 32 nodes to close the component of the drift flux aligned with gravitational acceleration, i.e. φ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 output variable of their ANN, instead of the drift flux.In a later study, 24 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 physical properties.The prediction improvement offered by the addition of the filtered gas phase pressure gradient was also confirmed by Ouyang et al. 50through interpretable ML metrics.Zhang et al. 25 proposed a convolutional neural netwok (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 that with the MLP architecture since information from neighboring cells is already embedded in the model.It is however expected that such CNN model will be more computationally intensive and therefore less tractable for practical fTFM simulations at industrial-scale.Zhu et al. 49 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. 26 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 that the isotropic counterpart, as previously noticed by Cloete et al. 30 , 54 for explicit models.These authors also added the Bond number as an additional marker to account for cohesive effects in gas-solid flows.

Artifical 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 sub-optimal neural networks, we developed our NN architecture, which is a feedforward MLP similar to Jiang et al. 23 .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 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 the nodes in the previous layer.
The outcome of a single neuron i within a layer n is described by where z (n−1) is the output vector of layer n − 1, 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 numbers 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 datasets generated for dilute systems in the present study, To that end, we investigated the following strategies: where N is the number of data points used in a training batch.The dataset for the reference case is made of 1.792 × 10 6 entries (accounting for the 7 filter widths).Training and testing phases use subsets corresponding to 80% and 20% of the dataset, respectively.Among the training dataset, 20% of the data are further preserved for validation and 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 . . . . . . . . . . . .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), ŷi is the corresponding ANN model output value and ⟨y⟩ is the mean value of y over the dataset.The quality of the model can also be quantified by the Probability Density Function (PDF) of the relative error e, defined as 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.
where the characteristic length scale L c is set equal to d p Fr 1/3 p , as suggested by Radl and Sundaresan 55 .
Rauchenzauner and Schneiderbauer 43 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 fluctua-tions (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 43 and Jiang et al. 24 employed different sets of input variables to model the drift velocity.To test whether the dataset generated by Rauchenzauner and Schneiderbauer 43 could have been captured by the ANN model supplied by Jiang et al. 24 , we compared the predictions of the Jiang et al. 24 model with the data generated by Rauchenzauner and Schneiderbauer 43 .(Specifically, we tested the results from case 2 of Rauchenzauner and Schneiderbauer 43 which closely corresponded to case 2 of Jiang et al. 24 .) Figure 3 reveals fairly 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. 24 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. 24 ; 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. 24 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 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 and In line with Jiang et al. 23 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. 25 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. 23 found that FD-ANN model performed poorly for dense flows, which is clearly not the case for dilute flows.It can be inferred from Figure 7 that the 4-marker ANN model described by  is capable of predicting the drift flux in the vertical direction with a 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 SV for this reference case.
Zhang et al. 25 argued that the flow information of 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 machine learning approach, we should identify which differential quantities adds 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 dense and dilute regimes.

Generalized ANN model
For a 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 datasets to cover a large range of these dimensionless input quantities.This approach has been investigated by Jiang et al. 24 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 8 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 the filtered drag force for constant physical parameters, a more general DF-ANN model can expressed as 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. 24 in the dense regime to reduce the number of independent variables is given by Eq. ( 28 for testing purpose (left side of Figure 9).It can be observed that the generalized ANN model is able to predict the filtered drag force with decent accuracy (R 2 = 0.88).The model is also tested on Case 9 (right side of Figure 9), which was not used for the training process.
This test case aims to verify that the scaling proposed in Eq. ( 28) is justified in the dilute regime and that the range of parameters spanned by the training dataset (Cases 1 to 8) is sufficient to build a robust generalized model.It is shown that the scatter of the model increases slightly (R 2 = 0.77).This observation may have several causes: • the set of dimensionless markers in Eq. ( 28) is imperfect or incomplete; • the range of parameters covered by the dataset formed by Cases 1 to 8 is not sufficient to train a robust, generalized DF-ANN model and additional data should be fed to the neural network; or • the neural network architecture is sub-optimal.
Though it is hard to make a definitive conclusion, our attempts to increase the complexity of the network (number of layers and number of nodes) did not yield improved predictions.
By changing the definition of the characteristic length scale, using L c = L II , we observed similar performances on Cases 1 to 8, but a dramatic decline of the model performance on Case 9.Then, considering that the characteristic length scale could be more generally written as where 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.
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 the best option for now, also in the dilute or moderately dense regime.
Nevertheless, future studies should clarify whether a more advanced ANN architecture or a different set of dimensionless markers could improve the predictive capacity.We also note that we used 8 different cases to train our model while Jiang et al. 24 trained their model with more than 20 different cases.Thus, enlarging the dilute regime datasets in future studies to include more cases spanning a wider range of parameters 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 DF-ANN model for drag correction trained with a comprehensive set of data can bridge both dilute and dense flow conditions.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 where µ s,meso is the so-called meso-scale viscosity, which is usually estimated from filtered fine-grid data as So far, ML-based models for the filtered solid stresses sought to describe the meso-scale pressure and meso-scale 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, 57 and they do not capture accurately anisotropic stresses, even in simple shear flows. 589][60] One of the most promising ideas has been proposed by Ling et al. 33 with a special neural network architecture referred to as a Tensor Basis Neural Network (TBNN).This architecture, sketched in Figure 11, satisfies the Galilean invariance by relying on the decomposition of the deviatoric stress tensor as a function of a 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 models development.
The TBNN approach developed by Ling et al. 33 was inspired by the work of Pope 61 , 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: . . . . . . . . . . . . . . . . . .

Invariant input layer
Hidden layers Last hidden layer
where a is the deviatoric stress tensor normalized by the turbulent kinetic energy, while λ 1 ,... λ 5 (the scalar basis) and T (1) ...T (10) (the tensor basis) are functions of the filtered strainand rotation-rate tensors.As mentioned by Ling et al. 33 , any tensor who can be expressed as Eq. ( 38) will satisfy Galilean invariance by construction.The simple eddy-viscosity model is recovered by limiting 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 61 , the tensor and scalar bases can be expressed as: T (2) = Ss Rs − Rs Ss T (7) = Rs Ss R2 s − R2 s Ss Rs and where Ss and Rs 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. 33 and Pope 61 , the tensor a is identified to the deviatoric part of the filtered solid stress scaled by the meso-scale pressure (or, equivalently, the sub-grid turbulent kinetic energy): The input tensors Ss and Rs are scaled using the time-scale U t g .The scaled deviatoric part of the solid stress tensor and its eigenvalues (ξ 1 > ξ 2 > ξ 3 ) must satisfy the following realizability conditions: 62 Ling et al. 33 suggest to add a post-processing step after the TBNN model to iteratively enforce conditions given by Eq. (42).Beside the TBNN model for the deviatoric part of the stress, a distinct ANN model is needed to predict the meso-scale pressure.Once P s,meso is known, the anisotropic stress τ s can be retrieved from Eq. ( 41) and the full stress tensor σ 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 ∆ is also added as an extra marker to account for the variation of the mesh size.In the present study, we only consider the reference case (Case 1) to examine the merits of TBNN-based stress modeling for fTFM analysis.Additional markers accounting for variation in physical properties (such as Re p , Fr p ) will enter the network (as was done for the drift flux model development) when the model is generalized to cover a wide a range of gas-particle systems.
A Priori Benchmark Results on ANN Solid Subgrid Stress Modeling MLP Models for Meso-Scale 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 meso-scale viscosity and meso-scale pressure.In both cases, the 3-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 meso-scale pressure and viscosity are scaled as follows: In line with the work of Ouyang et al. 34 , we compare two ANN models for µ * s,meso and P * s,meso with different input markers: • A 3-marker model with λ 1 = tr( S2 s ), φs and ∆ as inputs, • A 14-marker model with φs , ∆ and the components of ũs and ∇ũ s as inputs.
Ouyang et al. 34 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 τ s and Ss (see Eq. ( 36)).Nevertheless, Figures 12 and 13 confirm that the more complete 14-marker ANN model fed with the individual components of ũs and ∇ũ s shows reduced scatter with respect to the simpler 3-marker one, both for the meso-scale pressure and viscosity.
Figure 12: Prediction of the scaled meso-scale solid phase viscosity by the 3-marker (left) and 14-marker (right) ANN models for Case 1.
Figure 13: Prediction of the scaled meso-scale solid pressure by the 3-marker (left) and 14marker (right) ANN models for Case 1.

TBNN model for the subgrid stresses
It was shown in previous section that the meso-scale 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 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 6 independent components of τ * s : 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 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 previous section for comparison.First, it comes out that, even though the 14-marker eddy-viscosity model was shown to be more accurate than the 3-marker model to predict the meso-scale 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) , T (3) and the scalars λ 1 and λ 3 , i.e. discarding quantities that involve the rotation-rate tensor Rs .This model shows reduced RMSE with respect to the two eddy-viscosity models.However, models that account for the first four tensors in the basis display much better performance.Among these models, we observe that using only the first three scalars λ 1 to λ 3 defined by Eq. ( 40) appear sufficient.The last two scalars λ 4 and λ 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 input marker only marginally reduces the RMSE value.Finally, using the complete tensor basis (T (1) to 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 a practical use in fTFM simulations, we suggest to use the TBNN model with the first four basis tensors (T (1) to T (4) ), and, as an input to the network, the first three scalars of Pope's basis (λ 1 to λ 3 ), the filter size ∆ and the filtered solid volume fraction φ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 τ * s and for the square of its norm τ * s : τ * s .It appears that the scatter increases substantially in this second case., T (3) λ 1 , λ 3 , φs , ∆ 0.2076 TBNN model -T (1) to T (4) λ 1 to λ 3 , φs , ∆ 0.0697 TBNN model -T (1) to T (4) λ 1 to λ 3 , φs , ũslip , ∆ 0.0695 TBNN model -T (1) to T (4) λ 1 to λ 5 , φs , ∆ 0.0697 TBNN model -T (1) to T (10) λ 1 to λ 5 , φs , ∆ 0.0692 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 datasets are as follows: The code reads the filtered dataset 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 Figure 14: Prediction of the TBNN model built with four basis tensors (T (1) to T (4) ) and 5 input markers (λ 1 , λ 2 , λ 3 , φs , ∆) for the 6 independent components of the solid phase deviatoric subgrid stress tensor (left) and the square of its norm (right) in Case 1.
dataset of filtered simulation results generated by neptune_cfd , 46 we have uploaded about 3% 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: 63 1) a JSON file that contains the neural network structure; 2) a 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, which 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 ANN models.With ANN models and implementation approach, it is possible to assess its accuracy in a posteriori simulations.

Conclusion
Jiang et al. 23 , 24 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 to the drag force for fTFM simulations, which are feasible for industrial-scale gas-solid flows.
Rauchenzauner and Schneiderbauer 21 used the filtered particle volume fraction, the scalar variance of the sub-grid 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 dataset 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 to the drag force needed for fTFM models.It was found that the approach adopted by Jiang et al. 24 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 datasets.
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 a closure for 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 the 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.

Figure 1 :
Figure 1: Instantaneous solid volume fraction ϕ s in a TFM simulation of a tri-periodic fluidized bed in Case 1 (see Table1).

1 .
develop an ANN model for the filtered drag force (FD-ANN): the target quantity y is the filtered drag Ī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 φ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 different aforementioned strategies for the same level of complexity.The current ANN is made of 3 hidden layers of respectively 128, 32 and 8 nodes, with a ReLU (Rectified Linear Unit) activation function.The output layer has a linear activation function to return the regression result.The loss function used to trained the network is the mean absolute error (MAE), defined as

g d 3 p g µ 2 g
23 and Rauchenzauner and Schneiderbauer43 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 non-uniqueness 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 correctionJiang 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 to 0.4.24The 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 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 a 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 = (ρ s − ρ g )ρ as an additional input variable.They reported a model that employedthe Reynolds number as the extra input, which we consider here.Their drift flux ANN model takes the form:

Figure 3 :
Figure 3: Drift velocity predicted by the ANN of Jiang et al. 24 , compared to the filtered fine-grid simulation data produced by a different research group 43 for a dimensionless filter size ∆f = 9.The filter-size was made dimensionless using the characteristic length-scale L c = d p Fr 1/3p .55

Figure 8 8 PDF 3 -Figure 8 :
Figure 8 shows the PDF of the relative error on 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 that 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.
).They could not discriminate between the three different definitions of L c usually found in the literature: L c,I = d p Fr 1Fr p and L c,III = d p , and the authors set L c = L c,I .We start by adopting the same characteristic length in our analysis.Figure9shows the predictions of the generalized DF-ANN model described by Eq. (28) with L c = L c,I .The training has been performed using 80% of the data points from Cases 1 to 8 while the remaining 20% have been preserved

Figure 9 :
Figure 9: Assessment of the generalized DF-ANN model prediction with characteristic length scale L c = L c,I : testing on Cases 1 to 8 (left) and on Case 9 (right)

Figure 10 :
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

Future work should strive
to generate a comprehensive drift flux model that combines the datasets 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.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.

Table 1 :
Physical and flow parameters of the fine-grid TFM simulations of a gas-solid unbounded fluidized bed in a tri-periodic domain.dp (µm) ρ s (kg m −3 ) ⟨ϕ s ⟩ Fr p Re p Fr ∆

Table 2 :
Comparison between linear eddy-viscosity and various TBNN models to predict the solid phase deviatoric subgrid stress tensor.