# High-Precision Surrogate Modeling for Uncertainty Quantification in Complex
Slurry Flows

M. Elkarii  MINES ParisTech, PSL Research University, CES, Paris, France
Mohammed VI Polytechnic University (UM6P), College of Computing, Benguerir,
Morocco  R. Boukharfane  Mohammed VI Polytechnic University (UM6P), College of
Computing, Benguerir, Morocco<EMAIL_ADDRESS>N. El Moçayd
Mohammed VI Polytechnic University (UM6P), College of Agriculture and
Environmental Sciences, Benguerir, Morocco

###### Abstract

Slurry transportation via pipelines is essential for global industries,
offering efficiency and environmental benefits. Specifically, the precise
calibration of physical parameters for transporting raw phosphate material to
fertilizer plants is crucial to minimize energy losses and ensure secure
operations. Computational fluid dynamics (CFD) is commonly employed to
understand solid concentration, velocity distributions, and flow pressure
along the pipeline. However, numerical solutions for slurry flows often entail
uncertainties from initial and boundary conditions, emphasizing the need for
quantification. This study addresses the challenge by proposing a framework
that combines proper orthogonal decomposition and polynomial chaos expansions
to quantify uncertainties in two-dimensional phosphate slurry flow
simulations. The use of surrogate modeling methods, like polynomial chaos
expansion, proves effective in reducing computational costs associated with
direct stochastic simulations, especially for complex flows with high spatial
variability, as observed in phosphate slurries. Numerical results demonstrate
the accuracy of the non-intrusive reduction method in reproducing mean and
variance distributions. Moreover, the uncertainty quantification analysis
shows that the reduced-order model significantly reduces computational costs
compared to the full-order model.

##  I Introduction

The use of numerical modeling to simulate flows of slurry in pipelines has
significantly increased, enabling effective design through a thorough
understanding of slurry flow characteristics. Various mixtures have undergone
computational fluid dynamics (CFD) simulations. Dixit  _et al._ (2022)
investigated sand-water slurry flow in a pipe bend, focusing on assessing head
loss. Sontti  _et al._ (2023) studied a four-phase slurry used in oil sand ore
processing, examining the impact of secondary phases on flow characteristics.
Wan  _et al._ (2023) analyzed hydrodynamic characteristics and flow patterns
of large particles in deep-sea hydraulic lifting pipes, considering factors
such as feed concentration, initial velocity, and particle gradation. Zhou
_et al._ (2023) explored transport in confined spaces, implementing a bimodal
size distribution to enhance accuracy. Joshi, Parkash, and Krishan (2023)
investigated the settling patterns of solid particles in a slurry pipeline,
focusing on the influence of Prandtl numbers. These studies highlight the
growing application of numerical modeling in analyzing slurry pipelines,
providing valuable insights into complex flow behavior. Yet, uncertainty
remains ubiquitous in all those simulations. It is worth to recognize the
diverse mechanical perspectives employed in modeling slurry transport. Some
studies adopt an approach wherein the slurry is treated as a cohesive fluid
continuum, often characterized as non-Newtonian; an illustrative example can
be found in El Moçayd and Seaid (2021). Conversely, alternative investigations
conceptualize the flow as solid suspended particles evolving within a carrier
Newtonian fluid. Simulations conducted based on both modeling strategies are
inherently susceptible to uncertainties. However, it is imperative to
emphasize that the present study exclusively concentrates on simulations
executed utilizing the latter strategy. All simulations performed in this
study utilize the well-known open-source solver, OpenFOAM. This choice is
motivated by its demonstrated efficacy in simulating such flows, as evidenced
by our prior research endeavors (see Elkarii  _et al._ (2022)). Furthermore,
the software’s widespread use across diverse contexts, including industrial
applications, supports its selection. In industrial scenarios where data are
limited, practical experiences are unfeasible, and flow configurations are
seldom predictable, OpenFOAM proves to be a valuable tool. The primary
objective of the present work is to assess uncertainties and conduct
sensitivity analyses in such contexts, where decisions must be made prior to
the validation of the simulation. There exist various sources of uncertainty
when conducting numerical simulations. For slurry flow, uncertainties arise
from limited knowledge of initial and boundary conditions, material
properties, and model parameters. Assumptions made in physical approximations
and numerical parameterization also introduce uncertainties. Current research
focuses on evaluating uncertainty in numerical models (Ghanem  _et al._ ,
2017; Poëtte, Birolleau, and Lucor, 2015). UQ can be achieved through
intrusive approaches that modify deterministic numerical algorithms (Després,
Poëtte, and Lucor, 2013), or non-intrusive techniques employing black box
numerical models (Gilli  _et al._ , 2013). Since solving two-phase Navier-
Stokes equations remains a challenging field (Diwan and Mohamed, 2019), this
study focuses on the non-intrusive method to establish a generic framework for
UQ in computational slurry flows. Stochastic simulations, based on a Monte-
Carlo (MC) approach (Cheng  _et al._ , 2019; Vrugt and Robinson, 2007), are
commonly used to quantify uncertainty without modifying the primary model.
Those are based on the central limit theorem (Geyer, 1994). They rely on a
large number of well-characterized random variable samples in a probabilistic
space, resulting in the convergence of distributions according to their
probability density functions. While Monte-Carlo sampling methods guarantee
convergence, they require a significant number of simulations (Faivre  _et
al._ , 2013). Which may be unfeasible for a large range of applications due to
the computational cost. To overcome the limitations of traditional MC
simulations, numerous approaches have been introduced to perform UQ. Their aim
is to be at least as accurate as MC methods in estimating uncertainty while
requiring less computational time. Meta-models, such as Kriging and Polynomial
Chaos Expansion (PCE), are very useful within this context, by emulating the
behavior of the full model with reduced computational effort. Our past results
(Elkarii  _et al._ , 2023) demonstrate that PCE is particularly efficient for
quantifying uncertainties for one dimensional slurry flows simulations.
Accordingly this surrogate model is chosen in the current study as well.

The PCE process initiates by building surrogate models known as Polynomial
Chaos expansions (PCEs). PCE theory utilizes orthogonal polynomials to
establish the basis (Xiu and Karniadakis, 2002), enabling the computation of
polynomial coefficients that define a reduced model. The intrusive method
modifies the existing code by approximating the model response in the relevant
partial differential equations (Askey and Wilson, 1985). Alternatively, non-
intrusive methods treat the code as a black box, obtaining response solutions
through model evaluations (Askey and Wilson, 1985). Both intrusive and non-
intrusive approaches derive polynomial approximations using ensemble-based
strategies (Ghanem and Spanos, 1991). Originally proposed as an alternative to
Monte Carlo methods, PCE has proven effective in uncertainty quantification
and sensitivity analysis (Wang  _et al._ , 2019). Generalized Polynomial Chaos
(gPC), an extension of PCE, handles probability distributions beyond Gaussian
rules using the Askey scheme (Xiu and Karniadakis, 2002). For convenience, we
use the term PCE throughout this study.

Addressing the difficulties posed by the complex numerical models as part of
the development of meta-models to lower the computing cost in the UQ has
received a lot of interest. This is mostly because the construction time
needed for PCE rises exponentially with the stochastic dimension. One of the
solutions to tackle the problem of stochastic dimensions was to introduce the
adaptive sparse PCE. The objective is to estimate the PCE using a regression
approach using least angle regression, as specified in (Efron  _et al._ ,
2004). Another way to counter this curse of dimensionality was to combine
Kriging with PCE to offer a novel meta-model (Yamazaki  _et al._ , 2018).

This method enabled (Elkarii  _et al._ , 2023) to conduct a UQ and a global
sensitivity analysis (GSA) for one-dimensional quantities of interests (QoI)s
related to the transportation of phosphate slurry through pipelines. The
objective of this study is to expand the dimensionality of our QoIs from one
dimension to two dimensions. Expanding our study to consider two-dimensional
QoIs in UQ and GSA holds significant importance. While one-dimensional
analyses provide valuable insights into specific local aspects of slurry flow
in pipelines, incorporating the second dimension allows for a more
comprehensive understanding of the system’s behavior. By considering both
spatial dimensions, we can capture additional information about variables such
as particle distribution, velocity profiles, which play crucial roles in
optimizing slurry transportation processes. The inclusion of two-dimensional
QoIs enhances our ability to identify and address potential challenges, such
as anisotropic behavior, localized flow variations, or concentration gradients
that may impact the system’s performance. Additionally, a higher-dimensional
analysis facilitates a more accurate representation of the flow phenomena,
allowing us to refine the surrogate models.

The adaptive sparse PCE previously discussed was used in (Elkarii  _et al._ ,
2023), to alleviate the dimension in the model input parameters, however, they
neglect the spatial dimension of the output variables. Classically when the
output solution is a physical field discretized in a mesh, the approach
requires PCE calculation at each node of the mesh, as shown in (El Moçayd  _et
al._ , 2018). In such circumstances, the numerical resolution of complicated
slurry flow characteristics necessitates the use of fine meshes, which incurs
large computing costs. In addition, it should be emphasized that the physical
outputs are stochastic processes with spatial and temporal descriptions, not
only stochastic random variables. Therefore, outputs are computed for each
node in the computational mesh and at every time step in the time interval.
Which necessitates the development of less time-intensive methods for
quantifying uncertainty when the output is a stochastic process. Various
approaches, such as the Principal Component Analysis (PCA) (Blatman and
Sudret, 2011a), the Proper Generalized Decomposition (PGD) (Chevreuil, 2012),
and the Proper Orthogonal Decomposition POD (Raisee, Kumar, and Lacor, 2015),
might be used to address this problem. Indeed, the majority of these
techniques emanate from the domain of data science. Specifically, they fall
under the category of unsupervised learning techniques, facilitating
clustering analysis. In the context of uncertainty quantification, these
clusters represent distinct physical patterns driven by stochastic modes—those
that are the drivers of uncertainty. While their utilization has witnessed a
growing trend within the field, the systematic assessment of their application
to intricate physical problems, such as slurry flows, remains an unexplored
territory until now. It is also important to recognize that these strategies
operate under the assumption that uncertainties can be captured by a few non-
physical factors manifested in the response of random vectors. By resampling
the newly constructed surrogate model, the relevant statistical quantities can
be obtained.

Due to the complex physical characteristics of slurry flow, including steep
gradients and localized structures in response to the two-phase Navier-Stokes
equations, this study employs the POD technique. Our research introduces a
non-intrusive reduced-order approach for simulating the two-dimensional
transport of slurries. The main focus of this study is to analyze the
uncertainty associated with initial and boundary conditions, as well as solid
properties, and its impact on the energy efficiency of slurry pipe flows.
Random environmental factors often affect the slurry transport energy
consumption, leading to uncertainties in measured results. Quantifying these
uncertainties is important to validate deterministic solutions. The specific
energy consumption (SEC), which measures the energy required to move a slurry
batch over a given distance, is used to quantify the transport efficiency. It
is directly proportional to the hydraulic gradient, representing the drop in
static pressure due to friction in the slurry flow per unit length of the
pipe. The proposed methodology focuses on the UQ in terms of pressure drop and
delivered solid concentration, resulting from uncertain modeling and physical
parameters. Additionally, a GSA is conducted using a variance-based method,
specifically Sobol’ indices. To assess the performance of the proposed
methodology, a comparison is made with a class of extensive Monte Carlo
simulations in terms of mean and variance fields.

The remaining sections of this paper are organized as follows. Section III
presents the slurry transport model, the CFD modeling procedure, and the
numerical validation conducted. In Section II, the model reduction methodology
is discussed, along with the PCE and POD techniques employed. Section IV
presents the results obtained for UQ and sensitivity analysis of the QoIs.
Finally, Section V provides concluding remarks based on the findings of this
study.

##  II Uncertainty quantification in stochastic processes

The overall concept is described with a view toward addressing stochastic
slurry pipe flows. Polynomial chaos expansions and stochastic proper
orthogonal decomposition techniques are briefly discussed. The suggested UQ
approaches are also employed to calculate sensitivity indices in the present
study.

###  II.1 Polynomial Chaos Expansion

Within the framework of uncertainty quantification, the PCE has seen extensive
application as a surrogate model, as its objective is to recreate the global
behavior of a simulation following a polynomial decomposition. The latter are
multivariate orthogonal polynomials and serve as the basis functions selected
with respect to the joint probability distributions of the stochastic input
variables following the so-called Askey scheme of polynomial (Xiu and
Karniadakis, 2002). In the present study, if our variables are stored in the
vector $\mathbf{X}\in\mathbb{R}^{d}$ where $d$ is the number of input
variables in $\mathbf{X}$, and $\mathcal{Y}$ represents the model responses
discussed earlier, then one may write, according to PCE:

$\mathcal{Y}=\mathcal{M}(\mathbf{X})=\sum_{\mathbf{\alpha}\in\mathbb{N}^{d}}\mathbf{a}_{\mathbf{\alpha}}\psi_{\mathbf{\alpha}}(\mathbf{X}),$
(1)

where $\mathbf{\alpha}=\left\\{\alpha_{1},\ldots\alpha_{d}\right\\}$ denotes a
$d$-dimensional index,
$\left\\{\mathbf{a}_{\alpha},\mathbf{\alpha}\in\mathbb{N}^{d}\right\\}$ are
the PCE coefficients to be computed. The $\psi_{\mathbf{\alpha}}$ are the
multivariate polynomial basis elements of degree
$|\mathbf{\alpha}|=\sum_{i=1}^{d}\alpha_{i}$ that are orthonormal with regards
to the joint pdf, $f_{\mathbf{X}}$, of $\mathbf{X}$, i.e.,
$\mathbb{E}\left[\psi_{\mathbf{\alpha}}(\mathbf{X})\psi_{\mathbf{\beta}}(\mathbf{X})\right]=1$
if $\mathbf{\alpha}=\mathbf{\beta}$ and 0 otherwise. For practical purposes,
the determination of a PCE is dependent upon the calculation of the spectral
coefficients $\left\\{\mathbf{a}_{\alpha},\alpha\in\mathcal{A}\right\\}$ where
$\mathcal{A}\subset\mathbb{N}^{d}$ is a finite set of multi-indices such that

$\mathcal{Y}\approx\mathcal{M}^{\text{PC}}\left(\mathbf{X}\right)=\sum_{\bm{\alpha}\in\mathcal{A}}\mathbf{a}_{\bm{\alpha}}\psi_{\bm{\alpha}}\left(\mathbf{X}\right).$
(2)

To construct such a set, a standard truncation scheme, which corresponds to
all polynomials in the $d$ input variables of total degree less than or equal
to $p$ is defined by the set $\mathcal{A}^{d,p}$. This set is traditionally
defined as $\mathcal{A}^{d,p}=\left\\{\alpha\in\mathbb{N}^{d}:|\alpha|\leq
p\right\\}$. It is worth noting that when this method is used to define the
set of multi-indices, $\mathcal{A}^{d,p}$ comprises $(d+p)!/d!p!$ items
following the binomial rule for a given value of the maximum polynomial degree
$p$. If the model’s response is strongly nonlinear with numerous input
parameters (requiring a large value of $p$) or if the input vector $X$ is
large (e.g., $M>10$), implementing such a truncation scheme would result in
challenging problems. This issue is commonly referred to as the curse of
dimensionality, we refer the reader to (Xiu and Karniadakis, 2003; Xiu, 2010)
for further details. In the current study, for the purpose of simplicity and
without losing generality, just the regression approach is studied, and the
most efficient sparse decomposition may be used when the stochastic dimension
is very large (Blatman and Sudret, 2011b). It must be emphasized that the
well-established quadrature methods are not examined in the present work since
they cannot be utilized to construct POD modes. To estimate the coefficients
$\mathbf{a}=\left\\{\mathbf{a}_{\mathbf{\alpha}}\in\mathcal{A}\right\\}$, the
regression approach relies on the solution of a least-square (LS) minimization
problem in any $\ell_{2}$-norm (Choi  _et al._ , 2004; Berveiller, Sudret, and
Lemaire, 2006). To put this into application, the expansion coefficients
$\mathbf{a}$ are calculated as:

$\mathbf{a}=\operatorname*{arg\,min}_{\mathbf{a}\in\mathbb{R}^{\mathtt{M}}}\mathbb{E}\left[\left(\mathcal{Y}-\sum_{\mathbf{\alpha}\in\mathcal{A}}\mathbf{a}_{\mathbf{\alpha}}\psi_{\mathbf{\alpha}}(\mathbf{X})\right)^{2}\right],$
(3)

where $\mathtt{M}$ the cardinal of the set $\mathcal{A}$. The error is defined
as the difference between the (exact) model’s evaluation

$\mathcal{Y}=\left\\{\mathcal{M}\left(X_{1}\right),\ldots,\mathcal{M}\left(X_{\mathtt{N}}\right)\right\\},$
(4)

and the PCE surrogate estimations for a finite training set or experimental
design (DE) of $\mathtt{N}$ randomly sampled input varfiables
$\mathbf{X}=\left\\{X_{1},\ldots,X_{\mathtt{N}}\right\\}$ where
$X_{i}\in\mathbb{R}^{d}$. The set of $\mathtt{N}$ realization of the input
vector, $\mathbf{X}$, is then needed, called design of experiment (DE). The
generation of the DE is obtained in the present work using the Latin Hypercube
sampling approach, which helps to reduce the computational cost of the
classical Monte Carlo sampling (Iman, 2008). The LS error is transformed into
a discretized mean-square problem as:

$\widehat{\mathbf{a}}=\operatorname*{arg\,min}_{\mathbf{a}\in\mathcal{R}^{\mathtt{M}}}\frac{1}{\mathtt{N}}\sum_{i=1}^{\mathtt{N}}\left(\mathcal{Y}^{(i)}-\sum_{\mathbf{\alpha}\in\mathcal{A}}\mathbf{a}_{\mathbf{\alpha}}\psi_{\mathbf{\alpha}}\left(X_{i}\right)\right)^{2}.$
(5)

Following the work of Blatman (2009), the regression estimates of the
polynomial chaos coefficients are given by:

$\widehat{\mathbf{a}}=\left(\bm{\Psi}^{\intercal}\bm{\Psi}\right)^{-1}\bm{\Psi}^{\intercal}\mathcal{Y},$
(6)

where $\bm{\Psi}$ is the information matrix of size
$\mathtt{N}\times\mathtt{M}$ give by $\Psi_{ij}=\psi_{j}\left(X_{i}\right)$
assembling the values of all orthonormal polynomials at the ED nodes. In case
of high dimensional problem an attractive modification on the conventional
least squares formulation could be adopted, that limits the sum of the
absolute regression coefficients. In fact, the least-angle regression (LAR),
offers a sparse PCE representation with fewer regressors than the traditional
complete representation. Here, we apply the LAR approach as suggested by
Blatman and Sudret (2011b). For further information on the LAR approach and
its use in the context of adaptive sparse PCE, see (Efron  _et al._ , 2004)
and (Blatman and Sudret, 2011b), respectively. Finally, in order to evaluate
the quality of the PCE, Leave-One-Out (LOO) error is often used as an
excellent measure of accuracy since it enables error estimates without
acquiring further forward evaluation of the numerical model (Blatman and
Sudret, 2010). The definition of the relative LOO error is as follows:

$\epsilon_{\text{LOO}}=\frac{\sum_{i=1}^{\mathtt{N}}\left(\frac{\mathcal{M}\left(X_{i}\right)-\mathcal{M}^{PC}\left(X_{i}\right)}{1-h_{i}}\right)^{2}}{\sum_{i=1}^{\mathtt{N}}\left(\mathcal{M}\left(X_{i}\right)-\hat{\mu}_{\mathcal{Y}}\right)^{2}},$
(7)

where $h_{i}$ is the $i^{th}$ diagonal term of
$\mathbf{\Psi}\left(\mathbf{\Psi}^{\intercal}\mathbf{\Psi}\right)^{-1}\mathbf{\Psi}^{\intercal}$
and
$\hat{\mu}_{Y}=\frac{1}{\mathtt{N}}\sum_{i=1}^{\mathtt{N}}\mathcal{M}\left(X_{i}\right)$
According to a research summarized by Molinaro, Simon, and Pfeiffer (2005),
the leave-one-out method has a low estimate bias and a low mean-square error
and could be comparable with standard relative $\ell_{2}$ error in some cases,
see for example (El Moçayd  _et al._ , 2018).

###  II.2 Stochastic proper orthogonal decomposition

The POD method, enables a high-dimensional system to be addressed by a low-
dimensional one. This approach involves establishing a set of orthogonal
eigenvalues that are representative of the simulated physics. Solving the
integral of Fredholm yields the eigenvectors, while the kernel of this
integral is derived from a series of simulations performed using an
experimental design. Specifically, the eigenfunctions associated with the
problem are optimal within the context of the dynamic representation
(explained below), allowing them to be used to build a simplified
representation of physics. In uncertainty quantification, the POD serves to
decrease the size of a random vector produced by the model. As a result, the
uncertainty is applied for each direction given by the eigenvectors
$\lambda_{i}$. The concept is based on the projection of the model’s solution
$\mathcal{M}\left(X_{i}\right)$ onto a finite and orthonormal basis
$\left\\{\phi_{i},i\in\mathcal{I}_{\text{POD}}\right\\}$, where
$\mathcal{I}_{\text{POD}}$ is a discrete finite collection of indices. Thus,
the procedure $\mathcal{M}\left(X_{i}\right)$ is decomposed as follows:

$\mathcal{M}\left({X}\right)=\sum_{i\in\mathcal{I}_{\text{POD}}}\lambda_{i}({X})\phi_{i},$
(8)

where $\phi_{i}$ is estimated by decomposing the spatial covariance matrix
constructed as:

$\mathbf{C}=\frac{1}{\mathtt{N}}\mathcal{M}\mathcal{M}^{\intercal}.$ (9)

In addition, we construct a POD-truncated error $\varepsilon$ so that only the
$k$ most valuable eigenvectors are maintained as the solution following the
estimate:

$\frac{\sum_{i=0}^{k}\lambda_{i}}{\sum_{i=0}^{\mathtt{N}}\lambda_{i}}>1-\varepsilon,$
(10)

where $\widehat{\lambda}_{i}$ is the mean value of $\lambda_{i}({X}_{i})$.
More details about the POD technique are available in El Moçayd  _et al._
(2020).

###  II.3 POD-PCE meta-model

The eigenvalues are handled as random variables after rebuilding the
stochastic POD. This suggests that we may preserve the spatial dependence
suggested by the eigenvectors $\phi_{i}$ while generating a PCE for each
eigenvalue in the manner described in the section on polynomial chaos
expansions before.

$\lambda_{i}({X})=\sum_{j=0}^{\mathtt{M}}\gamma_{j}\psi_{j}\left({X}_{i}\right).$
(11)

The equation (1) becomes,

$\mathcal{M}\left({X}\right)=\sum_{i\in\mathcal{I}_{POD}}\left(\sum_{j=0}^{\mathtt{M}}\gamma_{j}\psi_{\alpha}\left({X}\right)\right)\phi_{i}.$
(12)

A summary of both algorithms retained in the present study for quantification
of uncertainties in stochastic slurry flows is depicted in Fig. 17.

Figure 1: Diagram illustrating the distinction between the traditional PCE-
based surrogate model and the POD-based surrogate model.

It should be noted that the traditional method of dealing with an uncertainty
quantification issue utilizing the PCE is conducted when a decomposition is
completed for each node in the computational mesh. This procedure is shown in
Fig. 17 flowchart by steps 1 through 3. However, as previously stated, one way
to ease the spatial distribution is to do a first reduction using the POD. The
stochasticity is assumed to be incorporated in the related eigenvalues, this
latter will aid in distinguishing the spatial dependence from the stochastic
dependency. Once the POD is performed, just a few eigenmodes are preserved,
allowing for less PCE to be generated than its traditional equivalent. This
procedure is shown in Fig. 17 flowchart as stages 1, 2, 4, and 5. The accuracy
of the approaches under consideration is evaluated by comparing their findings
to large MC simulations (with N=1000 samples) for the quantities of interest
in the mean and variance fields. This is generally achieved using a
$l_{2}$-error following this:

$l_{2}(\theta)=(\theta_{\text{MC}}-\theta_{\text{POD-PCE}})^{2}$ (13)

where the $\theta$ is the required quantity of interest (Mean value or
variance value), $\theta_{\text{MC}}$ refer to its estimation using Monte-
Carlo simulations and $\theta_{\text{POD-PCE}}$ refer to its estimation using
the surrogate model.

###  II.4 Aggregated Sobol’ indices

The sensitivity indices in this paper are derived using the PCE approach with
a non-intrusive sampling procedure. Therefore, once the PCE is built, the mean
$\mu$ and the total variance $D$ can be obtained using the value of the
spectral coefficients such that $\mu=a_{0}$ and
$D=\sum_{\alpha\in\mathcal{A}\setminus 0}a_{\mathbf{\alpha}}^{2}$. The
computation of Sobol’ indices of the output variable
$\mathcal{M}_{j}\left({X}_{i}\right)$ concerning the input variable ${X}_{i}$
of any order is a straightforward process. The first order Sobol’ indices can
be expressed as follows
$\mathtt{S}_{i}^{j}=\sum_{\mathbf{\alpha}\in\mathcal{A}_{i}}a_{\mathbf{\alpha}}^{2}/D$
with
$\mathcal{A}_{i}=\left\\{\mathbf{\alpha}\in\mathcal{A}:\alpha_{i}>0,\alpha_{j\neq
i}=0\right\\}$, while the total Sobol’ indices is computed as
$\mathtt{S}^{T}_{i,j}=\sum_{\mathbf{\alpha}\in\mathcal{A}_{i}^{T}}a_{\mathbf{\alpha}}^{2}/D$
with
$\mathcal{A}_{i}^{T}=\left\\{\mathbf{\alpha}\in\mathcal{A}:\alpha_{i}>0\right\\}$.
Since the QoIs under consideration in the present studies have spatial
dependencies, sensitivity indices should be defined accordingly. The
aggregated first Sobol’ index $GS_{i}$ as introduced in (Blatman, 2009) are
defined here:

$\mathtt{GS}_{i}=\frac{\sum_{j=1}^{p}\operatorname{Var}\left(\mathcal{M}\left({X}_{j}\right)\right)S_{i}^{j}}{\sum_{j=1}^{p}\operatorname{Var}\left(\mathcal{M}\left({X}_{j}\right)\right)}$
(14)

The aggregated total Sobol’ index are calculated as:

$\mathtt{GS}^{T}_{i}=\frac{\sum_{j=1}^{p}\operatorname{Var}\left(\mathcal{M}\left({X}^{(j)}\right)\right)S^{T}_{i,j}}{\sum_{j=1}^{p}\operatorname{Var}\left(\mathcal{M}\left({X}^{(j)}\right)\right)}$
(15)

Notice that $\mathtt{GS}_{i}\in[0,1]$. As in the case of the scalar indices,
the main advantage of the aggregated indices is their easy interpretation: a
high value means the input is important, a value close to zero means it is
not, and $\sum_{i=1}^{N}\mathtt{GS}_{i}=1$.

##  III Computational Model (CFD modeling procedure)

The system under consideration is a three-dimensional pipeline with a diameter
of $D=0.103$ m and a length of $L=140D$ to ensure fully developed flow. The
transported slurry consists of a mixture of phosphate and water, with
densities of $\rho_{s}=2600$ kg/m3 and $\rho_{l}=1000$ kg/m3, respectively.
The actual flow of phosphate slurry operates under specific conditions.
Regarding the flow rate, the pipeline has a capacity of approximately $4400$
t/h, depending on the properties of the slurry. This corresponds to a mixture
inlet velocity of $u_{m}=2.16$ m/s, and it can accommodate a maximum deposit
limit velocity of ${u}_{dl}=1.75$ m/s. The delivered phosphate pulp has a
volumetric solid concentration of approximately $\phi_{s}=33.3\%$, and it
contains rigid solid particles with random shapes, having equivalent spherical
diameters ranging from $44~{}\mathrm{\mu}$m to $250~{}\mathrm{\mu}$m. For
further information about the transportation process and the operating limits
of the pipeline, refer to the work of Rusconi, Lakhouaja, and Kopuz (2016).

###  III.1 CFD numerical modeling

The two-phase flow is modeled using the Eulerian-Eulerian approach implemented
in the twoPhaseEulerFoam solver available in OpenFOAM (Weller  _et al._ ,
1998). The mass and momentum balance equations for the liquid ($k\equiv l$)
and solid ($k\equiv s$) phases are solved with closure models for the
turbulent stress tensor and interphase forces. The mass and momentum balance
equations are the following:

$\frac{\partial\rho_{k}\phi_{k}}{\partial
t}+\nabla\cdot\left(\rho_{k}\phi_{k}\mathbf{u}_{k}\right)=0,$ (16)
$\frac{\partial\rho_{k}\phi_{k}\mathbf{u}_{k}}{\partial
t}+\nabla\cdot\left(\rho_{k}\phi_{k}\mathbf{u}_{k}\mathbf{u}_{k}\right)=-\nabla\cdot\left(\phi_{k}\mathbf{\tau}_{k}\right)-\phi_{k}\nabla
p+\phi_{k}\rho_{k}\mathbf{g}+\mathbf{F}_{k},$ (17)

where $\rho$ is the density, $\alpha$ is the phase fraction, $\mathbf{u}$ is
the average velocity, $p$ is the average pressure field shared by all the
phases, and $\mathbf{g}$ is the gravity field. The inter-phase momentum
transfer term is denoted as $\mathbf{F}_{k}$, accounting for the drag, lift,
virtual mass, and turbulent drag forces. For the dispersed solid phase, the
inter-phase momentum transfer term takes the following form:

$\mathbf{F}_{s}=\mathbf{F}_{d,s}+C_{l}\rho_{l}\mathbf{u}_{r}\times\left(\nabla\times\mathbf{u}_{c}\right)+C_{vm}\rho_{l}\left(\frac{D_{l}\mathbf{u}_{l}}{Dt}-\frac{D_{s}\mathbf{u}_{s}}{Dt}\right)-C_{td}\alpha_{s}\rho_{c}\kappa_{m}\nabla\alpha_{s},$
(18)

Here, $\kappa_{m}$ represents the multiphase mixture of turbulent kinetic
energy. The interfacial momentum balance implies that
$\mathbf{F}_{s}=-\mathbf{F}_{l}$. In Eq. (18), $\mathbf{F}_{d,s}$ is the drag
force acting in the solid phase. The symmetric drag formulation proposed by
Weller (Weller, 2002) is applied:

$\mathbf{F}_{d,s}=\alpha_{s}\alpha_{l}\left(f_{s}\frac{C_{d,l}\rho_{s}}{d_{s}}+f_{l}\frac{C_{d,s}\rho_{l}}{d_{l}}\right)|\mathbf{u}_{r}|\mathbf{u}_{r},$
(19)

where

$f_{k}=\min\left[\max\left(\frac{\phi_{k}-\phi_{k}^{\max,f}}{\phi_{k}^{\max,p}-\phi_{k}^{\max,f}},0\right),1\right],$

where $\alpha_{l}^{\max,p}=\alpha_{s}^{\max,p}=0.63$. In this formulation,
$\mathbf{u}_{r}=\mathbf{u}_{s}-\mathbf{u}_{l}$ represents the relative
velocity between phases. The coefficients $C_{d}$, $C_{l}$, $C{vm}$, and
$C{td}$ correspond to the drag, lift, virtual mass, and turbulent drag
coefficients, respectively, with values set as $C_{l}=C_{vm}=C_{td}=0.5$. The
drag coefficients are determined using the Schiller and Naumann correlation
(Schiller, 1933a, b):

$C_{d,k}=\max\left[0.44,\frac{24}{\mathrm{Re}_{k}}\left(1+0.15\mathrm{Re}_{k}^{0.687}\right)\right],~{}$

where the Reynolds numbers based on the solid and liquid phase diameters,
$\mathrm{Re}_{s}$ and $\mathrm{Re}_{l}$ respectively, are defined as
$\mathrm{Re_{s}}=\rho_{l}\mathbf{u}_{r}d_{s}/\mu_{l}$ and
$\mathrm{Re_{l}}=\rho_{l}\mathbf{u}_{r}d_{l}/\mu_{l}$, where $d$ represents
the dispersed phase diameter and $\mu$ is the dynamic viscosity. The total
stress tensor per unit mass, denoted as $\mathbf{\tau}_{k}$, is calculated
using a Reynolds-Averaged Navier-Stokes approach for the turbulence modeling
of the two-phase mixture, as described in (Behzadi, Issa, and Rusche, 2004):

$\mathbf{\tau}_{k}=\mu_{\mathrm{eff},k}\left[\nabla\mathbf{u}_{k}+\left(\nabla\mathbf{u}_{k}\right)^{\intercal}-\frac{2}{3}\mathbf{I}\nabla\cdot\mathbf{u}_{k}\right],$
(20)

where $\mu_{\mathrm{eff},k}$ represents the effective phase viscosity,
combining the molecular and turbulent dynamic viscosity of the phase, given by
$\mu_{\mathrm{eff},k}=\mu_{k}+\mu_{t,k}$. However, in relation to the
viscosity of the solid phase, it’s necessary to use closures to represent how
the behavior of this phase, encompassing aspects like contact pressure and
velocity collisions, translates into continuous pressure and viscosity. This
study employs the Kinetic Theory of Granular Flow (KTGF), which establishes a
relationship between the overall kinetic energy of a group of particles moving
randomly within the carrier fluid and their varying velocity, employing the
concept of granular temperature. The viscosity coefficients encompass a
combination of three distinct elements:

$\mu_{s}=\mu_{s,\mathrm{col}}+\mu_{s,\mathrm{kin}}+\mu_{s,\mathrm{fric}},$
(21)

where $\mu_{s,\mathrm{col}}$, $\mu_{s,\mathrm{kin}}$, and
$\mu_{s,\mathrm{fric}}$ represent the collisional Gidaspow, Bezburuah, and
Ding (1991), kinetic Syamlal, Rogers, and O’Brien (1993), and frictional
Schaeffer (1987) components of the total shear stress, respectively, as
derived from KTGF theory. Further elaboration on the formulation of each term
can be found in Ding and Gidaspow (1990); Elkarii  _et al._ (2022).

The turbulent viscosity is modeled using the $k-\omega$ SST model (Wilcox,
1988) in order to guarantee mixture’s pseudo-homogeneity. For large scale
slurry flows, the $k-\omega$ SST turbulence model should be natural as it
properly predicts turbulent variables in the vicinity of walls thanks to the
$k-\omega$ model and it allows to achieve turbulent regime in the bulk is
achieved by means of $k-\varepsilon$ model. It is worth noting that the
$k-\omega$ SST turbulence model is employed with an enhanced wall function
treatment that blends laminar and turbulent laws-of-the-wall to provide
continuous and asymptotically correct values for all $y^{+}$. A mesh
refinement is performed in the near-wall regions. A priori estimate of shear
stresses are used to determine the wall-adjacent mesh cell heights required to
obtain $y^{+}\approx 1$ (cf. Fig. LABEL:fig:Section_mesh), where
$y^{+}={u}^{\ast}x_{2}/\nu$ is a function of the friction velocity
${u}^{\ast}=\sqrt{\tau_{\mathrm{wall}}/\rho}$. The wall shear stress
$\tau_{\mathrm{wall}}$ is estimated with the friction factor approximated by
the Swamee and Jain (1976)’s equation and considering the root-mean-squared
value of the axial and circumferential Reynolds numbers with zero wall
roughness and the average axial velocity. To enable an effective capturing of
the viscous sublayer, the value of $y^{+}$ is kept less than 1. The growth
rate of the mesh is set at $1.1$ until it reaches the log-law layer,
equivalent to $y^{+}=60$. A total of $700\times 10^{3}$ grid meshes are used
to discretize the whole three-dimensional computational domain, as illustrated
in Fig. LABEL:fig:3D_domain.

(a)

(b)

Figure 2: (a) Three dimensional view of the computational domain and boundary
conditions, (b) detail of the discretization of the pipe section and local
magnification of the mesh close to the near-wall cells.

The solver has been able to predict particle profiles for very concentrated
sand-water mixtures (Elkarii, Bouallou, and Ratnani, 2020). Computational
meshes are generated using Open$\nabla$FOAM®’s blockMesh utility. The
simulations are performed utilizing an implicit-Euler scheme to discretize in
time; Gradient terms are discretized using a cell limited linear Gauss scheme;
The second order corrected linear Gauss scheme is employed for Laplacian
terms; The interpolation scheme is linear; For the surface-normal gradient
scheme, a corrected scheme is used, and for the wall distance computation, a
mesh-wave is used. Table 1 summarizes the major solvers employed in the
generation of the present database. The reader is referred to (Elkarii  _et
al._ , 2022) for more details of the CFD modeling approach and numerical
validation of the CFD model.

Table 1: Solver scheme used in the simulations Term |  TwoPhaseEulerFoam scheme
---|---
Convection & Momentum transport | Gauss Van Leer
Turbulence & Artificial compression | Gauss upwind

##  IV Numerical results

This section assesses the computational efficiency of the proposed POD
approach for analyzing slurry flows. For each test case, we compare the
results of the mean and variance of the potential solution using two methods:
MC and POD with POD-PCE. The two-dimensional plan considered in the present
study is the mid-plane of the pipe that is a squared domain $\Omega\equiv
L_{x_{1}}\times L_{x_{2}}=[0,140D]\times[-0.5D,0.5D]$. Consequently, we
extract the two-dimensional two-phase flow solution by employing the
deterministic solution obtained from the three-dimensional two-phase Navier-
Stokes equations. In this context, our objective is to evaluate the
uncertainty associated with the QoIs, which include the solid concentration,
velocity, and pressure distribution over the mid-plane pipe.

###  IV.1 Two-dimensional solid concentration

Initially, we executed 1000 three-dimensional CFD simulations based on a
designed experiment, where input parameters were identified and sampled using
Latin hypercube sampling to ensure uniform coverage across the designated
parameter space for each input. Subsequently, we obtained the solution on a
structured mesh consisting of a total of $18894$ nodes by extracting a mid-
plane section. Next, surrogate models were constructed using POD-PCE
methodology. It is important to highlight that, unlike the conventional PCE
method applied to the complete model, the POD-PCE method solves a reduced
model. The accuracy of the POD-PCE method for each eigenvalue is evaluated and
presented in Table 2.

Table 2: Optimal polynomial degrees with LOO errors for the POD modes in the slurry particles concentration ${\phi_{s}}$ Number of modes | LOO error | Polynomial degree
---|---|---
Mode #1 | $1.2758\times 10^{-5}$ | $6$
Mode #2 | $5.7266\times 10^{-4}$ | $8$
Mode #3 | $3.1746\times 10^{-3}$ | $7$
Mode #4 | $1.0955\times 10^{-2}$ | $10$
Mode #5 | $3.3131\times 10^{-2}$ | $10$
Mode #6 | $1.3435\times 10^{-1}$ | $10$
Mode #7 | $2.9421\times 10^{-1}$ | $10$
Mode #8 | $1.0350\times 10^{-1}$ | $10$
Mode #9 | $4.7471\times 10^{-1}$ | $10$
Mode #10 | $7.8437\times 10^{-2}$ | $9$
Mode #11 | $6.7515\times 10^{-1}$ | $10$
Mode #12 | $1.1624\times 10^{-1}$ | $10$

In this particular scenario, employing a truncation criterion of
$\varepsilon=10^{-4}$, the model selectively retains only the initial $12$
eigenvalues, and the PCE process is carried out individually for each mode. As
a result, the surrogate model significantly reduces computational
requirements, necessitating only $12$ decompositions instead of the original
$18894$ associated with the total number of nodes in the two-dimensional
numerical model. Table 2 showcases the optimal polynomial degrees employed to
approximate each spectral mode using the POD technique, alongside an
evaluation of the Leave-One-Out (LOO) error generated by the PCE method when
estimating the numerical solution. Upon analyzing the outcomes, it is evident
that the problem exhibits nonlinearity, as high polynomial degrees are
required to effectively capture the uncertainty. Furthermore, the findings
indicate satisfactory magnitudes for the LOO errors. Moving forward, in Fig.
3, we compare the mean solid concentration solutions obtained from the
stochastic simulation with the deterministic exact solution, while also
presenting the discrepancy between the two solutions for easy visual
comparison.

(a)

(b)

(c)

Figure 3: Mean solid concentration $\phi_{s}$ obtained for the stochastic
simulation (a), the deterministic exact solution (b) and the difference
between the two solutions (c) obtained for a slurry flow problem.

The results depicted in Fig. 3 illustrate that, considering the level of
uncertainty examined in the problem, the mean solutions computed by both
methods exhibit similar trends and patterns. The stochastic approach
successfully addresses the distribution of solid concentration, producing
numerical outcomes that are devoid of any nonphysical oscillations. In
general, these findings imply that the model reduction has minimal impact on
the accuracy of UQ. Furthermore, Fig. 4 showcases the obtained results for the
variance.

(a)

(b)

(c)

Figure 4: Same as Fig. 3 but for the variance.

The results for mean and variance make it evident that the two methods under
consideration exhibit identical flow patterns and accurately capture the solid
concentration over the two-dimensional pipe. To provide a more comprehensive
understanding, Fig LABEL:fig:L2_MEAN_C and Fig. LABEL:fig:L2_VAR_C illustrate
the logarithmic $\ell_{2}$-error in estimating the statistical moments
considered. The remarkably low magnitude of the error $(<10^{-10})$ clearly
indicates that the surrogate models proficiently capture the uncertainty
manifested over the physical domain for the slurry flow problem under
consideration. An additional benefit is that the proposed POD-PCE method
attains this high level of accuracy while incurring significantly lower
computational costs in comparison to the MC and PCE methods. For a better
insight, we display in Fig. 5 vertical cross-sections of the mean and variance
at $x_{1}=126D$.

Figure 5: Cross-sections of the mean and variance solutions of the solid
concentration distribution using PCE-POD

The distribution of the mean solid concentration reflects the general behavior
of the flow in the case of phosphate slurry. The findings indicating that the
average concentration of particles in the phosphate slurry flow is stratified
into two distinct layers. This confirms that the slurry flow is of the
settling type. In the core and lower wall of the pipe where particles move
slowly, the solid concentration variance is insignificant. However, the top
wall exhibits the highest variance values due to the relatively sparse flow
and increased particle movement and kinetic energy.

The ability of the surrogate model to estimate the whole Probability Density
Function (PDF) of the solids concentration is also investigated. This
comparison is shown in Fig. 6 where the estimation of the whole PDF using
Monte-Carlo simulation are compared to those obtained using the PCE-POD
surrogate model. PDF estimation is done here using a non-parametric estimation
following a kernel smoothing methodology.

Figure 6: Visualization of solid concentration $(\phi_{s})$ PDF through Kernel
smoothing (KS) method, estimated using MC and PCE-POD techniques.

It is clear that the results demonstrate that there is a good estimation of
the probability density function using the tools presented in this study. The
results demonstrate a robust convergence in the decomposition, rendering them
suitable for UQ in the distribution of solid concentration in slurry flows.
The general shape of the PDF exhibits also the non-linearity effects, as it is
not suited to a Gaussian-like shape. One can also see the presence of two
modes. This may be due to segregation effects. The slurry flow experiences
segregation phenomena during transportation and the PDF of the slide
concentration indicates the formation of two distinct modes in the
concentration distribution.

The sensitivity of the solid concentration field to the variability of the
random parameters $\left\\{\phi_{s},u_{m},d_{p},\theta,\mathcal{SC}\right\\}$,
can be assessed by means of bar-plots of the first order and total aggregated
Sobol’ indices (cf. Fig. 7). It is important to mention that a sensitivity
analysis using only the PCE method has already been conducted for one-
dimensional and scalar QoIs (Elkarii  _et al._ , 2023). The corresponding
Sobol’ indices results will be presented alongside the outcomes obtained
through the POD-PCE method, with the intention of facilitating a comprehensive
comparison. The objective of this comparison is to assess whether increasing
the dimensionality of our QoIs has any impact on the obtained results.

(a)

(b)

Figure 7: Barplots of first and total Sobol’ indices for solid concentration
field calculated based on: (a) One-dimensional data, (b) Two-dimensional data

The reduced model generates a comparable ranking of the most influential
parameters; however, there are variations in the magnitude of the Sobol’
indices. The Sobol’ indices analysis conducted indicates that the distribution
of solid concentration along the vertical diameter of the pipe remains
primarily influenced by the initial concentration introduced into the
pipeline, followed by the particle sizes. Thus, these two parameters are
responsible for regulating the behavior and pattern of the slurry flow.
Indeed, the inlet concentration of a slurry directly impacts the mixture
density, and it plays a significant role in the behavior of particles,
including their collisions and friction, ultimately leading to their
arrangement inside the pipeline. The behavior of particles is also induced by
their size, as larger particles generally cause sedimentation, resulting in a
significant change in the in situ solid concentration in the pipeline.

###  IV.2 Two-dimensional solid velocity

Our second QoI consists of a two-dimensional velocity field of particles.
Similar to the previous example, we have provided the optimal polynomial
degree distribution for approximating each spectral mode using the POD
technique (cf. Table 3). Additionally, we have included the Leave-One-Out
(LOO) error associated with the PCE method for estimating the numerical
solution for each mode.

Table 3: Optimal polynomial degrees with LOO errors for the POD modes in the slurry particles velocity ${u_{s}}$ Number of modes | LOO error | Polynomial degree
---|---|---
Mode #1 | $4.4253\times 10^{-5}$ | $7$
Mode #2 | $2.0990\times 10^{-3}$ | $9$
Mode #3 | $1.1500\times 10^{-2}$ | $6$
Mode #4 | $2.0464\times 10^{-2}$ | $7$
Mode #5 | $3.9424\times 10^{-2}$ | $8$

By using a set of $5$ eigenmodes, the physical model can be significantly
reduced. This means that the surrogate model for two-dimensional solid
velocity can be represented by only the first $5$ modes, rather than the
$18894$ required to represent the total number of nodes in the two-dimensional
numerical model.

The ability of the proposed POD-PCE method for recovering the particle
velocity of a plane pipe is examined. We showcase in Fig. 8 the results
estimated using MC and the proposed surrogate model for the mean solid
velocity, as well as , the difference between them.

(a)

(b)

(c)

Figure 8: Mean solid velocity ${u}_{s}$ obtained for the stochastic simulation
(a), the deterministic exact solution (b) and the difference between the two
solutions (c) obtained for a slurry flow problem.

Under the considered conditions of the slurry flow, both methods show
comparable trends in the velocity profiles of particles. Fig. 9 presents the
comparative results obtained for the variance.

(a)

(b)

(c)

Figure 9: Same as Fig. 8 but for the variance.

The magnitude of logarithmic $\ell_{2}$-error being less than $(<10^{-6})$
indicates that both methods show comparable velocity profile trends for the
solid phase. Once again, the results that have been presented demonstrate a
strong convergence in the decomposition process, making them trustworthy for
UQ purposes. This confirms that the proposed POD-PCE approach is effective in
resolving the issue of the two-dimensional velocity distribution of particles.
Fig. 10 depicts the vertical cross-sections of the mean and variance at
$x_{1}=126D$.

Figure 10: Cross-sections of the mean and variance solutions of the solid
velocity distribution using PCE-POD

The presence of an asymmetrical profile in the distribution of mean solid
velocity is likely a contributing factor to the observed differences in solid
concentration variance across the various walls of the pipe. Specifically, the
skewness of the velocity profile towards the top wall indicates that the flow
in that region is relatively more diluted, and the particles have a higher
kinetic energy. This is consistent with the highest variance values observed
in that area.

The surrogate model ability to estimate the whole Probability Density Function
(PDF) of the solids velocity is also checked here. This comparison is shown in
Fig. 11 where the estimation of the whole PDF using Monte-Carlo simulation are
compared to those obtained using the PCE-POD surrogate model. Once again, the
nonlinear behavior of the model and its impacts on the velocity distribution
is emphasized. This considers a further assessment of the ability of the
surrogate model to perform UQ correctly.

Figure 11: Visualization of solids velocity $(u_{s})$ PDF through Kernel
smoothing (KS) method, estimated using MC and PCE-POD techniques.

Similar to $\phi_{s}$ the sensitivity of the solid velocity field to the
variability of our random parameters can also be evaluated through barplots of
the first-order and total Sobol’ indices.

(a)

(b)

Figure 12: Barplot of first and total Sobol’ indices for solid velocity
calculated based on: (a) One-dimensional data, (b) Two-dimensional data

The velocity distribution of solids in the pipe, as depicted in Fig. 12, is
primarily influenced by the inlet flow rate. On the other hand, the
utilization of the POD-PCE technique reveals slightly dissimilar outcomes,
indicating that the velocity distribution of particles is relatively more
affected by their sizes and also slightly by the pipe inclination. Based on
their values, the observed discrepancies could be attributed to the stochastic
nature of the methodology. Normally, smaller particles have higher velocities
than larger particles at the same temperature. This is because smaller
particles have a higher surface area-to-volume ratio, which means they have
more contact with the surrounding fluid, allowing them to exchange energy and
momentum more efficiently. Additionally, smaller particles experience less
drag and resistance from the surrounding fluid, which can also contribute to
their higher velocity.

###  IV.3 Two-dimensional pressure distribution

The pressure distribution within the slurry flow is the third QoI that was
addressed in this study. Table 4 summarizes the results obtained for PCE-POD
modes.

Table 4: Optimal polynomial degrees with LOO errors for the POD modes in the slurry flow pressure $P$ Number of modes | LOO error | Polynomial degree
---|---|---
Mode #1 | $6.0655\times 10^{-3}$ | $6$
Mode #2 | $3.1614\times 10^{-3}$ | $6$

It can be observed that the model size has significantly decreased, as the
pressure distribution of the slurry flow can be reconstructed using only two
modes instead of the $18894$ that were required for obtaining the
deterministic solution.

Similar to previous analyses, we compare the average pressure distribution
obtained from both MC and POD-PCE methods to evaluate the capability of the
proposed POD-PCE approach in reconstructing the pressure of the two-
dimensional pipe flow. One should note that the results displayed here for the
pressure are calculated in bar unit. The outcomes are depicted in Fig. 13.

(a)

(b)

(c)

Figure 13: Mean flow pressure $P$ obtained for the stochastic simulation (a),
the deterministic exact solution (b) and the difference between the two
solutions (c) obtained for a slurry flow problem.

There are no significant changes between the two solutions, it is further
confirmed by the logarithmic $\ell_{2}$-error estimation. Since, the error’s
magnitude depicted here is less than $10^{-8}$.

Fig. 14 presents the comparative results obtained for the variance.

(a)

(b)

(c)

Figure 14: Same as Fig. 13 but for the variance.

The surrogate models effectively capture the uncertainty observed across the
physical domain in the slurry flow problem under consideration. Given the
relatively low computational effort required to implement the PCE-POD method,
it can be deemed an ideal algorithm for conducting uncertainty quantification
in computational slurry flows involving stochastic inputs. By analyzing the
distribution of the mean flow pressure, it can be deduced that the pressure of
the slurry typically decreases within the range of variability of the random
inputs. This pressure evolution along the pipeline is evident from the
stratified flow pattern discussed in subsection IV.1. Such slurry flows
exhibit fluctuations in particle concentration and velocity, resulting in
frictional resistance and subsequent pressure decrease. The highest variance
values are observed at the inlet of the pipe because the pressure at the inlet
varies in conjunction with the velocity.

The ability of the surrogate model to estimate the entire Probability Density
Function (PDF) of the flow’s pressure drop is examined. This comparison is
depicted in Fig. 15, where the estimation of the complete PDF using Monte
Carlo simulation is compared to that obtained using the PCE-POD surrogate
model.

Figure 15: Visualization of the flow’s pressure drop PDF through Kernel
smoothing (KS) method, estimated using MC and PCE-POD techniques.

Sobol’ indices are utilized to evaluate the sensitivity of slurry flow
pressure. The first and total aggregated Sobol’ indices are shown in Fig. 16.

(a)

(b)

Figure 16: Barplots of first and total Sobol’ indices for the flow’s pressure
drop calculated based on: (a) One-dimensional data, (b) Two-dimensional data

The reduced model exhibits the same ranking of the most influential parameters
as the full-order model. The results clearly indicate that the velocity and
pipe inclination have the greatest impact on the pressure drop, followed by
particle size, particle concentration, and the specularity coefficient. When
considering the range of pipe inclination variations, it is concluded that it
only affects the magnitude of the pressure drop and does not result in a
pressure gain when the pipe is inclined downward. The particle size influences
the friction along the lower wall, and as the slurry flows through the
pipeline, the different layers of particles interact with the component walls
in varying ways. Moreover, the layer with the highest particle concentration
experiences greater frictional resistance against the wall, leading to a
larger pressure drop compared to the layer with lower particle concentration.
Consequently, a pressure gradient forms across the slurry, where the pressure
decreases as the particle concentration increases.

##  V Conclusion

In this study, we have presented and discussed an efficient analysis of UQ for
slurry flows problems. Slurry flows are complex and highly nonlinear, making
UQ studies valuable for gaining physical insights into the underlying
processes. Our work focuses on two-dimensional QoIs problems, where we have
applied a reduced-order model to the two-dimensional Two-Fluid model Navier-
Stokes equations. The adopted strategy combines PCE with the POD technique,
which allows to selectively retain the most relevant PCEs that best capture
the underlying flow physics. Consequently, computational time is significantly
reduced, enabling computations to be executed within a few minutes. A summary
of the proposed methodology developed in this work is presented in Fig. 17.

Figure 17: The workflow of the UQ analyses based on POD-PCE technique.

The numerical results demonstrate that our developed reduced-order model is
capable of producing acceptable results for statistical quantities of
interest. We compare the mean and variance distributions obtained from the
reduced-order model with those computed using deterministic simulations,
revealing that the variability of each QoI is well represented. This
comparison provides valuable information about the distribution of the
outputs.

Our analysis also evaluates the influence of uncertain parameters on the QoIs,
shedding light on the energy efficiency of slurry transport. Notably, we
observed a noticeable impact on the results when increasing the dimensionality
of the QoIs from one to two dimensions. Incorporating two-dimensional QoIs in
our UQ and GSA study allows for a more robust and comprehensive assessment of
uncertainties and sensitivities in complex slurry flow problems. Specifically,
the additional information gained regarding the sensitivity of our QoIs to
input dimensions is crucial. Parameters with high Sobol’ indices play a
significant role in ensuring accurate measurements and should be prioritized.
By determining appropriate measurement uncertainty for these parameters and
understanding the level of relaxation in measurement sensitivity for others,
engineers can mitigate nonlinear effects caused by stochastic measurement
errors and ultimately reduce uncertainty. Furthermore, the inherent
uncertainty in field measurements can introduce uncertainty into CFD
simulations aimed at predicting the specific energy consumption (SEC).
Prioritizing simulation parameters with high uncertainty levels helps reduce
measurement uncertainty and enhances simulation reliability. Therefore,
controlling parameters that have a significant impact on the overall
uncertainty of the simulation is important for improved control and
optimization of slurry flow, leading to more reliable and accurate
simulations.

Our statistical approach surpasses the classical "one-at-a-time" sensitivity
analyses by employing GSA and UQ. This approach considers the joint effects of
multiple parameters, providing a more comprehensive understanding of slurry
flows in pipelines, a complex system. By demonstrating the effectiveness of
these techniques, our study serves as a prototype, paving the way for future
research in the field of slurry flow analysis and optimization.

### Conflict of Interest

The authors have no conflicts to disclose.

### Data availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

### Acknowledgement

This work was supported by OCP Group (Morocco). The authors gratefully
acknowledge the support and computing resources from the African
Supercomputing Center (ASCC) at UM6P (Morocco).

##  Appendix A Solver validation and verification

###  A.1 Numerical validation

The present modeling approach undergoes validation against experimental data
obtained from a sand-water slurry outlined in the literature Gillies, Shook,
and Xu (2004). Comparisons between experimental and simulation data for
particle concentrations at $\phi_{{s}}=0.19$ and $\phi_{{s}}=0.4$ along the
vertical diameter direction across the horizontal pipe cross-section (cf. Fig.
18(b) and Fig. 19(b)) are depicted in Fig. 18(a) and Fig. 19(a), respectively.

(a)

(b)

Figure 18:  Comparison of radial distributions of the volume fraction of sand
particles $\varepsilon_{{s}}$: (a) $\mathrm{d}_{{p}}=0.09$ mm and ${u_{m}}=3$
m/s, (b) Cross section at the pipe outlet.

(a)

(b)

Figure 19:  Comparison of radial distributions of the volume fraction of sand
particles $\phi_{{s}}$: (a) $\mathrm{d}_{{p}}=0.27$ mm and ${u_{m}}=5.4$ m/s,
(b) Cross-section at the pipe outlet

The numerical model’s performance is further analyzed by comparing how the
pressure drop per unit length changes with respect to the inlet velocity of
the mixture, ${u}_{m}$, ranging from $2$ m/s to $8$ m/s. The results are
illustrated in Fig. 20.

(a) $\phi_{{s}}=0.19$, $\mathrm{d}_{{p}}=0.09$ mm

(b) $\phi_{{s}}=0.40$, $\mathrm{d}_{{p}}=0.27$ mm

Figure 20:  Predicted pressure drop per unit of length $\Delta{p}/L$ versus
experimental data. Differences in relative errors between predictions and
experiments are also displayed.

Our CFD results show a significant alignment with the experimental data
provided. The overall mean relative errors in pressure drops remained below
10$\%$, a level widely deemed acceptable. This proves the validity of the
numerical model within the specified range of physical conditions.

## References

  * Dixit  _et al._ (2022) S. Dixit, S. Kumar, A. R. Pradhan, S. Kumar, K. Kumar, N. I. Vatin, T. Miroshnikova, and K. Epifantsev, “Numerical simulation of sand–water slurry flow through pipe bend using CFD,” International Journal on Interactive Design and Manufacturing (IJIDeM) , 1–13 (2022).
  * Sontti  _et al._ (2023) S. G. Sontti, M. Sadeghi, K. Zhou, E. Zheng, and X. Zhang, “Computational fluid dynamics investigation of bitumen residues in oil sands tailings transport in an industrial horizontal pipe,” Physics of Fluids 35, 013340 (2023).
  * Wan  _et al._ (2023) C. Wan, S. Xiao, D. Zhou, H. Zhu, Y. Bao, K. Kakanda, and Z. Han, “Numerical simulation on transportation behavior of gradated coarse particles in deep-sea vertical pipe transportation,” Physics of Fluids (2023).
  * Zhou  _et al._ (2023) H. Zhou, J. Guo, T. Zhang, M. Li, T. Tang, and H. Gou, “Eulerian multifluid simulations of proppant transport with different sizes,” Physics of Fluids (2023).
  * Joshi, Parkash, and Krishan (2023) T. Joshi, O. Parkash, and G. Krishan, “Slurry flow characteristics through a horizontal pipeline at different Prandtl number,” Powder Technology 413, 118008 (2023).
  * El Moçayd and Seaid (2021) N. El Moçayd and M. Seaid, “Data-driven polynomial chaos expansions for characterization of complex fluid rheology: Case study of phosphate slurry,” Reliability Engineering & System Safety 216, 107923 (2021).
  * Elkarii  _et al._ (2022) M. Elkarii, R. Boukharfane, S. Benjelloun, and C. Bouallou, “A CFD-based surrogate model for predicting slurry pipe flow pressure drops,” Particulate Science and Technology , 1–11 (2022).
  * Ghanem  _et al._ (2017) R. Ghanem, D. Higdon, H. Owhadi, _et al._ ,  _Handbook of uncertainty quantification_ , Vol. 6 (Springer, 2017).
  * Poëtte, Birolleau, and Lucor (2015) G. Poëtte, A. Birolleau, and D. Lucor, “Iterative polynomial approximation adapting to arbitrary probability distribution,” SIAM Journal on Numerical Analysis 53, 1559–1584 (2015).
  * Després, Poëtte, and Lucor (2013) B. Després, G. Poëtte, and D. Lucor, “Robust uncertainty propagation in systems of conservation laws with the entropy closure method,” Uncertainty quantification in computational fluid dynamics , 105–149 (2013).
  * Gilli  _et al._ (2013) L. Gilli, D. Lathouwers, J. L. Kloosterman, T. H. J. J. Van der Hagen, A. J. Koning, and D. Rochman, “Uncertainty quantification for criticality problems using non-intrusive and adaptive polynomial chaos techniques,” Annals of Nuclear Energy 56, 71–80 (2013).
  * Diwan and Mohamed (2019) G. C. Diwan and M. S. Mohamed, “Pollution studies for high order isogeometric analysis and finite element for acoustic problems,” Computer Methods in Applied Mechanics and Engineering 350, 701–718 (2019).
  * Cheng  _et al._ (2019) J. Cheng, J. Wang, X. Wu, and S. Wang, “An improved polynomial-based nonlinear variable importance measure and its application to degradation assessment for high-voltage transformer under imbalance data,” Reliability Engineering & System Safety 185, 175–191 (2019).
  * Vrugt and Robinson (2007) J. A. Vrugt and B. A. Robinson, “Treatment of uncertainty using ensemble methods: Comparison of sequential data assimilation and bayesian model averaging,” Water Resources Research 43 (2007).
  * Geyer (1994) C. J. Geyer, “On the convergence of monte carlo maximum likelihood calculations,” Journal of the Royal Statistical Society: Series B (Methodological) 56, 261–274 (1994).
  * Faivre  _et al._ (2013) R. Faivre, D. Makowski, S. Mahévas, and B. Iooss, “Analyse de sensibilité et exploration de modèles: application aux sciences de la nature et de l’environnement,” Analyse de Sensibilité et Exploration de Modèles , 1–352 (2013).
  * Elkarii  _et al._ (2023) M. Elkarii, R. Boukharfane, S. Benjelloun, C. Bouallou, and N. El Moçayd, “Global sensitivity analysis for phosphate slurry flow in pipelines using generalized polynomial chaos,” Physics of Fluids 35 (2023).
  * Xiu and Karniadakis (2002) D. Xiu and G. E. Karniadakis, “The Wiener–Askey polynomial chaos for stochastic differential equations,” SIAM journal on Scientific Computing 24, 619–644 (2002).
  * Askey and Wilson (1985) R. Askey and J. A. Wilson,  _Some basic hypergeometric orthogonal polynomials that generalize Jacobi polynomials_ , Vol. 319 (American Mathematical Soc., 1985).
  * Ghanem and Spanos (1991) R. G. Ghanem and P. D. Spanos, “Stochastic finite element method: Response statistics,” Stochastic finite elements: a spectral approach , 101–119 (1991).
  * Wang  _et al._ (2019) Y. Wang, K. Hu, L. Ren, and G. Lin, “Optimal observations-based retrieval of topography in 2D shallow water equations using PC-EnKF,” Journal of Computational Physics 382, 43–60 (2019).
  * Efron  _et al._ (2004) B. Efron, T. Hastie, I. Johnstone, and R. Tibshirani, “Least angle regression,” The Annals of Statistics 32, 407–499 (2004).
  * Yamazaki  _et al._ (2018) W. Yamazaki, T. Kato, T. Homma, K. Shimoyama, and S. Obayashi, “Stochastic tsunami inundation flow simulation via polynomial chaos approach,” Journal of Fluid Science and Technology 13, JFST0025–JFST0025 (2018).
  * El Moçayd  _et al._ (2018) N. El Moçayd, S. Ricci, N. Goutal, M. C. Rochoux, S. Boyaval, C. Goeury, D. Lucor, and O. Thual, “Polynomial surrogates for open-channel flows in random steady state,” Environmental Modeling & Assessment 23, 309–331 (2018).
  * Blatman and Sudret (2011a) G. Blatman and B. Sudret, “Principal component analysis and least angle regression in spectral stochastic finite element analysis,” in  _Proc. 11th Int. Conf. on Applications of Stat. and Prob. in Civil Engineering (ICASP11), Zurich, Switzerland_ (2011).
  * Chevreuil (2012) A. Chevreuil, M. Nouy, “Model order reduction based on proper generalized decomposition for the propagation of uncertainties in structural dynamics,” International Journal for Numerical Methods in Engineering 89, 241–268 (2012).
  * Raisee, Kumar, and Lacor (2015) M. Raisee, D. Kumar, and C. Lacor, “A non-intrusive model reduction approach for polynomial chaos expansion using proper orthogonal decomposition,” International Journal for Numerical Methods in Engineering 103, 293–312 (2015).
  * Xiu and Karniadakis (2003) D. Xiu and G. E. Karniadakis, “Modeling uncertainty in flow simulations via generalized polynomial chaos,” Journal of Computational Physics 187, 137–167 (2003).
  * Xiu (2010) D. Xiu,  _Numerical methods for stochastic computations: a spectral method approach_ (Princeton university press, 2010).
  * Blatman and Sudret (2011b) G. Blatman and B. Sudret, “Adaptive sparse polynomial chaos expansion based on least angle regression,” Journal of Computational Physics 230, 2345–2367 (2011b).
  * Choi  _et al._ (2004) S.-K. Choi, R. V. Grandhi, R. A. Canfield, and C. L. Pettit, “Polynomial chaos expansion with latin hypercube sampling for estimating response variability,” AIAA journal 42, 1191–1198 (2004).
  * Berveiller, Sudret, and Lemaire (2006) M. Berveiller, B. Sudret, and M. Lemaire, “Stochastic finite element: a non intrusive approach by regression,” European Journal of Computational Mechanics/Revue Européenne de Mécanique Numérique 15, 81–92 (2006).
  * Iman (2008) R. L. Iman, “Latin hypercube sampling,” Encyclopedia of Quantitative Risk Analysis and Assessment 3 (2008).
  * Blatman (2009) G. Blatman,  _Adaptive sparse polynomial chaos expansions for uncertainty propagation and sensitivity analysis_ , Ph.D. thesis, Clermont-Ferrand 2 (2009).
  * Blatman and Sudret (2010) G. Blatman and B. Sudret, “An adaptive algorithm to build up sparse polynomial chaos expansions for stochastic finite element analysis,” Probabilistic Engineering Mechanics 25, 183–197 (2010).
  * Molinaro, Simon, and Pfeiffer (2005) A. M. Molinaro, R. Simon, and R. M. Pfeiffer, “Prediction error estimation: a comparison of resampling methods,” Bioinformatics 21, 3301–3307 (2005).
  * El Moçayd  _et al._ (2020) N. El Moçayd, M. S. Mohamed, D. Ouazar, and M. Seaid, “Stochastic model reduction for polynomial chaos expansion of acoustic waves using proper orthogonal decomposition,” Reliability Engineering & System Safety 195, 106733 (2020).
  * Rusconi, Lakhouaja, and Kopuz (2016) J. Rusconi, A. Lakhouaja, and M. Kopuz, “The design and engineering of the 187 km khouribga to Jorf Lasfar phosphate slurry pipeline,” Procedia Engineering 138, 142–150 (2016).
  * Weller  _et al._ (1998) H. G. Weller, G. Tabor, H. Jasak, and C. Fureby, “A tensorial approach to computational continuum mechanics using object-oriented techniques,” Computers in physics 12, 620–631 (1998).
  * Weller (2002) H. G. Weller, “Derivation, modelling and solution of the conditionally averaged two-phase flow equations,” Nabla Ltd, No Technical Report TR/HGW 2, 9 (2002).
  * Schiller (1933a) V. L. Schiller, “Uber die grundlegenden berechnungen bei der schwerkraftaufbereitung,” Z. Vereines Deutscher Inge. 77, 318–321 (1933a).
  * Schiller (1933b) L. Schiller, “A drag coefficient correlation,” Zeit. Ver. Deutsch. Ing. 77, 318–320 (1933b).
  * Behzadi, Issa, and Rusche (2004) A. Behzadi, R. I. Issa, and H. Rusche, “Modelling of dispersed bubble and droplet flow at high phase fractions,” Chemical Engineering Science 59, 759–770 (2004).
  * Gidaspow, Bezburuah, and Ding (1991) D. Gidaspow, R. Bezburuah, and J. Ding, “Hydrodynamics of circulating fluidized beds: kinetic theory approach,” Tech. Rep. (Illinois Inst. of Tech., Chicago, IL (United States). Dept. of Chemical …, 1991).
  * Syamlal, Rogers, and O’Brien (1993) M. Syamlal, W. Rogers, and T. J. O’Brien, “MFIX documentation: Volume 1, theory guide. national technical information service, springfield, va,” Tech. Rep. (DOE/METC-9411004, NTIS/DE9400087, 1993).
  * Schaeffer (1987) D. G. Schaeffer, “Instability in the evolution equations describing incompressible granular flow,” Journal of Differential Equations 66, 19–50 (1987).
  * Ding and Gidaspow (1990) J. Ding and D. Gidaspow, “A bubbling fluidization model using kinetic theory of granular flow,” AIChE journal 36, 523–538 (1990).
  * Wilcox (1988) D. C. Wilcox, “Reassessment of the scale-determining equation for advanced turbulence models,” AIAA Journal 26, 1299–1310 (1988).
  * Swamee and Jain (1976) P. K. Swamee and A. K. Jain, “Explicit equations for pipe-flow problems,” Journal of the Hydraulics Division 102, 657–664 (1976).
  * Elkarii, Bouallou, and Ratnani (2020) M. Elkarii, C. Bouallou, and A. Ratnani, “Towards modelling a diphasic flow using the CFD technique to achieve a digital twin of a phosphate slurry piping process,” Chemical Engineering Transactions (2020).
  * Gillies, Shook, and Xu (2004) R. G. Gillies, C. A. Shook, and J. Xu, “Modelling heterogeneous slurry flows at high velocities,” The Canadian Journal of Chemical Engineering 82, 1060–1065 (2004).