# Physics-informed neural network for modelling force and torque fluctuations
in a random array of bidisperse spheres

Zihao Cheng  Anthony Wachs<EMAIL_ADDRESS>Department of Mathematics,
University of British Columbia, Vancouver, BC, Canada.  Department of
Mechanical Engineering, University of British Columbia, Vancouver, BC, Canada.
Department of Chemical & Biological Engineering, University of British
Columbia, Vancouver, BC, Canada.

###### Abstract

We present a physics-informed neural network (PINN) model to predict the
hydrodynamic force and torque fluctuations in a random array of stationary
bidisperse spheres. The PINN model is formulated based on two hypotheses: (i)
pairwise interaction assumption that approximates the total force/torque
exerted on a target sphere by linear superposition of individual contributions
from a finite number of influential neighbors; (i) unified function
representation that suggests a single functional form to describe the
contribution from different neighbors based on the observation of probability
distribution maps obtained with various binary interaction modes in bidisperse
particle-laden flows. On this basis, we accordingly establish a compact PINN
architecture to evaluate individual force/torque contribution of influential
neighbors through the same neural network block which tremendously reduces the
number of unknown parameters, and ultimately compute the total force/torque
exerted on target sphere by their weighted sum. We compare the model
predictions to Particle-Resolved Direct Numerical Simulation (PR-DNS) data of
eight different cases in a range of Reynolds number $1\leq{\mathcal{R}e}\leq
100$, solid volume fraction $10\%\leq\phi\leq 40\%$, sphere diameter ratio
$1.5\leq d_{l}^{*}/d_{s}^{*}\leq 2.5$ and volume ratio $1\leq
V_{l}^{*}/V_{s}^{*}\leq 4$, which demonstrates excellent performance with an
optimal $R^{2}\approx 0.9$ for both force and torque predictions. We establish
a universal model that is applicable within the aforementioned input space,
and examine its interpolation capability to the unseen data with multiple
additional datasets. Finally, we extract and illustrate the interpretable
information of our PINN model in binary and trinary interactions, and discuss
its potential extensions to other particle-laden flow problems with more
complicated scenarios and Eulerian-Lagrangian simulations as a superior
alternative to the classic average drag laws.

###### keywords:

Physics-informed Neural network; Bidisperse particle-laden flows; Force and
torque fluctuations; Pairwise interaction assumption; Unified function
representation

††journal: International Journal of Multiphase Flow

##  1 Introduction

A great portion of numerical efforts in the multiphase flow community over the
last three decades has been dedicated to particle-laden flows that are
commonly encountered in nature (river bank erosion, volcanic eruption, etc.)
and industry (fluid catalytic cracking, coal combustion, etc.). Such a
dispersed multiphase flow typically involves a massive amount of particles
with distinct properties immersed in a homogeneous fluid. The system scale is
often much larger than the characteristic length of the comparatively tiny
particles. With the rapid development of computational methods in fluid
mechanics, we are now able to complement Eulerian-Eulerian (EE) and Eulerian-
Lagrangian (EL) simulations where the boundary layers are unresolved at
particle length scales (even though the particles can be individually tracked
in EL simulations) with high fidelity Particle-Resolved Direct Numerical
Simulation (PR-DNS) that is capable of accessing the hydrodynamic force/torque
exerted on each individual particle and hence of providing abundant microscale
details to reflect the non-negligible force/torque fluctuations induced by the
local microstructure that are not accessible with non-fully resolved
approaches [1, 2, 3, 4].

Even though we have been equipped with several versatile and robust numerical
methods for large-scale computations and benefitted from the expeditious
advance in computing power, PR-DNS of such a complicated and non-linear
particle-laden flow problem is nevertheless expensive in terms of computing
time and resource consumption. Therefore, to provide a theoretical insight
into the neighbor-induced disturbance in particle-laden flows and avoid
repeating the intrinsically same simulation of any given flow configuration,
particular emphasis has been placed on developing a micro/mesoscale model to
account for the force/torque fluctuations. Some early efforts attempted to
apply modifications based on the classical average drag closures to
incorporate the effect of local microstructure heterogeneity with different
factors, e.g. standard deviation of solid concentration fluctuation [5], local
drift velocity [6], gradient of solid concentration [7] and its angle with the
flow direction [8]. In recent years, a so-called stochastic approach drew much
attention and several formulations of drag distributions have been proposed
based on the statistical analysis of PR-DNS results, with each claiming
satisfactory agreement compared to the simulation data obtained in freely
evolving liquid-solid and gas-solid suspensions [9, 10, 11, 12, 13]. With a
distinct perspective, Akiki et al. proposed a pairwise interaction extended
point-particle (PIEP) model that decomposes the microscale flow around a
target/reference particle into the perturbation flows induced by its
neighbors, and the influence due to each neighbor can be separately evaluated
and linearly superposed to obtain the total force exerted on the target
particle. A similar idea was presented in a recent work of Seyed-Ahmadi and
Wachs [14], in which a microstructure-informed probability-driven point-
particle (MPP) model was established to predict the force/torque fluctuations
based on conditioned probability distributions.

Regardless of the categories that distinguish different micro/mesoscale models
for force/torque fluctuations in particle-laden flows, they are all deduced
from and/or compared with PR-DNS data. Thanks to the collective contribution
of the entire CFD community, we have been accumulating an expanding database
obtained from PR-DNS of many particle-laden flow configurations. Considering
the growing demand for processing big data in our problems of interest and the
unprecedented development of data science, we come forward naturally to seek
the help of machine learning (ML) techniques to tackle the still unresolved
problem of accurate force/torque point-particle models. The utilization of ML
concepts in fluid mechanics can even be dated back to the early 1940s for
turbulence study [15]. We refer the interested reader to the review paper by
Burton et al. [16] and the literature therein for a comprehensive
understanding of various branches of ML-related works in flow mechanics. Among
the extensive applications, physics-informed learning, especially the family
of physics-informed neural networks (PINN) pioneered by Raissi et al. [17, 18,
19, 20] has stirred a lively discussion in the scientific computing community.
The conventional ML approaches that leverage a purely data-driven fitting
process may not be able to understand or exploit the interpretable information
and physical patterns underlying the data deluge, thus leading to inconsistent
predictions or even poor performance. On the contrary, the physics-informed
learning machines take into account the fundamental physics that governs the
observation data, e.g., symmetries, invariance and conservation principles,
and are hence capable of fulfilling robust generalization performance even
with sparse and noisy data of complex physical systems. Therefore, an
increasing number of applications are inspired to incorporate the existing
knowledge and constraints to the corresponding learning machines for different
real-world problems [21, 22, 23, 24]. Apart from the original PINN model
proposed by Raissi et al. which successfully reproduced the hydrodynamics of
the flow past a fixed cylinder and vortex-induced vibrations in the case of a
moving cylinder [20, 19], numerous PINN variants have been developed to
increase the overall accuracy and efficiency to solve the Navier-Stokes
equations, e.g., cPINN with domain decomposition technique [25], hp-VPINN
based on the sub-domain Petrov-Galekin method [26], can-PINN with coupled
automatic differentiation (AD) and numerical differentiation (ND) methods [27]
and NSFnets [28]. The interested reader is referred to [29] for a non-
exhaustive list of various PINN architectures to model assorted flow
configurations.

Back to our problem of particle-laden flows, the mission is not to find an
alternative to the conventional CFD methods to solve the forward and inverse
problems involving various partial differential equations, but to utilize a
physics-informed data-driven approach to find the statistical relationship
between hydrodynamic force/torque and different independent variables
throughout innumerable simulation results, and presumably substitute this
relationship in the form of an advanced model into EL simulations with
satisfactory accuracy and efficiency. This has been an active research topic
in the last five years, about which we can find a variety of related works in
the literature, e.g., drag predictions of 2D [30] and 3D non-spherical
particles [31] in both compressible and incompressible flows, rough
surfaces[32] and translational particles in viscoelastic fluids [33]. Several
attempts of leveraging the power of ML have been devoted to multi-particle
systems as well. He and Tafti [34] considered the relative positions of up to
15 nearest neighbors together with the Reynolds number ${\mathcal{R}e}$ and
the solid volume fraction $\phi$, which are the only independent variables in
classic average drag laws, to improve the drag predictions in dense fluid-
particle flows. Muralidhar et al. [35, 36] introduced a PhyDNN model that
demonstrated a better performance by learning with several physical
intermediate variables under the same consideration of 15 nearest neighbors as
inputs. Moore et al. [37] extended the physical PIEP model to a hybrid version
that combines it with a data-driven approach to improve the model performance
at higher $\phi$. Seyed-Ahmadi and Wachs [38] proposed a new PINN architecture
in which the NN block (of fully connected layers) was shared by all the
influential neighbors based on a unified function representation, thus greatly
reducing the number of model parameters and enabling to take more neighbors
into account without changing the PINN framework. Siddani and Balachandar [39]
presented an equivariant NN and trained their universal model based on
datasets from two sources in a range of $0.25\leq{\mathcal{R}e}\leq 250$ and
$0\leq\phi\leq 0.4$. Furthermore, there exists a few PINN models for
predicting forces exerted on non-spherical particles in dense particle-laden
flows, e.g., ellipsoidal particles in 2D [40] and 3D [41].

Even though we are embracing a new era of physics-informed ML with an ever-
increasing number of its applications in the fluid mechanics field, including
particle-laden flow problems that are of particular interest to us, it is
notable that the vast majority of those studies have exclusively focused on
flow configurations with monodisperse particles, i.e., particles of the same
size. To the best of our knowledge, only Hwang et al. [42] presented a deep
learning-based approach to model the drag force in polydisperse particle-laden
flows, but only a dilute system with $\phi=0.5\%$ was considered in that work.
Polydisperse particle-laden flows are not uncommon but rather highly prevalent
in the real world, and of great research interest for their rich
hydrodynamics. The existing gap in the extension of ML techniques to dense
polydisperse systems may be attributed to the data paucity as PR-DNS of such
complicated flows requires more computing resources in order to properly
resolve the smallest particles. In fact, without the comprehensive
understanding of force/torque distributions in each and every particle class
(particles in each class have the same size), it is impossible to provide
insightful knowledge to learning machines and optimize the architectures
accordingly. As a preliminary study, we conducted a series of PR-DNS of the
flow past a random array of stationary bidisperse and polydisperse spheres to
investigate the microstructure-induced force/torque fluctuations in these flow
configurations [43]. The simulation results indicated that every component of
the hydrodynamic force and torque exerted on each particle class follows a
Gaussian distribution. We also confirmed that different binary interaction
modes in bidisperse particle-laden flows contain the same qualitative
characteristics based on the comparison of corresponding probability
distribution maps, which ultimately determines the framework of our PINN
model. Consequently, we select the flow past stationary bidisperse spheres as
a starting point of our "proof-of-concept" study. In this work, we establish a
PINN model similar to that proposed by Seyed-Ahmadi and Wachs [38] in terms of
common NN block for all influential neighbors, and account for the effect of
both the local flow and particle microstructure through an element-wise
multiplication operation. The individual contribution of each neighbor is
evaluated separately but processed by the same NN block according to the
pairwise interaction assumption and unified function representation, and
linearly superposed to acquire the total force/torque exerted on a target
sphere. Apart from the great accuracy and generalization performance of the
prediction of force/torque fluctuations in bidisperse particle-laden flows, we
can also see the potential of extending such a compact yet robust PINN model
to more complicated polydisperse spheres to later plug it in EL simulations.

The rest of paper is organized as follows. In Section 2, we describe the flow
configurations and key parameters of the selected PR-DNS cases, followed by
the dataset preparation procedures. In Section 3, we start with the model
formulation based on the physical fundamentals provided by the previous
research and our preliminary numerical findings, then transition to the PINN
architecture used in the present work and corresponding hyper-parameters
tuning process. In Section 4, we demonstrate the performance of two variants
of the model for force/torque predictions, the former variant corresponding to
a PINN model per set of
$\left({\mathcal{R}e},\phi,d_{l}^{*}/d_{s}^{*},V_{l}^{*}/V_{s}^{*}\right)$
values and the later variant corresponding to a universal model trained with
all available datasets. We evaluate the interpolation capability of the
universal model, and explore the interpretability of our PINN models in terms
of binary and trinary interactions. Finally, we share our outlook for future
research in Section 5.

##  2 Dataset construction

###  2.1 Numerical method and PR-DNS cases

To provide a sufficiently large and reliable database for the subsequent
neural network training procedures, we implement a series of PR-DNS of the
flow past a random array of stationary bidisperse spheres with an Immersed
Boundary-lattice Boltzmann (IB-LB) method that was described comprehensively
in our prior work [44] (we refer the interested reader to this work for
details). We implemented the solver on the open-source platform Basilisk [45]
to enable MPI parallel computing on an adaptive octree-grid, which facilitates
the current large-scale simulations on up to 384 processors with satisfactory
scalability. Apart from the high efficiency, our previous numerical results
demonstrate very satisfactory accuracy of the simulations of the flow past a
random array of moderately to strongly bidisperse spheres as well as of
strongly polydisperse spheres. Relying on the same framework, we extend these
bidisperse array simulations and the flow system can be characterized as
follows.

We consider a pressure-driven flow with average velocity $\bm{u_{0}^{*}}$
(superscript $\dagger^{*}$ describes dimensional quantities) around a
stationary sphere assembly of two groups with different diameter $d^{*}$ in a
tri-periodic cubic box with edge length $L^{*}$. The solid volume fraction
$\phi$ of the system can be corresponding expressed as:

$\phi=\frac{V_{s}^{*}+V_{l}^{*}}{\left(L^{*}\right)^{3}}=\frac{\pi\left[N_{s}\left(d_{s}^{*}\right)^{3}+N_{l}\left(d_{l}^{*}\right)^{3}\right]}{6\left(L^{*}\right)^{3}}$
(1)

where the subscripts $\dagger_{s}$ and $\dagger_{l}$ represent the small and
large spheres respectively, and $V^{*}$ and $N$ denote the total volume and
number of spheres in each group. Here we introduce the most commonly used
characteristic length for bidisperse spheres, the Sauter mean diameter:

$\langle
d^{*}\rangle=\frac{N_{s}\left(d_{s}^{*}\right)^{3}+N_{l}\left(d_{l}^{*}\right)^{3}}{N_{s}\left(d_{s}^{*}\right)^{2}+N_{l}\left(d_{l}^{*}\right)^{2}}$
(2)

which leads us to define the particle Reynolds number ${\mathcal{R}e}$ as
follows:

${\mathcal{R}e}=\frac{\left(1-\phi\right)\left|\bm{u_{0}^{*}}\right|\langle
d^{*}\rangle}{\nu^{*}}=\frac{\left|\bm{U^{*}}\right|\langle
d^{*}\rangle}{\nu^{*}}$ (3)

where $\nu^{*}$ is the fluid kinematic viscosity and
$\bm{U^{*}}=\left(1-\phi\right)\bm{u_{0}^{*}}$ is the superficial velocity. We
constantly impose the flow along the $x$-direction in this work, and record
the steady state data after a sufficiently long simulation time.

${\mathcal{R}e}$ | $\phi$ | $d_{l}^{*}/d_{s}^{*}$ | $N_{l}$ | $N_{s}$ | $V_{l}^{*}/V_{s}^{*}$ | NOR
---|---|---|---|---|---|---
1 | 15% | 2 | 192 | 768 | 2 | 20
10 | 15% | 2 | 192 | 768 | 2 | 20
100 | 15% | 2 | 192 | 768 | 2 | 20
1 | 30% | 2 | 384 | 1536 | 2 | 10
10 | 30% | 2 | 384 | 1536 | 2 | 10
100 | 30% | 2 | 384 | 1536 | 2 | 10
10 | 20% | 1.5 | 453 | 1528 | 1 | 10
10 | 40% | 2.5 | 313 | 1222 | 4 | 13
Table 1: Parameters of PR-DNS cases for generating datasets in the present
study. NOR: number of realizations.

Table 1 lists the eight PR-DNS cases selected in this work, which cover
$1\leq{\mathcal{R}e}\leq 100$, $15\%\leq\phi\leq 40\%$, diameter and volume
ratios of bidisperse spheres $1.5\leq d_{l}^{*}/d_{s}^{*}\leq 2.5$ and $1\leq
V_{l}^{*}/V_{s}^{*}\leq 4$, respectively. We design each and every case with
$L^{*}=20d_{s}^{*}$ and seed the flow domain with $O(10^{3})$ spheres. We
perform $10\sim 20$ realizations to generate approximately 20000 data points
(among which around 4000 are related to large spheres) for each case.

###  2.2 Data post-processing

We organize the dataset of each simulation case in the following manner: we
find the $P$ ($P=100$ in this work) nearest neighbors of each sphere and sort
them by their center-to-center distance to the target sphere. The tri-periodic
boundary conditions have been considered during the search procedure. Each
neighbor features four quantities: its diameter and three components of inter-
center distance. Following the neighbor information, we attach six key
parameters of interest to each target sphere $i$: three components of both the
total fluid-particle interaction force $\bm{F}_{i}^{*}$ and torque
$\bm{T}_{i}^{*}$. Furthermore, we include the target sphere diameter
$d_{i}^{*}$, the local volume-averaged fluid velocity
$\overline{\bm{u}}_{i}^{*}$ and the local solid volume fraction $\phi_{i}$ as
additional information to depict the corresponding local microstructure of
each target sphere as suggested in [38, 39]. In this work, we consider as
"local" a cubic box centered at each target sphere with an edge length of
$5d_{s}^{*}$ and $6d_{s}^{*}$ for small and large spheres respectively with
$\phi\leq 20\%$, and $4d_{s}^{*}$ and $5d_{s}^{*}$ otherwise.

Now we collect the data for all spheres and adopt the following scaling
factors for the dimensional quantities:

  * •

length scale: $d_{s}^{*}$

  * •

velocity scale: $\left|\bm{U^{*}}\right|$

  * •

force scale:
$3\pi\mu^{*}d_{i}^{*}\left|\bm{\overline{\bm{u}}_{i}^{*}}\right|$, similar to
the Stokes drag, where $\mu^{*}$ is the fluid dynamic viscosity.

  * •

torque scale:
$\mu^{*}\left(d_{i}^{*}\right)^{2}\left|\bm{\overline{\bm{u}}_{i}^{*}}\right|$

After the aforementioned non-dimensionalization procedures, we constructed
eight datasets corresponding to eight simulation cases. To characteristically
distinguish them, we label each dataset in the following manner
$\left({\mathcal{R}e},\phi,d_{l}^{*}/d_{s}^{*},V_{l}^{*}/V_{s}^{*}\right)$
hereafter (please note that since $d_{s}=d_{s}^{*}/d_{s}^{*}=1$, we can also
use $d_{l}^{*}/d_{s}^{*}=d_{l}$ for conciseness, but prefer the former
notation to emphasize the fact that the third parameter in the label of each
dataset is indeed a diameter ratio).

Figure 1: Comparison of $F_{x}$ distribution of large and small spheres in
each case obtained with PR-DNS.

Figure 1 compares the distribution of streamwise force $F_{x}$ exerted on
bidisperse spheres in each computed case. Please note that $F_{x}$ is the only
force component to have a non-zero mean value, while $F_{y}$, $F_{z}$ and all
torque components have a zero mean value. We remind the reader that we use a
different scaling for the $y$ axis in each case of Figure 1. While the shaded
areas are visually not equal in size, they are all equal to $1$ as expected.
From the plot, we attempt to convey two messages: (i) there exists a certain
pattern of the force and torque distribution in bidisperse particle-laden flow
problems (and also polydisperse cases), i.e., hydrodynamic force and torque
exerted on each particle class both follow a Gaussian distribution, which
constitute a foundation to successfully train a neural network; (ii) the
distribution of each group exhibits different mean values, standard deviations
and other statistical characteristics, and each case shows different extents
of overlapping between two groups of data. These features collectively add
complexity to our problem in comparison to the monodisperse particle-laden
flows, especially when aiming to establish a universal model without the need
to separate the small and large spheres.

##  3 Neural network architecture and training process

###  3.1 Model formulation

The force and torque distributions in a random array of stationary spheres are
dependent on the unique local microstructure surrounding each sphere as well
as the macroscopic fluid properties, hence of great difficulty to identify and
predict. The pairwise interaction assumption, which decomposes the total
force/torque exerted on a target sphere into the individual contribution from
only a few influential neighbors, significantly simplifies the analysis of
force and torque fluctuations in dense particle-laden flows, and paves the way
for the development of several deterministic models, e.g. PIEP [46, 47, 37,
48], MPP [14], PINN [38] and equivariant NN models [39], which all demonstrate
good predictive performance in monodisperse particle-laden flow problems
spanning a range from low to moderate ${\mathcal{R}e}$. In this work, we also
formulate the model based on the pairwise interaction assumption and merely
consider the binary interactions as a starting point.

Inspired by Seyed-Ahmadi and Wachs [14, 38], the force $\bm{F}_{i}$ and torque
$\bm{T}_{i}$ exerted on the target sphere $i$ can be approximated as:

$\displaystyle\bm{F}_{i}$
$\displaystyle\approx\sum_{j=1}^{M}\Delta\bm{F}_{j\rightarrow
i}\approx\sum_{j=1}^{M}\alpha_{j}f_{j\rightarrow
i}\left(\left\\{\bm{r}_{j},d_{j}\right\\},\left\\{\bm{\mathrm{{\mathcal{R}e}}}_{i},\phi_{i},d_{i}\right\\}\right)$
(4a) $\displaystyle\bm{T}_{i}$
$\displaystyle\approx\sum_{j=1}^{M}\Delta\bm{T}_{j\rightarrow
i}\approx\sum_{j=1}^{M}\beta_{j}g_{j\rightarrow
i}\left(\left\\{\bm{r}_{j},d_{j}\right\\},\left\\{\bm{\mathrm{{\mathcal{R}e}}}_{i},\phi_{i},d_{i}\right\\}\right)$
(4b)

where $\Delta\bm{F}_{j\rightarrow i}$ and $\Delta\bm{T}_{j\rightarrow i}$ are
the binary interaction force and torque induced by the $j$th neighbor with
diameter $d_{j}$ and distance $\bm{r}_{j}$ to the target sphere $i$. $M$
denotes the number of influential neighbors considered in this work.
$\bm{\mathrm{{\mathcal{R}e}}}_{i}=d_{i}\overline{\bm{u}}_{i}/\nu$ denotes the
local ${\mathcal{R}e}$ of the mesoscale flow around target sphere $i$. It
should be noted that $\Delta\bm{F}_{j\rightarrow i}$ and
$\Delta\bm{T}_{j\rightarrow i}$ are typically interpreted as the perturbation
influences of neighboring particles, thus for the streamwise force $F_{x}$
that has a non-zero mean value, an additional term due to the mesoscale flow
should be added, i.e.,
$\bm{F}_{i}\approx\langle\bm{F}\rangle\left(\bm{\mathrm{{\mathcal{R}e}}}_{i},\phi_{i},d_{i}\right)+\sum_{j=1}^{M}\Delta\bm{F}_{j\rightarrow
i}$. However, we remove the mesoscale contribution term in Equation 4a and
maintain a universal form for all components of both force and torque. In this
regard, $\Delta\bm{F}_{j\rightarrow i}$ and $\Delta\bm{T}_{j\rightarrow i}$
can be interpreted as the contribution of $j$th neighbor through the "local"
microstructure as well as the mesoscale flow, which has incorporated the
collective effect of other spheres in this local volume in an average sense.
Back to Equation 4, this means we evenly distribute the unary term
$\langle\bm{F}\rangle$ to each binary interaction term
$\Delta\bm{F}_{j\rightarrow i}$.

Our previous simulation results [43] provide us with another insight to
further simplify this formulation: the probability distribution maps of
different binary interaction modes (small-to-small, large-to-small, small-to-
large and large-to-large spheres) exhibit identical qualitative features in
terms of the critical regions of neighboring sphere that result into a certain
condition on the force/torque exerted on the target sphere. Therefore, we are
able to use a unified function representation for contributions from different
neighbors as suggested by Seyed-Ahmadi and Wachs [38] and consequently
substitute the giant fully connected NN by a compact share-block NN (we will
present this NN architecture in the next subsection). The unified functions
for binary force and torque in this work, $f_{j\rightarrow i}$ and
$g_{j\rightarrow i}$, are slightly more complicated than those in [38] as they
are not only dependent on the neighbor properties
$\left\\{\bm{r}_{j},d_{j}\right\\}$, but also on the mesoscale flow and local
microstructure properties
$\left\\{\bm{\mathrm{{\mathcal{R}e}}}_{i},\phi_{i},d_{i}\right\\}$. This is
also reflected in the design of our NN architecture. As the final operation,
the total force and torque exerted on a target sphere are evaluated through
the weighted sums with $\alpha_{j}$ and $\beta_{j}$ based on the ordering of
neighbors.

We would like to point out that we exclude the local volume ratio
$V_{i,l}^{*}/V_{i,s}^{*}$ when describing the local microstructure while it
would seem natural to include it. Actually, the global volume ratio
$V_{l}^{*}/V_{s}^{*}$ was unexpectedly ignored in most average drag closures
despite the fact that it has a non-negligible impact on the forces and torques
exerted on bidispere spheres. To justify our choice of local quantities in the
input space, we compared the predictions provided by models trained with and
without $V_{i,l}^{*}/V_{i,s}^{*}$ as one input parameter, and it turned out
that the model performance constantly reduces by a small amount, at most
$2\%$, when taking $V_{i,l}^{*}/V_{i,s}^{*}$ into account. We consider this
slight model performance reduction as insignificant, although it is partially
unexpected. One plausible explanation is that the values of $V_{i,l}^{*}$ and
$V_{i,s}^{*}$ can be calculated by $\phi_{i}$ and
$\left\\{\bm{r}_{j},d_{j}\right\\}$ of neighbors, thus
$V_{i,l}^{*}/V_{i,s}^{*}$ introduces repeated information which eventually
leads to predictions that are slightly less satisfactory. It is also possible
that we fully exploited the potential of our PINN models with the current low
data density, and we can reasonably expect improved performance with
$V_{i,l}^{*}/V_{i,s}^{*}$ as an input if we access to a substantially larger
number of data points. Consequently, in the present work, we consider
$\left\\{\bm{\mathrm{{\mathcal{R}e}}}_{i},\phi_{i},d_{i}\right\\}$ only,
instead of
$\left\\{\bm{\mathrm{{\mathcal{R}e}}}_{i},\phi_{i},d_{i},V_{i,l}^{*}/V_{i,s}^{*}\right\\}$,
to characterize the mesoscale flow and local microstructure around target
sphere $i$.

Finally, we consider the same flow symmetry of pairwise interaction as
introduced in [37, 38]. We refer our reader to these works and a graphical
illustration can be found in Figure 1 of [38].

###  3.2 Neural network architecture

Figure 2 illustrates the NN architecture developed in this work. We can see in
figure 2 that the network contains three main modules: an input layer, two
blocks of fully connected layers and an output layer, each consisting of
several neurons which are the fundamental units in a NN to process the input
signals through a so-called activation function and assist the NN to learn the
non-linearity after multiple-layer transformations. Inspired by Siddani and
Balachandar [39], we introduce an intermediate step of element-wise
multiplication to incorporate the local flow effects, and a same projection
step for the lateral components of force and torque as adopted by Seyed-Ahmadi
and Wachs [38].

Figure 2: Physics-informed neural network architecture used in the present
work.

Now we elaborate on each module and intermediate step: for the input layer, we
need the properties of both the target sphere and its $M$ influential
neighbors, i.e.
$\left\\{\bm{\mathrm{{\mathcal{R}e}}}_{i},\phi_{i},d_{i}\right\\}$ and
$\left\\{\bm{r}_{j},d_{j}\right\\}$ respectively as mentioned above. Each
subset of input parameters is then fed into an individual block of fully
connected layers (also called hidden layers) with the same width, i.e., the
same number of neurons in one layer. Please note that as suggested by the
unified function representation assumption, information of different neighbors
is passed through the same block at this stage, thus significantly reducing
the width (and presumably the length, i.e., the number of layers) of fully
connected layers and hence the unknown parameters. Here we adopt the
exponential linear unit (elu) [49] as the activation function of fully
connected layers which leads to fast and precise learning of NN and avoids the
vanishing gradient problem during back-propagation. Next, the outputs
generated by two individual blocks experience an element-wise multiplication
step to combine the contribution from each neighbor with the effects of local
flow and microstructure of a target sphere. According to the flow symmetry,
each scalar value that corresponds to the unweighted contribution from each
neighbor is multiplied by appropriate direction vectors and then projected on
$y$ and $z$-axis to obtain the lateral force and torque. Eventually, all
scalars are processed by a linear unit (lu) at the output layer to perform a
weighted sum that yields the total force and torque on a target sphere.

It is also worth mentioning that due to the features of our problem, we select
a special loss function, the Huber loss to optimize the values of unknown
parameters (weights and bias in each neuron). In bidisperse particle-laden
flows, the streamwise forces exerted on large spheres are significantly larger
than those of small spheres, whereas the number of larger spheres is much less
than its counterpart. Therefore, the data points corresponding to large
spheres are prone to be treated as "outliers" during the regression process.
Our NN may neglect these "outliers" with a simple MAE (mean average error)
function, or be trapped by them and predict poorly for the small spheres with
a conventional MSE (mean square error) function. The Huber loss balances these
two functions in the following form:

$\mathcal{L}\left(\bm{w},\bm{b}\right)=\begin{cases}\frac{1}{2}\left(\bm{F}_{i,PINN}-\bm{F}_{i,DNS}\right)^{2}&\left|\bm{F}_{i,PINN}-\bm{F}_{i,DNS}\right|\leq\delta\\\
\delta\left(\left|\bm{F}_{i,PINN}-\bm{F}_{i,DNS}\right|-\frac{1}{2}\delta\right)&otherwise\end{cases}$
(5)

where $\bm{w}$ and $\bm{b}$ are the weight and bias matrices respectively. In
the present work, $\delta$ equals to the standard deviation of the
distribution of force/torque exerted on large spheres.

###  3.3 Hyper-parameters and implementation

We presented the NN architecture in the previous subsection, while we still
need to determine the values of several hyper-parameters, that significantly
affect the learning process and performance of our NN and need to be fine-
tuned in advance. Before introducing the hyper-parameter tuning process in
this work, we would like to briefly present two notions: (i) $5$-fold cross-
validation: we randomly shuffle the dataset first and then split it into 4
folds of training set and 1 fold of test set for regression and evaluation
purposes respectively. Each fold of data is retained as test set in turn to
evaluate the model performance on the rest folds after each round of training,
and the ultimate performance is averaged over all 5 rounds. In our
implementation, 20% of the training set is additionally separated and used as
validation set during the training process to validate the model
generalization. (ii) coefficient of determination $R^{2}$: a metric to
quantify the model performance, which typically has the following form:

$R^{2}=1-\frac{\sum_{i=1}^{N}\left(\bm{F}_{i,PINN}-\bm{F}_{i,DNS}\right)^{2}}{\sum_{i=1}^{N}\left(\bm{F}_{i,DNS}-\langle\bm{F}_{i,DNS}\rangle\right)^{2}}$
(6)

where $\langle\bm{F}_{i,DNS}\rangle$ is the average force/torque computed by
PR-DNS. Considering the problem of interest in this work, we use the following
formula to evaluate the overall performance:

$R^{2}=1-\frac{\sum_{i=1}^{N_{s}}\left(\bm{F}_{(s,i),PINN}-\bm{F}_{(s,i),DNS}\right)^{2}+\sum_{i=1}^{N_{l}}\left(\bm{F}_{(l,i),PINN}-\bm{F}_{(l,i),DNS}\right)^{2}}{\sum_{i=1}^{N_{s}}\left(\bm{F}_{(s,i),DNS}-\langle\bm{F}_{(s,i),DNS}\rangle\right)^{2}+\sum_{i=1}^{N_{l}}\left(\bm{F}_{(l,i),DNS}-\langle\bm{F}_{(l,i),DNS}\rangle\right)^{2}}$
(7)

where we separate the force/torque exerted on small and large spheres. Please
note that without special indication, the $R^{2}$ scores reported hereafter
are evaluated on test sets through 5-fold cross-validation.

We begin with the number of influential neighbors to be included in the input
space. Neighbor-truncation error analysis of both [38] and [39] shows that the
$R^{2}$ converges with $25\sim 30$ neighbors for force and 10 for torque
predictions in monodisperse particle-laden flows. We carry out a similar
convergence test for the streamwise force $F_{x}$, lateral force $F_{L}$ and
lateral torque $T_{\perp}$ in bidisperse particle-laden flows and compare the
training (hollow bar) and test (solid bar) $R^{2}$ scores for each component
of both sphere classes with a different number of neighbors $M$ in Figure 3.
Once again, we gain two observations from the plot: (i) the training and test
$R^{2}$ are highly close to each other in every case regardless of $M$,
demonstrating great generalization of our NN model, and (ii) $R^{2}$ reaches a
plateau with $M=40\sim 50$ for force and $15\sim 20$ for torque, and
consequently we select $M=50$ for the former and and $20$ for the latter.

Figure 3: $R^{2}$ scores of $F_{x}$, $F_{L}$ and $T_{\perp}$ for the
increasing number of neighbors $M$ included in the PINN model. We select
$M=50$ and $M=20$ for force and torque predictions respectively based on this
convergence study.

We construct and implement our NN model using TensorFlow [50] with Keras API
[51] in Python. The Nadam optimizer [52], which is an extension to the Adam
optimizer with Nesterov momentum, is adopted in our model training process.
Two callbacks are applied to enhance the efficiency and accuracy: (i)
ReduceLROnPlateau: the learning rate will be reduced by a factor of 2 if no
obvious improvement is achieved in terms of decreasing the total loss after 10
epochs; (ii) EarlyStopping: the training process will be terminated if no
obvious improvement is achieved after 40 consecutive epochs. Figure 4 exhibits
the evolution of Huber losses of both the training and validation sets as well
as the change in learning rate with increasing epochs. We can see that the
training and validation losses both decrease rapidly at the early stage, and
the adaptive learning rate assists the NN model to escape from the local
minimum and continue to search for the global minimum. All of the loss
functions reach the converged values and no over-fitting occurs during the
training process. In addition, a batch size of 128 is used in all cases, and 3
layers of 30 neurons for force and 20 neurons for torque predictions are
determined for the fully connected layer blocks after a grid-search within a
broad range of hyper-parameters.

Figure 4: Evolution of Huber loss during the model training process. From left
to right: $F_{x}$, $F_{L}$ and $T_{\perp}$.

##  4 Results and discussions

###  4.1 Performance over monodisperse cases

To validate the performance of our PINN model, we start with the simulation
cases of flow past a random array of stationary monodisperse spheres
implemented by Seyed-Ahmadi and Wachs [38]. These simulations cover broad
ranges of $10\%\leq\phi\leq 40\%$ and $0.2\leq{\mathcal{R}e}\leq 150$. We
perform the non-dimensionalization procedures described in Section 2.2 to the
original datasets. In our training process, $M=30$ and 10 are considered for
force and torque predictions which are the same as those used in [38]. Table 2
compares the $R^{2}$ scores obtained with the PINN model proposed in [38] and
the present model for $F_{x}$, $F_{L}$ and $T_{\perp}$ (denoted $T_{L}$ in
[38]) respectively. We can see from the table that almost all $R^{2}$ have
increased with the present model except the $R^{2}$ of $F_{L}$ in the case
$\left({\mathcal{R}e},\phi\right)=\left(150,40\%\right)$, which reaches the
upper limits of both parameters so that more data points should be expected to
acquire a satisfactory prediction for the lateral force that varies more
dramatically with higher ${\mathcal{R}e}$ and $\phi$. Unfortunately, we can
only access approximately the same amount of data points for every case due to
limited computing resources (in fact, the data points for the case
$\left(150,40\%\right)$ are even fewer than those in the other cases [38]).
Nevertheless, our PINN model exhibits superior performance in all the other
predictions, where more than half of them achieve a higher than 5% increase in
$R^{2}$. In addition, our results disclose a clear pattern that higher $R^{2}$
can be achieved at lower ${\mathcal{R}e}$. This improvement is physically
understandable since the problem is less non-linear but not being observed in
the previous results.

$\phi$ | ${\mathcal{R}e}$ | $F_{x}$ | $F_{L}$ | $T_{\perp}$
---|---|---|---|---
[38] | Present | [38] | Present | [38] | Present
10% | 2 | 0.81 | 0.88 | 0.77 | 0.84 | 0.89 | 0.90
10 | 0.82 | 0.84 | 0.78 | 0.80 | 0.86 | 0.87
40 | 0.70 | 0.72 | 0.73 | 0.77 | 0.74 | 0.81
20% | 0.2 | 0.68 | 0.79 | 0.57 | 0.78 | 0.74 | 0.89
2 | 0.73 | 0.78 | 0.68 | 0.76 | 0.88 | 0.90
40 | 0.71 | 0.74 | 0.67 | 0.72 | 0.79 | 0.81
150 | 0.61 | 0.63 | 0.56 | 0.65 | 0.56 | 0.71
40% | 2 | 0.61 | 0.79 | 0.51 | 0.61 | 0.72 | 0.77
40 | 0.66 | 0.75 | 0.58 | 0.59 | 0.69 | 0.73
150 | 0.52 | 0.69 | 0.53 | 0.51 | 0.59 | 0.62
Table 2: $R^{2}$ score comparison of the PINN model proposed in [38] and the
present model for monodisperse particle-laden flows.

###  4.2 Performance over bidisperse cases

Now we turn our attention to our problem of interest: flow past a random array
of stationary bidisperse spheres described in Section 2.1. Figure 5 shows the
correlation plot of $F_{x}$, $F_{L}$ and $T_{\perp}$ exerted on all spheres
obtained with our PR-DNS and the PINN model predictions in the two cases
$\left({\mathcal{R}e},\phi,d_{l}^{*}/d_{s}^{*},V_{l}^{*}/V_{s}^{*}\right)=\left(1,15\%,2,2\right)$
and $\left(100,30\%,2,2\right)$. These two cases correspond to the lower and
upper limits of the two key parameters ${\mathcal{R}e}$ and $\phi$ where we
indeed observe the best and worst model performance. To distinguish two sphere
classes, we plot the large sphere data in red and the small sphere data in
blue. On the $F_{x}$ plot where the large and small spheres have non-zero mean
values, we attach the probability distributions of both groups aside for
further comparison. The quantitative results, $R^{2}$ scores of each group
obtained with Equation 6, are also added to the corresponding subplot. We
observe that both the red and blue scatters are located closely around the
diagonal dashed line that represents the perfect prediction (i.e., $R^{2}=1$),
which reflects that the error distribution is not tremendously skewed to one
side and our PINN model gives ideal predictions in most cases. We also note
that the data bands are much narrower at low ${\mathcal{R}e}$ and low $\phi$
than at large ${\mathcal{R}e}$ and large $\phi$. The probability distributions
of $F_{x}$ obtained with PR-DNS and PINN model are remarkably similar in shape
and range, and visually symmetric about the dashed line in both cases.
Moreover, the converged shape at the two ends of each data band in all
subplots indicates that our PINN model possesses good learning
ability/generalization and is capable of providing acceptable predictions even
in regions with sparse data.

(a)
$\left({\mathcal{R}e},\phi,d_{l}^{*}/d_{s}^{*},V_{l}^{*}/V_{s}^{*}\right)=\left(1,15\%,2,2\right)$

(b)
$\left({\mathcal{R}e},\phi,d_{l}^{*}/d_{s}^{*},V_{l}^{*}/V_{s}^{*}\right)=\left(100,30\%,2,2\right)$

Figure 5: Comparison of $F_{x}$, $F_{L}$ and $T_{\perp}$ obtained with our PR-
DNS and PINN predictions.

Table 3 summarizes the individual and overall performance of $F_{x}$, $F_{L}$
and $T_{\perp}$ predictions in all eight cases. Several observations can be
made from the comparison: (i) The $R^{2}$ scores of the training and test sets
are close to each other in most cases (difference $<$ 0.03), especially in the
lateral torque predictions that exhibit the most satisfactory model
generalization. (ii) The $R^{2}$ scores of small spheres are generally better
than those of large spheres in force predictions, while they are highly
comparable in torque predictions. The reason is that large spheres typically
experience higher forces than their small counterparts, whereas they only
possess a small proportion in the entire dataset, thus compromising the
performance during the regression process. However, $T_{\perp}$ exerted on
different sphere classes are quite close to each other, and less sensitive to
the change of parameters ${\mathcal{R}e}$, $\phi$ and $d_{l}^{*}/d_{s}^{*}$.
We refer the interested reader to our previous work for more details on
simulation results of force and torque distributions in flows past a random
array of stationary bidisperse spheres, and we will come back to this later in
the statistical analysis. (iii) There exists the following clear pattern: if
we increase any key parameter ${\mathcal{R}e}$, $\phi$ and
$d_{l}^{*}/d_{s}^{*}$ while keeping the two other parameters constant, we will
see a significant performance drop in all three predictions $F_{x}$, $F_{L}$
and $T_{\perp}$. (iv) The $T_{\perp}$ prediction provides the best overall
performance even with a smaller number of influential neighbors than that used
for force predictions, whereas $F_{L}$ produces the worst overall performance.
This observation is in line with the results reported by Seyed-Ahmadi and
Wachs [38]. As a corollary, our PINN model generates great predictions in
bidisperse particle-laden flows and the performance is either at the same
level or even better than that obtained in monodisperse particle-laden flows
in [38] and [39]. The overall $R^{2}$ spans $0.706\sim 0.892$ for $F_{x}$,
$0.669\sim 0.871$ for $F_{L}$ and $0.741\sim 0.910$ for $T_{\perp}$
predictions when $1\leq{\mathcal{R}e}\leq 100$, $10\%\leq\phi\leq 40\%$,
$1.5\leq d_{l}^{*}/d_{s}^{*}\leq 2.5$ and $1\leq V_{l}^{*}/V_{s}^{*}\leq 4$.

$\left({\mathcal{R}e},\phi,d_{l}^{*}/d_{s}^{*},V_{l}^{*}/V_{s}^{*}\right)$ | Large spheres $R^{2}$ | Small spheres $R^{2}$ | Overall $R^{2}$
---|---|---|---
Training sets | Test sets | Training sets | Test sets
$\left(1,15\%,2,2\right)$ | 0.898 | 0.875 | 0.911 | 0.900 | 0.892
$\left(10,15\%,2,2\right)$ | 0.888 | 0.852 | 0.892 | 0.874 | 0.868
$\left(100,15\%,2,2\right)$ | 0.777 | 0.712 | 0.743 | 0.704 | 0.706
$\left(1,30\%,2,2\right)$ | 0.893 | 0.849 | 0.901 | 0.884 | 0.871
$\left(10,30\%,2,2\right)$ | 0.888 | 0.838 | 0.893 | 0.871 | 0.859
$\left(100,30\%,2,2\right)$ | 0.818 | 0.755 | 0.800 | 0.771 | 0.766
$\left(10,20\%,1.5,1\right)$ | 0.909 | 0.858 | 0.894 | 0.869 | 0.864
$\left(10,40\%,2.5,4\right)$ | 0.852 | 0.829 | 0.876 | 0.860 | 0.852
(a) $F_{x}$

$\left({\mathcal{R}e},\phi,d_{l}^{*}/d_{s}^{*},V_{l}^{*}/V_{s}^{*}\right)$ | Large spheres $R^{2}$ | Small spheres $R^{2}$ | Overall $R^{2}$
---|---|---|---
Training sets | Test sets | Training sets | Test sets
$\left(1,15\%,2,2\right)$ | 0.854 | 0.832 | 0.886 | 0.879 | 0.865
$\left(10,15\%,2,2\right)$ | 0.834 | 0.805 | 0.832 | 0.818 | 0.814
$\left(100,15\%,2,2\right)$ | 0.802 | 0.733 | 0.741 | 0.707 | 0.717
$\left(1,30\%,2,2\right)$ | 0.893 | 0.849 | 0.901 | 0.884 | 0.871
$\left(10,30\%,2,2\right)$ | 0.793 | 0.815 | 0.740 | 0.794 | 0.777
$\left(100,30\%,2,2\right)$ | 0.748 | 0.670 | 0.705 | 0.666 | 0.669
$\left(10,20\%,1.5,1\right)$ | 0.791 | 0.768 | 0.806 | 0.794 | 0.786
$\left(10,40\%,2.5,4\right)$ | 0.770 | 0.664 | 0.816 | 0.789 | 0.743
(b) $F_{L}$

$\left({\mathcal{R}e},\phi,d_{l}^{*}/d_{s}^{*},V_{l}^{*}/V_{s}^{*}\right)$ | Large spheres $R^{2}$ | Small spheres $R^{2}$ | Overall $R^{2}$
---|---|---|---
Training sets | Test sets | Training sets | Test sets
$\left(1,15\%,2,2\right)$ | 0.915 | 0.910 | 0.914 | 0.911 | 0.910
$\left(10,15\%,2,2\right)$ | 0.901 | 0.892 | 0.896 | 0.891 | 0.891
$\left(100,15\%,2,2\right)$ | 0.782 | 0.759 | 0.746 | 0.731 | 0.741
$\left(1,30\%,2,2\right)$ | 0.910 | 0.902 | 0.896 | 0.888 | 0.895
$\left(10,30\%,2,2\right)$ | 0.906 | 0.895 | 0.888 | 0.882 | 0.886
$\left(100,30\%,2,2\right)$ | 0.788 | 0.767 | 0.754 | 0.742 | 0.750
$\left(10,20\%,1.5,1\right)$ | 0.904 | 0.896 | 0.896 | 0.892 | 0.893
$\left(10,40\%,2.5,4\right)$ | 0.870 | 0.855 | 0.875 | 0.869 | 0.864
(c) $T_{\perp}$

Table 3: $R^{2}$ scores of $F_{x}$, $F_{L}$ and $T_{\perp}$ predictions for
large and small spheres, as well as the overall performance of each dataset.

Finally, we analyze some statistical features of the $F_{x}$, $F_{L}$ and
$T_{\perp}$ distributions obtained with PR-DNS and with the PINN model. We
first compare the mean value of $F_{x}$ that is the major concern of classic
drag models proposed in many previous works [53, 54, 55, 56, 57, 58, 59, 60,
61]. It can be seen through the comparison that the largest error between the
PR-DNS results and the PINN predictions is merely 0.18%. Such a small
discrepancy is undoubtedly better than that obtained with any existing drag
closure. As mentioned previously, although we are not able to obtain a perfect
prediction on each sphere, the errors are ideally centered around 0 and follow
a Gaussian distribution, which compensates the positive errors with negative
errors after summing over all spheres and ultimately yields extraordinary
predictions of average $F_{x}$ that are close to the ground truth. Then we
move our attention to the standard deviation of all three components. Overall,
the standard deviations of all distributions are accurately predicted by our
PINN model. The ratios of the largest to the smallest standard deviation of
$F_{x}$, $F_{L}$ and $T_{\perp}$ are $7.4\sim 7.5$, $6.3\sim 6.7$ and $2.7\sim
2.8$ respectively, confirming that $T_{\perp}$ is the least sensitive
component to the variation of the key parameters. Last but not the least, we
can observe that all standard deviations obtained with our PINN model are
smaller than that computed by PR-DNS, demonstrating that the present PINN
model is conservative in low data density regimes. This seems to imply that
our model is unsurprisingly well trained for interpolation but not for
extrapolation (please note that this is the common challenge for NN to deal
with the unseen data).

$\left({\mathcal{R}e},\phi,d_{l}^{*}/d_{s}^{*},V_{l}^{*}/V_{s}^{*}\right)$ | Large spheres | Small spheres
---|---|---
$\langle F_{x,DNS}\rangle$ | $\langle F_{x,PINN}\rangle$ | Error | $\langle F_{x,DNS}\rangle$ | $\langle F_{x,PINN}\rangle$ | Error
$\left(1,15\%,2,2\right)$ | 5.585 | 5.588 | 0.05% | 2.918 | 2.914 | 0.12%
$\left(10,15\%,2,2\right)$ | 6.763 | 6.772 | 0.14% | 3.509 | 3.512 | 0.10%
$\left(100,15\%,2,2\right)$ | 15.592 | 15.581 | 0.07% | 7.673 | 7.665 | 0.11%
$\left(1,30\%,2,2\right)$ | 15.771 | 15.756 | 0.10% | 6.712 | 6.711 | 0.00%
$\left(10,30\%,2,2\right)$ | 18.114 | 18.125 | 0.06% | 7.701 | 7.701 | 0.01%
$\left(100,30\%,2,2\right)$ | 37.809 | 37.741 | 0.18% | 15.784 | 15.795 | 0.07%
$\left(10,20\%,1.5,1\right)$ | 8.442 | 8.435 | 0.08% | 5.491 | 5.483 | 0.14%
$\left(10,40\%,2.5,4\right)$ | 34.065 | 34.076 | 0.03% | 10.085 | 10.081 | 0.04%
Table 4: Comparison of the mean values of $F_{x}$ for each case obtained with
our PR-DNS and PINN predictions.

$\left({\mathcal{R}e},\phi,d_{l}^{*}/d_{s}^{*},V_{l}^{*}/V_{s}^{*}\right)$ | $\sigma_{F_{x,DNS}}$ | $\sigma_{F_{x,PINN}}$ | $\sigma_{F_{L,DNS}}$ | $\sigma_{F_{L,PINN}}$ | $\sigma_{T_{\perp,DNS}}$ | $\sigma_{T_{\perp,PINN}}$
---|---|---|---|---|---|---
$\left(1,15\%,2,2\right)$ | 0.946 | 0.889 | 0.697 | 0.640 | 1.525 | 1.457
$\left(10,15\%,2,2\right)$ | 1.224 | 1.139 | 0.867 | 0.792 | 1.575 | 1.491
$\left(100,15\%,2,2\right)$ | 3.314 | 2.901 | 2.207 | 1.970 | 2.029 | 1.789
$\left(1,30\%,2,2\right)$ | 2.502 | 2.352 | 1.652 | 1.473 | 2.900 | 2.767
$\left(10,30\%,2,2\right)$ | 2.885 | 2.716 | 1.976 | 1.755 | 2.974 | 2.823
$\left(100,30\%,2,2\right)$ | 6.971 | 6.228 | 4.648 | 4.014 | 3.940 | 3.491
$\left(10,20\%,1.5,1\right)$ | 1.485 | 1.375 | 1.090 | 0.960 | 1.952 | 1.855
$\left(10,40\%,2.5,4\right)$ | 5.766 | 5.472 | 3.298 | 2.862 | 4.128 | 3.849
(a) Large spheres

$\left({\mathcal{R}e},\phi,d_{l}^{*}/d_{s}^{*},V_{l}^{*}/V_{s}^{*}\right)$ | $\sigma_{F_{x,DNS}}$ | $\sigma_{F_{x,PINN}}$ | $\sigma_{F_{L,DNS}}$ | $\sigma_{F_{L,PINN}}$ | $\sigma_{T_{\perp,DNS}}$ | $\sigma_{T_{\perp,PINN}}$
---|---|---|---|---|---|---
$\left(1,15\%,2,2\right)$ | 0.681 | 0.650 | 0.529 | 0.497 | 1.126 | 1.077
$\left(10,15\%,2,2\right)$ | 0.950 | 0.897 | 0.641 | 0.584 | 1.169 | 1.106
$\left(100,15\%,2,2\right)$ | 2.545 | 2.195 | 1.457 | 1.255 | 1.404 | 1.208
$\left(1,30\%,2,2\right)$ | 1.626 | 1.543 | 1.223 | 1.126 | 2.247 | 2.128
$\left(10,30\%,2,2\right)$ | 2.013 | 1.906 | 1.466 | 1.320 | 2.346 | 2.202
$\left(100,30\%,2,2\right)$ | 5.118 | 4.605 | 3.314 | 2.784 | 2.954 | 2.565
$\left(10,20\%,1.5,1\right)$ | 1.308 | 1.221 | 0.929 | 0.833 | 1.659 | 1.571
$\left(10,40\%,2.5,4\right)$ | 3.000 | 2.841 | 2.170 | 1.956 | 3.122 | 2.920
(b) Small spheres

Table 5: Comparison of the standard deviations of $F_{x}$, $F_{L}$ and
$T_{\perp}$ for each case obtained with our PR-DNS and PINN predictions.

###  4.3 Extension to a universal model

In the previous subsection, we trained and evaluated our PINN model on each
individual dataset as summarized in Table 3. It is natural to wonder if a
single model could provide predictions at the same level of accuracy in a
broad parameter space and would hence not be limited to a specific
configuration, as the ultimate objective of this study is to find a promising
alternative to the average drag closures that are typically applicable over a
range of ${\mathcal{R}e}$, $\phi$, $d_{l}^{*}/d_{s}^{*}$ and
$V_{l}^{*}/V_{s}^{*}$. Therefore, we implement the training procedures for the
$F_{x}$ prediction using the same hyper-parameters as those optimized in
Section 3.3 on multiple integrated datasets, and run additional simulations
with distinct sets of parameters that fall within the input space described in
Table 1. We feed the new datasets into the updated and so called universal
PINN model and evaluate the corresponding $R^{2}$ scores, thus examining the
interpolation capability of the present model within a given parameter regime.

(a) Interpolation of ${\mathcal{R}e}$

(b) Interpolation of $\phi$

(c) Interpolation of all particle properties

Figure 6: Prediction on the unseen data (green dots) using the PINN models
trained on integrated datasets (red and blue dots for large and small spheres
respectively).

$\left({\mathcal{R}e},\phi,d_{l}^{*}/d_{s}^{*},V_{l}^{*}/V_{s}^{*}\right)$ | Large spheres $R^{2}$ | Small spheres $R^{2}$
---|---|---
Training sets | Test sets | Training sets | Test sets
$\left(1,15\%,2,2\right)$ | 0.878 | 0.862 | 0.895 | 0.890
$\left(10,15\%,2,2\right)$ | 0.859 | 0.842 | 0.862 | 0.865
$\left(100,15\%,2,2\right)$ | 0.754 | 0.719 | 0.721 | 0.704
1-5 $\left(5,15\%,2,2\right)$ | - | 0.764 | - | 0.805
(a) Interpolation of ${\mathcal{R}e}$

$\left({\mathcal{R}e},\phi,d_{l}^{*}/d_{s}^{*},V_{l}^{*}/V_{s}^{*}\right)$ | Large spheres $R^{2}$ | Small spheres $R^{2}$
---|---|---
Training sets | Test sets | Training sets | Test sets
$\left(1,15\%,2,2\right)$ | 0.882 | 0.869 | 0.895 | 0.895
$\left(1,30\%,2,2\right)$ | 0.880 | 0.853 | 0.895 | 0.887
1-5 $\left(1,22.5\%,2,2\right)$ | - | 0.535 | - | 0.603
(b) Interpolation of $\phi$

$\left({\mathcal{R}e},\phi,d_{l}^{*}/d_{s}^{*},V_{l}^{*}/V_{s}^{*}\right)$ | Large spheres $R^{2}$ | Small spheres $R^{2}$
---|---|---
Training sets | Test sets | Training sets | Test sets
$\left(10,15\%,2,2\right)$ | 0.852 | 0.846 | 0.862 | 0.863
$\left(10,20\%,1.5,1\right)$ | 0.835 | 0.831 | 0.863 | 0.862
$\left(10,30\%,2,2\right)$ | 0.857 | 0.848 | 0.881 | 0.878
$\left(10,40\%,2.5,4\right)$ | 0.893 | 0.867 | 0.886 | 0.873
1-5 $\left(10,25\%,1.8,1.5\right)$ | - | 0.833 | - | 0.862
(c) Interpolation of all particle properties

Table 6: $R^{2}$ scores of $F_{x}$ predictions for large and small spheres on
both training and unseen datasets.

There are four changing parameters in our input space: flow property
${\mathcal{R}e}$ and particle properties $\phi$, $d_{l}^{*}/d_{s}^{*}$ and
$V_{l}^{*}/V_{s}^{*}$. We start with the simplest scenario of interpolating in
a single direction while fixing all the other parameters, and subsequently
proceed to a more challenging case of interpolating all particle properties:

  1. 1.

Interpolation of ${\mathcal{R}e}$: we train the PINN model on three datasets
with the same particle properties
$\left(\phi,d_{l}^{*}/d_{s}^{*},V_{l}^{*}/V_{s}^{*}\right)=\left(15\%,2,2\right)$
at ${\mathcal{R}e}=1$, 10 and 100, and generate an additional dataset with the
same particle properties at ${\mathcal{R}e}=5$. We exhibit the correlation
plot of all PR-DNS data and PINN predictions in Figure 6 where the large and
small particles with seen parameters (i.e., ${\mathcal{R}e}=1$, 10 and 100
that are used for training the model) are displayed in red and blue
respectively, and predictions with unseen parameter (i.e., ${\mathcal{R}e}=5$
that does not contribute to the training process) are represented by green
dots. We can see from Figure 6(a) that these green dots are distributed around
the black diagonal line that corresponds to $R^{2}=1$, indicating good
performance in predicting unseen data at ${\mathcal{R}e}=5$. The quantitative
results are presented in Table 6 where similar $R^{2}$ scores can be observed
for training datasets as those reported in Table 3 produced by the separate
models, and good predictions on the unseen dataset can be also claimed with
$R^{2}=0.764$ and $0.805$ for large and small spheres respectively.

  2. 2.

Interpolation of $\phi$: we train the PINN model on two datasets with the same
$\left({\mathcal{R}e},d_{l}^{*}/d_{s}^{*},V_{l}^{*}/V_{s}^{*}\right)=\left(1,2,2\right)$
but different $\phi=15\%$ and 30%, and generate an extra dataset at
$\phi=22.5\%$. Figure 6(b) reveals that green dots (predictions on unseen
data) are loosely distributed about the diagonal line, implying a higher level
of prediction error in general, with overestimation occurring most frequently
for those representing larger $F_{x}$. Quantitatively, the $R^{2}$ decreases
to $0.535$ for large spheres and $0.603$ for small spheres when only two
datasets (lower and upper bounds with a broad interval) are considered when
training the model, which is a compromised choice due to the limited number of
data available in the present study.

  3. 3.

Interpolation of all particle properties: we train the PINN model on four
datasets with the same ${\mathcal{R}e}=10$ but various particle properties,
and then generate a new dataset with
$\left(\phi,d_{l}^{*}/d_{s}^{*},V_{l}^{*}/V_{s}^{*}\right)=\left(25\%,1.8,1.5\right)$
such that all three parameters lie within the input ranges and differ to any
original inputs (please note that we add $V_{l}^{*}/V_{s}^{*}$ to the input
space to train this model). Figure 6(c) demonstrates astonishingly accurate
predictions on the unseen data and the corresponding $R^{2}=0.833$ and $0.862$
can be found in Table 6(c) for large and small spheres respectively, which are
even very close to the $R^{2}$ scores of the training datasets. One plausible
explanation is that the increments in the input parameters of these four
datasets are smaller than those of previous two tests
(${\mathcal{R}e}=1\rightarrow 10\rightarrow 100$ and $\phi=15\%\rightarrow
30\%$), resulting in a more uniform distribution of training data that leaves
fewer and smaller unexplored regions in the parameter space as shown in Figure
6(c).

(a) Universal model

(b) Predictions on unseen data

Figure 7: Training of a universal model. (a) Comparison of $F_{x}$ obtained with our PR-DNS and universal PINN predictions. Left: large spheres, and right: small spheres. (b) Interpolation of all parameters. $\left({\mathcal{R}e},\phi,d_{l}^{*}/d_{s}^{*},V_{l}^{*}/V_{s}^{*}\right)$ | Large spheres $R^{2}$ | Small spheres $R^{2}$
---|---|---
Training sets | Test sets | Training sets | Test sets
$\cellcolor{red!40}\left(1,15\%,2,2\right)$ | 0.868 | 0.859 | 0.885 | 0.880
$\cellcolor{orange!50}\left(10,15\%,2,2\right)$ | 0.861 | 0.851 | 0.865 | 0.867
$\cellcolor{yellow!40}\left(10,20\%,1.5,1\right)$ | 0.838 | 0.832 | 0.863 | 0.862
$\cellcolor{green!40}\left(1,30\%,2,2\right)$ | 0.870 | 0.858 | 0.891 | 0.887
$\cellcolor{cyan!50}\left(10,30\%,2,2\right)$ | 0.863 | 0.850 | 0.882 | 0.879
$\cellcolor{blue!30}\left(10,40\%,2.5,4\right)$ | 0.894 | 0.874 | 0.887 | 0.879
1-5 $\left(5,25\%,1.8,1.5\right)$ | - | 0.459 | - | 0.714
Table 7: $R^{2}$ scores of large and small spheres in various datasets with a
universal PINN model.

Given the excellent performance of our PINN modelling on multiple datasets and
its demonstrated interpolation capability to the unseen data, we are now in
the position to establish a universal model that incorporates both variable
fluid and particle properties. With the same logic of avoiding large
increments in parameter space, we omit the datasets with ${\mathcal{R}e}=100$
and combine the rest with ${\mathcal{R}e}=1$ and ${\mathcal{R}e}=10$
comprising over $110,000$ data points. We plot all the PR-DNS data with their
corresponding prediction values obtained by our universal PINN model in Figure
7(a), and list the $R^{2}$ scores for all the six training datasets in Table 7
(please find each case and its corresponding clustered scatters according to
the colors illustrated in the figure and table). Throughout the comparison of
$F_{x}$ predictions generated by the separate and universal models as shown in
Table 3 and 7 respectively, we can observe that the universal model
accomplishes great performance which is highly similar to the model
constructed with independent training. For the cases with relatively larger
mean values such as
$\left({\mathcal{R}e},\phi,d_{l}^{*}/d_{s}^{*},V_{l}^{*}/V_{s}^{*}\right)=\left(10,30\%,2,2\right)$
and $\left(10,40\%,2.5,4\right)$, our universal model produces even slightly
better predictions for both the large and small spheres, which is not
surprising given that (i) more high-quality data are provided into the PINN
model, and (ii) $\delta$ for evaluating the Huber loss is set as the largest
standard deviation among these six datasets, i.e., that of large spheres in
case of $\left(10,40\%,2.5,4\right)$. The success of such a universal model
justifies our previous implementation of training and testing on each
individual dataset, which significantly accelerates the fine-tuning process of
hyper-parameters and offers a fast way to determine the model performance by
merely examining the lower and upper bounds of parameter space.

We eventually generate a dataset with
$\left({\mathcal{R}e},\phi,d_{l}^{*}/d_{s}^{*},V_{l}^{*}/V_{s}^{*}\right)=\left(5,25\%,1.8,1.5\right)$
to examine the interpolation capability of our universal model in the entire
input space. Figure 7(b) and Table 7 show the distribution and $R^{2}$ scores
of predictions on this group of unseen data respectively, and it can be noted
that the predictive performance, especially in the low ${\mathcal{R}e}$
regime, is acceptable qualitatively and quantitatively considering the
complexity of this test, while we can reasonably expect much more accurate
predictions that may be provided by a universal model trained on uniformly
distributed data with higher density. Preliminary findings reported in the
above reveal that interpolating in the $\phi$ dimension of the parameter space
constitutes a major challenge and therefore high density of the training
dataset in the $\phi$ direction should be favored. This is not a surprising
recommendation from a Fluid Mechanics viewpoint since the level of
force/torque fluctuations depend on the strength of the flow disturbances
caused by neighbors, and the magnitude of these flow disturbances primarily
depends on $\phi$.

###  4.4 Interpretability of the PINN model

Apart from the generalization, interpretability is also a critical
consideration in a NN model. It is of great importance for the model to learn
the underlying mechanics and reproduce results that are physically meaningful.
We borrow the idea from Seyed-Ahmadi and Wachs [38] and generate the
contribution from the nearest and second nearest neighbors to the total
force/torque exerted on a target sphere. In fact, we can easily access the
individual contribution of each neighbor giving to that our NN architecture
separately evaluates every pairwise interaction and only sums over all
influential neighbors at the output layer with a set of weights. We first
present the contribution of $F_{x}$, $F_{L}$ and $T_{\perp}$ from the nearest
neighbor in Figures 8$-$10. Please note that in these plots, we put the target
sphere (solid black circle) at the center and the outer white space defines
the region where the neighbor cannot reside. We identify four different
interaction modes: $s-s$, $s-l$, $l-s$ and $l-l$ interactions in bidisperse
particle-laden flows with the first letter representing the size of the target
sphere and the second letter that of its neighbor ($l$: large or $s$: small).
We adopt the same notations as in [38] where $r_{x}$ stands for the streamwise
distance between two spheres and $r_{L}$ for the transverse distance. The
plots generated as follows are processed with an azimuthal average operation
over several planes containing the streamwise $x$ direction (which is also the
reason why we use $r_{L}$ instead of $r_{y}$ or $r_{z}$). The perfect symmetry
about the $x$-axis in each plot is an outcome of data augmentation (we mirror
the dataset about $x$-axis on the $xy$ and $xz$ planes due to the flow
symmetry) but not of this average operation. For the illustration purpose, we
select three cases
$\left({\mathcal{R}e},\phi,d_{l}^{*}/d_{s}^{*},V_{l}^{*}/V_{s}^{*}\right)=\left(1,15\%,2,2\right)$,
$\left(10,15\%,2,2\right)$ and $\left(100,15\%,2,2\right)$ while similar
patterns can be observed in all other cases.

(a)
$\left(Re,\phi,d_{l}^{*}/d_{s}^{*},V_{l}^{*}/V_{s}^{*}\right)=(1,15\%,2,2)$

(b)
$\left(Re,\phi,d_{l}^{*}/d_{s}^{*},V_{l}^{*}/V_{s}^{*}\right)=(10,15\%,2,2)$

(c)
$\left(Re,\phi,d_{l}^{*}/d_{s}^{*},V_{l}^{*}/V_{s}^{*}\right)=(100,15\%,2,2)$

Figure 8: $F_{x}$ contribution of the closest neighbor to the reference
particle in different interaction modes.

(a)
$\left(Re,\phi,d_{l}^{*}/d_{s}^{*},V_{l}^{*}/V_{s}^{*}\right)=(1,15\%,2,2)$

(b)
$\left(Re,\phi,d_{l}^{*}/d_{s}^{*},V_{l}^{*}/V_{s}^{*}\right)=(10,15\%,2,2)$

(c)
$\left(Re,\phi,d_{l}^{*}/d_{s}^{*},V_{l}^{*}/V_{s}^{*}\right)=(100,15\%,2,2)$

Figure 9: $F_{L}$ contribution of the closest neighbor to the reference
particle in different interaction modes.

(a)
$\left(Re,\phi,d_{l}^{*}/d_{s}^{*},V_{l}^{*}/V_{s}^{*}\right)=(1,15\%,2,2)$

(b)
$\left(Re,\phi,d_{l}^{*}/d_{s}^{*},V_{l}^{*}/V_{s}^{*}\right)=(10,15\%,2,2)$

(c)
$\left(Re,\phi,d_{l}^{*}/d_{s}^{*},V_{l}^{*}/V_{s}^{*}\right)=(100,15\%,2,2)$

Figure 10: $T_{\perp}$ contribution of the closest neighbor to the reference
particle in different interaction modes.

These contribution plots are not only similar to those presented in [37, 38,
39], but also reminiscent of the probability distribution maps in [14, 43]. As
an example, the $F_{x}$ contributions in case $\left(100,15\%,2,2\right)$
illustrated in Figure 8 are highly comparable to the combination of
$\tilde{p}\left(\Delta F_{x}>0\right)$ and $\tilde{p}\left(\Delta
F_{x}<0\right)$ in [43], where the high probability regions of a target sphere
experiencing lower than the average streamwise force are constantly located in
mostly upstream and somewhat downstream areas, meanwhile it is shown from the
contribution plots that a neighboring sphere induces a negative $\Delta F_{x}$
in the qualitatively identical regions (in blue). Similarly, the positive
$\Delta F_{x}$ contribution areas (in red), i.e., the high probability regions
of $\Delta F_{x}>0$ since the nearest neighbor plays a dominant role, are
constantly located on both sides of a target sphere. It should be noted that
the red regions are slightly drifted downstream-wise at both
${\mathcal{R}e}=10$ and $100$, while a fore-aft symmetry is captured at
${\mathcal{R}e}=1$ where the flow inertia has a minor effect and the main
source of force fluctuations comes from the neighbors, i.e., the
microstructure. Furthermore, the contributions from large spheres are
consistently larger than those from small spheres, which is also valid in the
comparison between $l-s$ and $s-l$ interactions. Both of these two
observations can be applied to $F_{L}$ and $T_{\perp}$ as shown in Figure 9
and 10 respectively, where the fore-aft symmetry only exists at
${\mathcal{R}e}=1$ and then the critical regions drift correspondingly toward
the $1-4$, and $2-3$ quadrants when ${\mathcal{R}e}$ increases. This trend
coincides with the probability distribution maps as well. In addition, the
upstream neighbors exhibit more effect on $F_{x}$ and less effect on $F_{L}$
at higher ${\mathcal{R}e}$, and vice versa, whereas the influential ranges of
$T_{L}$ are considerably smaller than those of forces, thus fewer neighbors
are required in torque predictions.

Finally, we carry out a similar analysis for the second nearest neighbor and
display its contribution to the total force/torque exerted on a target sphere
(solid black circle) by fixing the nearest neighbor (hollow black circle) at
an upstream position with a gap $=0.2$ (only for demonstration purpose).
Again, the outer white spaces enclosed by the black dashed circles define
regions where the neighbor cannot reside, and the outermost green dashed
circles reflect the range of influential neighbors, i.e., the position of the
50th and 20th neighbor in force and torque predictions respectively. We can
identify eight distinct modes for trinary interactions in bidisperse particle-
laden flows, and denote each mode with a format of target sphere $-$ nearest
and second nearest neighbors. We select the case
$\left({\mathcal{R}e},\phi,d_{l}^{*}/d_{s}^{*},V_{l}^{*}/V_{s}^{*}\right)=\left(100,15\%,2,2\right)$
as an example here, and similar observations can be made in any other case:
(i) Contributions of the second nearest neighbor are constantly smaller than
those of the nearest neighbor as shown by the legend of each subplot, which is
the result of prominently larger weight given to the nearest neighbor to
address its dominant role. (ii) Similar to the previous observation in binary
interactions, large spheres always produce a greater influence regardless of
the interaction mode, implying that our PINN model places particular emphasis
on the size of neighboring spheres throughout the learning process. (iii) The
contour shapes of all three contributions $F_{x}$ $F_{L}$ and $T_{\perp}$ of
the second nearest neighbor are identical to those of the nearest neighbor and
independent of the interaction mode. This is not surprising as in the design
of our PINN architecture, each input neighbor is separately sent to and
processed by the same NN block, thus generating the same influential pattern
with only different magnitude. This unified function/contribution
representation is the core feature of the pairwise interaction assumption. It
simplifies the multi-body interaction problem and enables to separate the
influence of a specific neighbor from that of the other neighbors and to
evaluate it based on the properties of each pair of spheres, whereas the order
of spheres in the local assembly only affects the weights of superposition but
not the general pattern.

(a) $F_{x}$

(b) $F_{L}$

(c) $T_{\perp}$

Figure 11: Contribution of the second closest neighbor to the reference
particle with its closest neighbor fixed at upstream position in different
interaction modes.

##  5 Conclusions

Physics-informed ML and data-driven approaches have been widely adopted to
tackle fluid mechanics problems and several works have demonstrated their
promising application to predict the hydrodynamic force/torque fluctuations in
monodisperse particle-laden flow problems. These fluctuations were usually
neglected in the past although PR-DNS data validated that the standard
deviation of drag distribution has a comparable magnitude as its mean value
[62, 14]. However, we still lack the knowledge of force/torque distributions
in polydisperse particle-laden flows due to the increased computing cost of
the corresponding PR-DNS, thus compromising further research including the
extension of ML techniques to this field. Consequently, we performed a series
of PR-DNS of the flow past a random array of stationary bidisperse and
polydisperse spheres in our previous work [43]. Results exhibited that the
force and torque exerted on each sphere class still follow a Gaussian
distribution with a non-negligible span, and the probability distribution maps
of different binary interaction modes in bidisperse particle-laden flows are
strikingly similar in terms of contour shapes and critical regions. On top of
these numerical findings, we are now able to extract the key characteristics
from this large dataset and leverage a PINN model to exploit the statistical
relationship between the force/torque fluctuations and relevant variables in
bidisperse particle-laden flows as a starting point.

Inspired by several example works on monodisperse particle-laden flows [37,
38, 39], we formulated our problem of interest based on the same pairwise
interaction assumption [46, 47] that approximates the hydrodynamic
force/torque of each sphere in a stationary random array as a linear
superposition of individual contributions from $M$ influential closest
neighbors. As another core assumption, the similar probability distribution
maps of different interaction modes imply a unified function to describe such
a force/torque contribution from each neighbor regardless of its size or
distance to its target sphere. Therefore, we construct a compact PINN
architecture that takes relative position and diameter
$\left\\{\bm{r}_{j},d_{j}\right\\}$ of neighbor $j$ as inputs to evaluate its
effect on a target sphere through the same NN block, which is a reflection of
unified function representation similar to that proposed in [38]. To
incorporate the effects of local mesoscale flow and microstructure, we select
the local average velocity, local solid volume fraction and target sphere
diameter $\left\\{\bm{\mathrm{{\mathcal{R}e}}}_{i},\phi_{i},d_{i}\right\\}$ as
additional inputs and process them with another NN block. The scalar outputs
of two NN blocks experience an element-wise multiplication to manifest the
contribution of each neighbor with the consideration of local flow conditions.
At the final stage, we compute the total force/torque exerted on a target
sphere using weighted sums according to the pairwise interaction assumption.

We determined the number of influential neighbors as $M=50$ and $20$ in force
and torque predictions respectively, and correspondingly the depth (number of
layers) and width of both NN blocks being $\left(3,30\right)$ and
$\left(3,20\right)$ according to the grid-search results. We trained our PINN
model with PR-DNS datasets of eight PR-DNS cases spanning
$1\leq{\mathcal{R}e}\leq 100$, $10\%\leq\phi\leq 40\%$, $1.5\leq
d_{l}^{*}/d_{s}^{*}\leq 2.5$ and $1\leq V_{l}^{*}/V_{s}^{*}\leq 4$. Through a
comparison of PINN model predictions and PR-DNS data, we demonstrated great
overall performance with $R^{2}$ in a range of $0.706\sim 0.892$ for $F_{x}$,
$0.669\sim 0.871$ for $F_{L}$ and $0.741\sim 0.910$ for $T_{\perp}$
predictions. The success of the current model does not only reflect in the
close performance of small and larger spheres with distinct force/torque
distributions, but also in the minor difference in $R^{2}$ of training and
test sets which indicates excellent model generalization to the data that it
has not trained on. We further compared the force/torque distribution
statistics of PINN predictions to those of PR-DNS data and obtained very
similar distributions with both approaches, including a good agreement of
their respective standard deviations and a remarkably small error of $0.18\%$
on the mean $F_{x}$ even in the worst case scenario. Subsequently, we repeated
the training procedures with the same hyper-parameters on multiple combined
datasets to construct a universal predictive model, and evaluated the model
performance on the unseen data again generated by PR-DNS with new sets of
input parameters. Satisfactory predictions were reported in most cases,
indicating good interpolation capability of our universal PINN model within
the considered range of fluid and particle properties. Finally, we showed the
contribution maps of the nearest and second nearest neighbors in different
binary and trinary interaction modes. We noticed that the critical regions,
i.e., positive or negative force/torque contributions, are qualitatively
identical to those shown in our previous probability distribution maps,
implying great interpretability provided by our PINN predictions.

As we mentioned previously, the bidisperse particle-laden flow cases examined
here are the starting point of our "proof-of-concept" study. Our ultimate goal
is to evolve from a problem-specific model to a versatile model which is
capable of coping with other particle-laden flow problems involving more
complicated configurations, e.g., flows laden with polydisperse spheres and/or
non-spherical particles. This cannot be fulfilled without the prior knowledge
we supplied in the present work and the fine-tuned design of the corresponding
PINN architecture. As pointed out by Seyed-Ahmadi and Wachs [38], the simple
fully-connected NN model even fails to deliver correct predictions with an
increasing number of influential neighbors. On the contrary, the PINN
architecture used in the present work possesses several key advantages that
will potentially facilitate its extension to more complex particle-laden flow
problems:

  1. 1.

Due to the common NN block shared by different neighbors, we are theoretically
able to incorporate an infinite number of influential neighbors as inputs
without adjusting the overall PINN architecture, whereas the size (both depth
and width) of hidden layers must be accordingly expanded in this case for a
fully-connected NN model. This is crucial for strongly bidisperse and
polydisperse particle-laden flows as the number of influential neighbors
should be increased when we consider a large sphere diameter ratio. (as an
example, $M=50$ and 20 in our bidisperse cases with only $1.5\leq
d_{l}^{*}/d_{s}^{*}\leq 2.5$, while $M\leq 30$ and 10 in monodisperse cases in
previous works).

  2. 2.

The unified function representation [38] (similar to the order invariance
assumption in [37]) has been proven to be suitable in bidisperse cases based
on our observations of probability distribution maps obtained with different
binary interaction modes, hence inspiring this particular architecture of
shared NN block for every neighbor-related input sets. Therefore, we assume
that this feature could be applied to any flow laden with polydisperse
spherical particles with the same inputs $\left\\{\bm{r}_{j},d_{j}\right\\}$
of neighbor $j$ (still needs to be trained and validated with PR-DNS data),
and similarly to flows laden with non-spherical particles with additional
independent variables to characterize neighbors, e.g., sphericity and aspect
ratio, to be inserted in the approximating contribution functions and
incorporated into the input space.

  3. 3.

We separately evaluated the effect of local mesoscale flow and neighborhood
based on $\left\\{\bm{\mathrm{{\mathcal{R}e}}}_{i},\phi_{i},d_{i}\right\\}$
with another NN block, and combined this effect with each neighbor
contribution throughout an element-wise multiplication as adopted in [39].
Such an operation is more physically reasonable than directly supplying these
inputs to the linear superposition which seems oversimplified in [38]. Thanks
to this second NN block, we are also allowed to consider many other
intermediate variables, e.g., pressure components as suggested in [35], to
improve the force/torque predictions.

It is worth highlighting another bonus promoted by such a compact PINN
architecture that is also pointed out in the previous section: the total
number of unknown parameters is tremendously decreased compared to a fully
connected NN. For example, the number of unknown parameters is 4100/4150 in
force predictions ($F_{L}$ has two components) and 1940 in torque predictions,
whereas it is stunningly 25850 for a (6,50) NN and 81500 for a (5,100) NN as
reported in [39] (the reason is that such a fully connected NN is assigned to
correlate a much larger number of input variables, i.e., inputs of all
neighbors versus inputs of a single neighbor in our PINN model). Consequently,
our PINN model is much more efficient not only in terms of training (i.e.,
optimizing the value of unknown parameters), but more importantly in terms of
evaluating the force/torque distributions for a random sphere assembly. In
addition, the elu activation function used in our model has a simpler form
than the hyperbolic tangent unit that also contributes to faster computation.
As an example, it takes $<1$ second for our model to evaluate the force or
torque component in an array of $4000$ spheres with 16 CPU-core parallel
computing. Therefore, it is promising to implement our PINN model in EL
simulations to provide the missing hydrodynamic particle force and torque
exerted by the locally complex fluid flow, which are usually oversimplified by
any average drag closure (as the other components have a mean value of $0$).
Undoubtedly, our PINN model is still slower than the average drag closure
computation since these closures typically have a concise functional form
dependent on mesoscale variables only, such as ${\mathcal{R}e}$, $\phi$ and
$d_{i}^{*}$. Eventually, there is no free lunch but there always exists an
appropriate trade-off between accuracy and computing cost, depending on the
level of fidelity expected from a simulation. We believe that our PINN model
that provides good performance in both aspects can arguably become a strong
alternative to enhance the EL simulations.

We would like to share final remarks about the fundamental assumption of
linear superposition for force/torque exerted on a targe sphere. In our work,
we consistently apply this linear superposition and achieve great predictions,
most probably because the upper bound of the selected range
$1\leq{\mathcal{R}e}\leq\ 100$ is relatively low, whereas we can reasonably
anticipate that a non-linear relation should be employed at higher
${\mathcal{R}e}$. Siddani and Balachandar [39] suggest a hierarchical
interpretation of hydrodynamic force/torque and expand it into unary, binary,
trinary and higher-order interaction terms. This is undoubtedly adequate but
training a NN model to learn trinary and higher order interactions is far more
challenging. However, if we are able to find the relation between binary and
trinary terms, and successively between $N$th and $\left(N+1\right)$th terms,
a non-linear relation between force/torque and the binary interaction terms
only might be equivalently deduced. Correspondingly in the PINN architecture,
we merely need to add an extra NN block to the output layers to approximate
this non-linear function, where each binary interaction contribution is taken
as input and force/torque as output. We leave this to our future research.

## References

## References

  * [1] M. A. van der Hoef, M. van Sint Annaland, N. Deen, J. A. M. Kuipers, Numerical simulation of dense gas-solid fluidized beds: A multiscale modeling strategy, Annual Review of Fluid Mechanics 40 (1) (2008) 47–70.  doi:10.1146/annurev.fluid.40.111406.102130.
URL https://doi.org/10.1146/annurev.fluid.40.111406.102130

  * [2] S. Balachandar, J. K. Eaton, Turbulent dispersed multiphase flow, Annu. Rev. Fluid Mech. 42 (1) (2009) 111–133.  doi:10.1146/annurev.fluid.010908.165243.
URL https://doi.org/10.1146/annurev.fluid.010908.165243

  * [3] S. Tenneti, S. Subramaniam, Particle-resolved direct numerical simulation for gas-solid flow model development, Annual Review of Fluid Mechanics 46 (1) (2014) 199–230.  arXiv:https://doi.org/10.1146/annurev-fluid-010313-141344, doi:10.1146/annurev-fluid-010313-141344.
URL https://doi.org/10.1146/annurev-fluid-010313-141344

  * [4] L. Brandt, F. Coletti, Particle-laden turbulence: Progress and perspectives, Annu. Rev. Fluid Mech. 54 (1) (2022) 159–189.  doi:10.1146/annurev-fluid-030121-021103.
URL https://doi.org/10.1146/annurev-fluid-030121-021103

  * [5] X. Wang, K. Liu, C. You, Drag force model corrections based on nonuniform particle distributions in multi-particle systems, Powder Technology 209 (1) (2011) 112–118.  doi:10.1016/j.powtec.2011.02.018.
URL https://www.sciencedirect.com/science/article/pii/S0032591011000738

  * [6] G. J. Rubinstein, A. Ozel, X. Yin, J. J. Derksen, S. Sundaresan, Lattice Boltzmann simulations of low-Reynolds-number flows past fluidized spheres: effect of inhomogeneities on the drag force, Journal of Fluid Mechanics 833 (2017) 599–630.  doi:10.1017/jfm.2017.705.
  * [7] M. Su, H. Zhao, Modifying the inter-phase drag via solid volume fraction gradient for CFD simulation of fast fluidized beds, AIChE J. 63 (7) (2017) 2588–2598.  doi:10.1002/aic.15573.
URL https://doi.org/10.1002/aic.15573

  * [8] T. Li, L. Wang, W. Rogers, G. Zhou, W. Ge, An approach for drag correction based on the local heterogeneity for gas-solid flows, AIChE J. 63 (4) (2017) 1203–1212.  doi:10.1002/aic.15507.
URL https://doi.org/10.1002/aic.15507

  * [9] S. Tenneti, M. Mehrabadi, S. Subramaniam, Stochastic Lagrangian model for hydrodynamic acceleration of inertial particles in gas-solid suspensions, Journal of Fluid Mechanics 788 (2016) 695–729.  doi:10.1017/jfm.2015.693.
  * [10] A. Esteghamatian, F. Euzenat, A. Hammouti, M. Lance, A. Wachs, A stochastic formulation for the drag force based on multiscale numerical simulation of fluidized beds, International Journal of Multiphase Flow 99 (2018) 363–382.  doi:https://doi.org/10.1016/j.ijmultiphaseflow.2017.11.003.
URL https://www.sciencedirect.com/science/article/pii/S0301932217304779

  * [11] A. M. Lattanzi, V. Tavanashad, S. Subramaniam, J. Capecelatro, Stochastic models for capturing dispersion in particle-laden flows, Journal of Fluid Mechanics 903 (A7) (2020) A7.  doi:10.1017/jfm.2020.625.
  * [12] A. M. Lattanzi, V. Tavanashad, S. Subramaniam, J. Capecelatro, Stochastic model for the hydrodynamic force in Euler-Lagrange simulations of particle-laden flows, Physical Review Fluids 7 (1) (2022) 014301.  doi:10.1103/PhysRevFluids.7.014301.
URL https://link.aps.org/doi/10.1103/PhysRevFluids.7.014301

  * [13] B. Hardy, O. Simonin, J. De Wilde, G. Winckelmans, Simulation of the flow past random arrays of spherical particles: Microstructure-based tensor quantities as a tool to predict fluid-particle forces, International Journal of Multiphase Flow 149 (2022) 103970.  doi:https://doi.org/10.1016/j.ijmultiphaseflow.2021.103970.
URL https://www.sciencedirect.com/science/article/pii/S0301932221003621

  * [14] A. Seyed-Ahmadi, A. Wachs, Microstructure-informed probability-driven point-particle model for hydrodynamic forces and torques in particle-laden flows, Journal of Fluid Mechanics 900 (2020) A21.  doi:10.1017/jfm.2020.453.
  * [15] A. N. Kolmogorov, The local structure of turbulence in incompressible viscous fluid for very large Reynolds number, Dokl. Akad. Nauk. SSSR 30 (1941) 301–303.
URL https://cir.nii.ac.jp/crid/1572824499981613312

  * [16] S. L. Brunton, B. R. Noack, P. Koumoutsakos, Machine learning for fluid mechanics, Annu. Rev. Fluid Mech. 52 (1) (2020) 477–508.  doi:10.1146/annurev-fluid-010719-060214.
URL https://doi.org/10.1146/annurev-fluid-010719-060214

  * [17] M. Raissi, P. Perdikaris, G. E. Karniadakis, Machine learning of linear differential equations using Gaussian processes, Journal of Computational Physics 348 (2017) 683–693.  doi:10.1016/j.jcp.2017.07.050.
URL https://www.sciencedirect.com/science/article/pii/S0021999117305582

  * [18] M. Raissi, G. E. Karniadakis, Hidden physics models: Machine learning of nonlinear partial differential equations, Journal of Computational Physics 357 (2018) 125–141.  doi:10.1016/j.jcp.2017.11.039.
URL https://www.sciencedirect.com/science/article/pii/S0021999117309014

  * [19] M. Raissi, P. Perdikaris, G. E. Karniadakis, Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations, Journal of Computational Physics 378 (2019) 686–707.  doi:10.1016/j.jcp.2018.10.045.
URL https://www.sciencedirect.com/science/article/pii/S0021999118307125

  * [20] M. Raissi, Z. Wang, M. S. Triantafyllou, G. E. Karniadakis, Deep learning of vortex-induced vibrations, Journal of Fluid Mechanics 861 (2019) 119–137.  doi:10.1017/jfm.2018.872.
  * [21] S. Cai, Z. Mao, Z. Wang, M. Yin, G. E. Karniadakis, Physics-informed neural networks (PINNs) for fluid mechanics: a review, Acta Mechanica Sinica 37 (12) (2021) 1727–1738.  doi:10.1007/s10409-021-01148-1.
URL https://doi.org/10.1007/s10409-021-01148-1

  * [22] G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, L. Yang, Physics-informed machine learning, Nature Reviews Physics 3 (6) (2021) 422–440.  doi:10.1038/s42254-021-00314-5.
URL https://doi.org/10.1038/s42254-021-00314-5

  * [23] S. Cuomo, V. S. Di Cola, F. Giampaolo, G. Rozza, M. Raissi, F. Piccialli, Scientific machine learning through physics-informed neural networks: where we are and what’s next, Journal of Scientific Computing 92 (3) (2022) 88.  doi:10.1007/s10915-022-01939-z.
URL https://doi.org/10.1007/s10915-022-01939-z

  * [24] L.-T. Zhu, X.-Z. Chen, B. Ouyang, W.-C. Yan, H. Lei, Z. Chen, Z.-H. Luo, Review of machine learning for hydrodynamics, transport, and reactions in multiphase flows and reactors, Ind. Eng. Chem. Res. 61 (28) (2022) 9901–9949.  doi:10.1021/acs.iecr.2c01036.
URL https://doi.org/10.1021/acs.iecr.2c01036

  * [25] A. D. Jagtap, E. Kharazmi, G. E. Karniadakis, Conservative physics-informed neural networks on discrete domains for conservation laws: Applications to forward and inverse problems, Computer Methods in Applied Mechanics and Engineering 365 (2020) 113028.  doi:10.1016/j.cma.2020.113028.
URL https://www.sciencedirect.com/science/article/pii/S0045782520302127

  * [26] E. Kharazmi, Z. Zhang, G. E. M. Karniadakis, hp-vpinns: Variational physics-informed neural networks with domain decomposition, Computer Methods in Applied Mechanics and Engineering 374 (2021) 113547.  doi:10.1016/j.cma.2020.113547.
URL https://www.sciencedirect.com/science/article/pii/S0045782520307325

  * [27] P.-H. Chiu, J. C. Wong, C. Ooi, M. H. Dao, Y.-S. Ong, Can-pinn: A fast physics-informed neural network based on coupled-automatic-numerical differentiation method, Computer Methods in Applied Mechanics and Engineering 395 (2022) 114909.  doi:10.1016/j.cma.2022.114909.
URL https://www.sciencedirect.com/science/article/pii/S0045782522001906

  * [28] X. Jin, S. Cai, H. Li, G. E. Karniadakis, NSFnets (Navier-Stokes flow nets): Physics-informed neural networks for the incompressible Navier-Stokes equations, Journal of Computational Physics 426 (2021) 109951.  doi:10.1016/j.jcp.2020.109951.
URL https://www.sciencedirect.com/science/article/pii/S0021999120307257

  * [29] S. A. Faroughi, N. Pawar, C. Fernandes, S. Das, N. K. Kalantari, S. K. Mahjour, Physics-guided, physics-informed, and physics-encoded neural networks in scientific computing (2022).  doi:10.48550/ARXIV.2211.07377.
URL https://arxiv.org/abs/2211.07377

  * [30] J. Viquerat, E. Hachem, A supervised neural network for drag prediction of arbitrary 2d shapes in laminar flows at low Reynolds number, Computers & Fluids 210 (2020) 104645\.  doi:10.1016/j.compfluid.2020.104645.
URL https://www.sciencedirect.com/science/article/pii/S0045793020302164

  * [31] D. Xin, J. Zeng, K. Xue, Surrogate drag model of non-spherical fragments based on artificial neural networks, Powder Technology 404 (2022) 117412.  doi:10.1016/j.powtec.2022.117412.
URL https://www.sciencedirect.com/science/article/pii/S0032591022003060

  * [32] S. Lee, J. Yang, P. Forooghi, A. Stroh, S. Bagheri, Predicting drag on rough surfaces by transfer learning of empirical correlations, Journal of Fluid Mechanics 933 (A18) (2022) A18.  doi:10.1017/jfm.2021.1041.
  * [33] S. A. Faroughi, A. I. Roriz, C. Fernandes, A Meta-model to predict the drag coefficient of a particle translating in viscoelastic fluids: A machine learning approach (2022).  doi:10.3390/polym14030430.
  * [34] L. He, D. K. Tafti, A supervised machine learning approach for predicting variable drag forces on spherical particles in suspension, Powder Technology 345 (2019) 379–389.  doi:10.1016/j.powtec.2019.01.013.
URL https://www.sciencedirect.com/science/article/pii/S0032591019300130

  * [35] N. Muralidhar, J. Bu, Z. Cao, L. He, N. Ramakrishnan, D. Tafti, A. Karpatne, Physics-guided design and learning of neural networks for predicting drag force on particle suspensions in moving fluids (2019).  doi:10.48550/ARXIV.1911.04240.
URL https://arxiv.org/abs/1911.04240

  * [36] N. Muralidhar, J. Bu, Z. Cao, L. He, N. Ramakrishnan, D. Tafti, A. Karpatne, Physics-guided deep learning for drag force prediction in dense fluid-particulate systems, Big Data 8 (5) (2020) 431–449.  doi:10.1089/big.2020.0071.
URL https://doi.org/10.1089/big.2020.0071

  * [37] W. Moore, S. Balachandar, G. Akiki, A hybrid point-particle force model that combines physical and data-driven approaches, Journal of Computational Physics 385 (2019) 187–208.  doi:https://doi.org/10.1016/j.jcp.2019.01.053.
URL https://www.sciencedirect.com/science/article/pii/S0021999119301111

  * [38] A. Seyed-Ahmadi, A. Wachs, Physics-inspired architecture for neural network modeling of forces and torques in particle-laden flows, Computers & Fluids 238 (2022) 105379.  doi:https://doi.org/10.1016/j.compfluid.2022.105379.
URL https://www.sciencedirect.com/science/article/pii/S004579302200055X

  * [39] B. Siddani, S. Balachandar, Point-particle drag, lift, and torque closure models using machine learning: Hierarchical approach and interpretability, Physical Review Fluids 8 (1) (2023) 014303.  doi:10.1103/PhysRevFluids.8.014303.
URL https://link.aps.org/doi/10.1103/PhysRevFluids.8.014303

  * [40] M. Shirzadi, T. Fukasawa, K. Fukui, T. Ishigami, Application of deep learning neural networks for the analysis of fluid-particle dynamics in fibrous filters, Chemical Engineering Journal 455 (2023) 140775.  doi:10.1016/j.cej.2022.140775.
URL https://www.sciencedirect.com/science/article/pii/S1385894722062556

  * [41] N. R. Ashwin, Z. Cao, N. Muralidhar, D. Tafti, A. Karpatne, Deep learning methods for predicting fluid forces in dense particle suspensions, Powder Technology 401 (2022) 117303.  doi:10.1016/j.powtec.2022.117303.
URL https://www.sciencedirect.com/science/article/pii/S0032591022001978

  * [42] S. Hwang, J. Pan, L.-S. Fan, Deep learning for drag force modelling in dilute, poly-dispersed particle-laden flows with irregular-shaped particles, Chemical Engineering Science 266 (2023) 118299.  doi:10.1016/j.ces.2022.118299.
URL https://www.sciencedirect.com/science/article/pii/S0009250922008843

  * [43] Z. Cheng, A. Wachs, Hydrodynamic force and torque fluctuations in a random array of polydisperse stationary spheres, International Journal of Multiphase Flow.
  * [44] Z. Cheng, A. Wachs, An immersed boundary/multi-relaxation time lattice Boltzmann method on adaptive octree grids for the particle-resolved simulation of particle-laden flows, Journal of Computational Physics 471 (2022) 111669.  doi:https://doi.org/10.1016/j.jcp.2022.111669.
URL https://www.sciencedirect.com/science/article/pii/S002199912200732X

  * [45] S. Popinet, A quadtree-adaptive multigrid solver for the Serre-Green-Naghdi equations, Journal of Computational Physics 302 (2015) 336–358.  doi:10.1016/j.jcp.2015.09.009.
URL https://www.sciencedirect.com/science/article/pii/S0021999115005902

  * [46] G. Akiki, W. Moore, S. Balachandar, Pairwise-interaction extended point-particle model for particle-laden flows, Journal of Computational Physics 351 (2017) 329–357.  doi:https://doi.org/10.1016/j.jcp.2017.07.056.
URL https://www.sciencedirect.com/science/article/pii/S0021999117306848

  * [47] G. Akiki, T. L. Jackson, S. Balachandar, Pairwise interaction extended point-particle model for a random array of monodisperseÂ spheres, Journal of Fluid Mechanics 813 (2017) 882–928.  doi:10.1017/jfm.2016.877.
  * [48] S. Balachandar, W. C. Moore, G. Akiki, K. Liu, Toward particle-resolved accuracy in Euler-Lagrange simulations of multiphase flow using machine learning and pairwise interaction extended point-particle (PIEP) approximation, Theoretical and Computational Fluid Dynamics 34 (4) (2020) 401–428.  doi:10.1007/s00162-020-00538-8.
URL https://doi.org/10.1007/s00162-020-00538-8

  * [49] D.-A. Clevert, T. Unterthiner, S. Hochreiter, Fast and accurate deep network learning by exponential linear units (elus) (2015).  doi:10.48550/ARXIV.1511.07289.
URL https://arxiv.org/abs/1511.07289

  * [50] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfellow, A. Harp, G. Irving, M. Isard, Y. Jia, R. Jozefowicz, L. Kaiser, M. Kudlur, J. Levenberg, D. Mane, R. Monga, S. Moore, D. Murray, C. Olah, M. Schuster, J. Shlens, B. Steiner, I. Sutskever, K. Talwar, P. Tucker, V. Vanhoucke, V. Vasudevan, F. Viegas, O. Vinyals, P. Warden, M. Wattenberg, M. Wicke, Y. Yu, X. Zheng, Tensorflow: Large-Scale Machine Learning on Heterogeneous Distributed Systems (2016).  doi:10.48550/ARXIV.1603.04467.
URL https://arxiv.org/abs/1603.04467

  * [51] L. Wood, Z. Tan, S. Ian, S. Zhu, F. Chollet, et al., Kerascv, https://github.com/keras-team/keras-cv (2022).
  * [52] T. Dozat, Incorporating Nesterov Momentum into Adam, in: Proceedings of the 4th International Conference on Learning Representations, 2016, pp. 1–4.
  * [53] M. A. van der Hoef, R. Beetstra, J. A. M. Kuipers, Lattice-Boltzmann simulations of low-Reynolds-number flow past mono- and bidisperse arrays of spheres: results for the permeability and drag force, Journal of Fluid Mechanics 528 (2005) 233–254.  doi:10.1017/S0022112004003295.
  * [54] R. Beetstra, M. A. van der Hoef, J. A. M. Kuipers, Drag force of intermediate Reynolds number flow past mono- and bidisperse arrays of spheres, AIChE Journal 53 (2) (2007) 489–501.  doi:https://doi.org/10.1002/aic.11065.
URL https://aiche.onlinelibrary.wiley.com/doi/abs/10.1002/aic.11065

  * [55] X. Yin, S. Sundaresan, Fluid-particle drag in low-Reynolds-number polydisperse gas-solid suspensions, AIChE Journal 55 (6) (2009) 1352–1368.  doi:https://doi.org/10.1002/aic.11800.
URL https://aiche.onlinelibrary.wiley.com/doi/abs/10.1002/aic.11800

  * [56] S. Sarkar, M. van der Hoef, J. Kuipers, Fluid-particle interaction from lattice Boltzmann simulations for flow through polydisperse random arrays of spheres, Chemical Engineering Science 64 (11) (2009) 2683–2691.  doi:https://doi.org/10.1016/j.ces.2009.02.045.
URL https://www.sciencedirect.com/science/article/pii/S0009250909001511

  * [57] W. Holloway, X. Yin, S. Sundaresan, Fluid-particle drag in inertial polydisperse gas-solid suspensions, AIChE Journal 56 (8) (2010) 1995–2004.  doi:https://doi.org/10.1002/aic.12127.
URL https://aiche.onlinelibrary.wiley.com/doi/abs/10.1002/aic.12127

  * [58] F. Cello, A. Di Renzo, F. P. Di Maio, A semi-empirical model for the drag force and fluid-particle interaction in polydisperse suspensions, Chemical Engineering Science 65 (10) (2010) 3128–3139.  doi:https://doi.org/10.1016/j.ces.2010.02.006.
URL https://www.sciencedirect.com/science/article/pii/S0009250910000722

  * [59] L. Rong, K. Dong, A. Yu, Lattice-Boltzmann simulation of fluid flow through packed beds of spheres: Effect of particle size distribution, Chemical Engineering Science 116 (2014) 508–523.  doi:https://doi.org/10.1016/j.ces.2014.05.025.
URL https://www.sciencedirect.com/science/article/pii/S0009250914002358

  * [60] M. Mehrabadi, S. Tenneti, S. Subramaniam, Importance of the fluid-particle drag model in predicting segregation in bidisperse gas-solid flow, International Journal of Multiphase Flow 86 (2016) 99–114.  doi:https://doi.org/10.1016/j.ijmultiphaseflow.2016.07.006.
URL https://www.sciencedirect.com/science/article/pii/S0301932216300738

  * [61] F. Duan, L. Zhao, X. Chen, Q. Zhou, Fluid-particle drag and particle-particle drag in low-Reynolds-number bidisperse gas-solid suspensions, Physics of Fluids 32 (11) (2020) 113311.  doi:10.1063/5.0023874.
URL https://doi.org/10.1063/5.0023874

  * [62] G. Akiki, T. L. Jackson, S. Balachandar, Force variation within arrays of monodisperse spherical particles, Phys. Rev. Fluids 1 (2016) 044202.  doi:10.1103/PhysRevFluids.1.044202.
URL https://link.aps.org/doi/10.1103/PhysRevFluids.1.044202

###### List of Figures

  1. 1 Comparison of $F_{x}$ distribution of large and small spheres in each case obtained with PR-DNS.
  2. 2 Physics-informed neural network architecture used in the present work.
  3. 3 $R^{2}$ scores of $F_{x}$, $F_{L}$ and $T_{\perp}$ for the increasing number of neighbors $M$ included in the PINN model. We select $M=50$ and $M=20$ for force and torque predictions respectively based on this convergence study.
  4. 4 Evolution of Huber loss during the model training process. From left to right: $F_{x}$, $F_{L}$ and $T_{\perp}$.
  5. 5 Comparison of $F_{x}$, $F_{L}$ and $T_{\perp}$ obtained with our PR-DNS and PINN predictions.
  6. (a) $\left({\mathcal{R}e},\phi,d_{l}^{*}/d_{s}^{*},V_{l}^{*}/V_{s}^{*}\right)=\left(1,15\%,2,2\right)$
  7. (b) $\left({\mathcal{R}e},\phi,d_{l}^{*}/d_{s}^{*},V_{l}^{*}/V_{s}^{*}\right)=\left(100,30\%,2,2\right)$
  8. 6 Prediction on the unseen data (green dots) using the PINN models trained on integrated datasets (red and blue dots for large and small spheres respectively).
  9. (a) Interpolation of ${\mathcal{R}e}$
  10. (b) Interpolation of $\phi$
  11. (c) Interpolation of all particle properties
  12. 7 Training of a universal model. (a) Comparison of $F_{x}$ obtained with our PR-DNS and universal PINN predictions. Left: large spheres, and right: small spheres. (b) Interpolation of all parameters.
  13. (a) Universal model
  14. (b) Predictions on unseen data
  15. 8 $F_{x}$ contribution of the closest neighbor to the reference particle in different interaction modes.
  16. (a) $\left(Re,\phi,d_{l}^{*}/d_{s}^{*},V_{l}^{*}/V_{s}^{*}\right)=(1,15\%,2,2)$
  17. (b) $\left(Re,\phi,d_{l}^{*}/d_{s}^{*},V_{l}^{*}/V_{s}^{*}\right)=(10,15\%,2,2)$
  18. (c) $\left(Re,\phi,d_{l}^{*}/d_{s}^{*},V_{l}^{*}/V_{s}^{*}\right)=(100,15\%,2,2)$
  19. 9 $F_{L}$ contribution of the closest neighbor to the reference particle in different interaction modes.
  20. (a) $\left(Re,\phi,d_{l}^{*}/d_{s}^{*},V_{l}^{*}/V_{s}^{*}\right)=(1,15\%,2,2)$
  21. (b) $\left(Re,\phi,d_{l}^{*}/d_{s}^{*},V_{l}^{*}/V_{s}^{*}\right)=(10,15\%,2,2)$
  22. (c) $\left(Re,\phi,d_{l}^{*}/d_{s}^{*},V_{l}^{*}/V_{s}^{*}\right)=(100,15\%,2,2)$
  23. 10 $T_{\perp}$ contribution of the closest neighbor to the reference particle in different interaction modes.
  24. (a) $\left(Re,\phi,d_{l}^{*}/d_{s}^{*},V_{l}^{*}/V_{s}^{*}\right)=(1,15\%,2,2)$
  25. (b) $\left(Re,\phi,d_{l}^{*}/d_{s}^{*},V_{l}^{*}/V_{s}^{*}\right)=(10,15\%,2,2)$
  26. (c) $\left(Re,\phi,d_{l}^{*}/d_{s}^{*},V_{l}^{*}/V_{s}^{*}\right)=(100,15\%,2,2)$
  27. 11 Contribution of the second closest neighbor to the reference particle with its closest neighbor fixed at upstream position in different interaction modes.
  28. (a) $F_{x}$
  29. (b) $F_{L}$
  30. (c) $T_{\perp}$