SlurryNet: Predicting Critical Velocities and Frictional Pressure Drops in Oilfield Suspension Flows

Abstract

Improving the accuracy of the slurry flow predictions in different operating flow regimes remains a major focus for multiphase flow research, and it is especially targeted at industrial applications such as oil and gas. In this paper we develop a robust integrated method consisting of an artificial neural network (ANN) and support vector regression (SVR) to estimate the critical velocity, the slurry flow regime change, and ultimately, the frictional pressure drop for a solid–liquid slurry flow in a horizontal pipe, covering wide ranges of flow and geometrical parameters. Three distinct datasets were used to develop machine learning models with totals of 100, 325, and 125 data points for critical velocity, and frictional pressure drops for heterogeneous and bed-load regimes respectively. For each dataset, 80% of the data were used for training and the rest 20% for evaluating the out of sample performance. The K-fold technique was used for cross-validation. The prediction results of the developed integrated method showed that it significantly outperforms the widely used existing correlations and models in the literature. Additionally, the proposed integrated method with the average absolute relative error (AARE) of 0.084 outperformed the model developed without regime classification with the AARE of 0.155. The proposed integrated model not only offers reliable predictions over a wide range of operating conditions and different flow regimes for the first time, but also introduces a general framework of how to utilize prior physical knowledge to achieve more reliable performances from machine learning methods.

1. Introduction

This paper addresses the application of machine learning (ML) methods to making accurate and relevant predictions of slurry flow behavior. Slurries are complex multi-phase systems studied actively from a physical perspective for >70 years. Flow regime prediction is inexact, generally relying on semi-empirical correlations that have been fitted to different data sets, which are expensive and non-trivial to gather. These regime predictions are used to make design decisions for pipelines and other transport applications where errors are costly. In applying ML methods in any mature industrial or scientific field one has two choices: (i) start from scratch with no prior knowledge; (ii) incorporate existing knowledge. This second approach is that used here. Thus below, we review both relevant slurry flow fundamentals and ML applications in this domain.

Pipe flows of slurries are commonly encountered in the mining industry (slurry transport) and in oil and gas well operations: hole cleaning, hydraulic fracturing and gravel packing. Here we deal with slurry applications in well drilling, where the relevance of multi-phase flow has long been recognized [1]. Horizontal slurry flow of sand and a Newtonian/non-Newtonian carrier liquid in a pipe geometry is encountered widely in horizontal wells, and also in gathering and transition lines. In drilling engineering, a variety of models are used for cuttings transport. At the lowest level, these simply compare the pumping velocity with a typical particle settling velocity. More sophisticated models [2,3,4] consider flows with a layered structure in the well cross-section, consisting of settled beds with mobile suspension flows above.

These mechanistic models are analogous to those developed earlier to predict slurry transport in the mining industry. Starting from the early work of Durand and Condolios [5], the importance of flow regimes was immediately recognized. Turian and Yuan [6] proposed four different semi-empirical correlations, based on more than 2800 experimental data points, to determine the slurry friction factor in the four flow regimes. This is probably the most comprehensive empirical correlation developed to date. However, an underlying criticism of [6] is that the frictional pressure does not represent the underlying physical balance leading to a solids bed. Transition between the four flow regimes has historically formed one major axis of research on slurry transport. These are typically known as transition velocities, the most important one being the deposition or critical velocity, which defines the onset of a stationary bed at the bottom of the pipe. In applications it is crucial to determine the most efficient pipe size to handle a variety of flow conditions throughout lifetime of the field/pipeline. Accurate prediction of the key flow parameters such as critical velocity, flow regime, and pressure drop has significant impact on such design decisions [7,8].

There are many empirical and semi-empirical models and correlations for predicting the critical velocity, e.g., that of Oroskar and Turian [9] and Kokpinar et al. [10]. One of the significant features of the critical velocity is that it corresponds to the minimum frictional pressure drop in the slurry flow, which has even motivated predictions based on this feature, e.g., see [11]. Likewise, many models and correlations have been developed for prediction of frictional pressure drops, e.g., [12,13] for heterogeneous flow, and [14,15] for the bed-load regime. Additionally, many comprehensive layered models were developed more recently to predict both critical velocity and frictional pressure drop in different regimes, e.g., the two layer model of Gillies et al. [3] and the modified three layer model of Sarraf Shirazi and Frigaard [16].

All the above models are to some extent a combination of phenomenological and mechanistic approaches: once the flow regime is phenomenologically defined this enables a more accurate mechanistic model. However, even targeted mechanistic models are limited by the physical complexity of the actual flows. This suggests that a more data-driven (ML) approach might be effective. In this context, Osman and Aggour [17] developed an artificial neural network for prediction of the frictional pressure drop of a slurry flow in horizontal and near horizontal pipes. The accuracy of their model outperformed that of the existing correlations compared. Ulker and Sorgun [18] used four different machine learning algorithms including k-nearest neighbor (kNN), support vector regression (SVR), linear regression, and ANN to estimate the sedimentation bed height inside a wellbore with and without drill pipe rotation. They found that ANN provided slightly better performance compared to other models. Azamathulla et al. [19] used adaptive neuro-fuzzy interference system (ANFIS) and gene-expression programming (GEP) for prediction of the pressure drop. Their results showed that the ANFIS model led to better performance compared to GEP and existing correlations. Lahiri and Ghanta [20] developed a hybrid SVR and genetic algorithm (GA) technique for prediction of the slurry frictional pressure drop, where GA was used for efficient tuning of SVR hyper-parameters. Their developed model accuracy outperformed that of all the existing correlations.

While the above ML methods have produced positive results for slurry transport in the past 2 decades, the picture is incomplete. First, the estimation of pressure drops only cover the heterogeneous regime. This is a practical drawback: not only are the methods limited to prediction in one regime, but one needs prior knowledge of the flow regime, which is not always the case in practice. Secondly, no dimensional analysis was performed before feeding the parameters as inputs to the algorithm. This necessarily means that there is significant redundancy methodologically. In this study, we address both issues and give a complete model. The key novelty of our approach is that we work with the known physical structure of slurry flows. First we use dimensional analysis to eliminate redundancy in variables. Second we integrate 2 models to mimic the physical studies: (a) a model to predict the regimes and transition; (b) knowing the regime, we predict pressure drop. This improves the accuracy in a physically consistent way.

An outline of the paper is as follows. Below in Section 2 we outline the dimensional analysis and the development of the features as inputs to our models for critical velocity and frictional pressure drop. Section 3 provides a brief background on the ANN and SVR models, and discuss the corresponding important hyper-parameters for each model that need to be tuned for the propose of training. In Section 4 we introduce our modeling and training approach in detail, specially for developing the integrated model for prediction of the slurry friction factor using our knowledge of the flow regime. Section 5 provides the acquired experimental data from the literature, and the detailed results produced by our model with the comparison against the well-known correlations in the literature.

2. Dimensional Analysis and Feature Selection

For a solid–liquid Newtonian slurry flowing through a horizontal pipe, we may assume that the steady flow depends on at least the following parameters: the pipe diameter, $\widehat{D}$, the liquid phase density, ${\widehat{\rho }}_l$, the solids phase density, ${\widehat{\rho }}_s$, the liquid phase viscosity, ${\widehat{\mu }}_l$, the particle diameter in the solids phase, ${\widehat{d}}_p$, gravitational acceleration, $\widehat{g}$, the mean slurry velocity, ${\widehat{U}}_s$, and the mean volumetric concentration of solids in pipe cross section, $C_v$. The last mentioned parameter is dimensionless, whereas the rest are dimensional. Throughout this paper we write all dimensional quantities with a $\widehat{·}$ symbol and dimensionless parameters without.

2.1. Critical Velocity

The deposition velocity, also referred to as the critical velocity, ${\widehat{V}}_c$, is one of the key design parameters for most of the slurry transport systems. It is defined as the velocity, lower than which there exists a stationary bed at the bottom of the pipe. Over the past decades, many researchers have developed empirical and/or semi-empirical correlations and models to predict the critical velocity in pipe geometry.

Table 1. Proposed correlations for the critical velocity.

Researcher	Proposed Correlation
Kokpinar et al. [10] ${}^a$	${\widehat{V}}_c/\sqrt{\widehat{g}\widehat{D}}=0.055{\delta }^{-0.60}{C}_v^{0.27}{(s-1)}^{0.07}R{e}_p^{0.30}$
Durand [21] ${}^b$	${\widehat{V}}_c=F_L\sqrt{2\widehat{g}\widehat{D}(s-1)}$
Zandi et al. [13] ${}^c$	${\widehat{V}}_c={\{\left[40C_v\widehat{D}\widehat{g}(s-1)\right]/\sqrt{C_D}\}}^{0.5}$
Yufin [22] ${}^d$	${\widehat{V}}_c=14.23{\widehat{d}}_p^{0.65}{\widehat{D}}^{0.54}exp\left(1.36{\left[C_v(s-1)\right]}^{0.5}{\widehat{d}}_p^{-0.13}\right)$
Oroskar & Turian [9] e	${\widehat{V}}_c/\sqrt{\widehat{g}{\widehat{d}}_p(s-1)}=1.85C_v^{0.1536}{(1-C_v)}^{0.3564}{\delta }^{-0.378}R{e}_p^{0.09}{x}^{0.30}$

${}^a$$R{e}_p={\widehat{\rho }}_l\widehat{\omega }{\widehat{d}}_p/{\widehat{\mu }}_l$. ${}^b$$F_L$ is a constant. ${}^c$ $C_D$ is the drag coefficient. ${}^d$ Lengths are measured in feet. ^e $x=\frac{2}{\sqrt{\pi }}\{\frac{2}{\sqrt{\pi }}\gamma exp(-4{\gamma }^2/\pi )+{\int }_{\gamma }^{\infty }exp(-4{\gamma }^2/\pi )d\gamma \}$ where $\gamma ={\widehat{V}}_p/{\widehat{V}}_c$.

lists the suggested correlations of Durand [21], Zandi et al. [13], Yufin [22], Oroskar and Turian [9], and Kokpinar et al. [10].

For the prediction of critical velocity, ${\widehat{U}}_s$ is replaced with ${\widehat{V}}_c$ whose value is to be derived. The critical depends on the following parameters:

$${\widehat{V}}_c=f(\widehat{D},{\widehat{\rho }}_l,{\widehat{\rho }}_s,{\widehat{\mu }}_l,{\widehat{d}}_p,\widehat{g},\widehat{\omega },C_v).$$

(1)

Some researchers proposed predictive correlations for obtaining critical velocity in which the particle settling velocity in the mixture, ${\widehat{\omega }}_m$ and the viscosity of the mixture ${\widehat{\mu }}_m$ are involved. However, we know that ${\widehat{\omega }}_m$ is a function of ${\widehat{\mu }}_l$, ${\widehat{\rho }}_l$, $C_v$, and s, and ${\widehat{\mu }}_m$ is a function of ${\widehat{\mu }}_l$ and $C_v$ [19]. By performing dimensional analysis on the parameters in (1) we derive the following dimensionless parameters based on which the critical velocity could be predicted:

$$\frac{{\widehat{V}}_c}{\sqrt{\widehat{g}\widehat{D}}}=f(\delta ,s,R{e}_{pw},C_v),$$

(2)

where $\delta ={\widehat{d}}_p/\widehat{D}$, $s={\widehat{\rho }}_s/{\widehat{\rho }}_l$, and $R{e}_{pw}={\widehat{\rho }}_l\widehat{\omega }{\widehat{d}}_p/{\widehat{\mu }}_l$ are the diameter ratio, density ratio, and the particle Reynolds number respectively. It should be noted that the $R{e}_{pw}$ is based on the settling velocity in clear water $\widehat{\omega }$. The functional relationship (2) among the dimensionless parameters, which has four inputs (features) and one output (target), is used to develop the predictive machine learning algorithms.

2.2. Frictional Pressure Drop

For the prediction of frictional pressure drop ($\frac{d\widehat{p}}{d\widehat{z}}$) of the slurry flow, we can introduce two important dimensionless groups, the Froude number ($Fr$) and Reynolds number ($Re$), which relate the balances of the representative forces and stresses in a slurry pipe flow:

$$\begin{array}{cc}\hfill Re& =\frac{{\widehat{\rho }}_l\widehat{D}{\widehat{U}}_s}{{\widehat{\mu }}_l},\end{array}$$

(3)

$$\begin{array}{cc}\hfill Fr& =\frac{{\widehat{U}}_s^2}{\widehat{g}\widehat{D}(s-1)}.\end{array}$$

(4)

It is possible (and necessary) to define and utilize $Re$ and $Fr$ numbers for prediction of the frictional pressure drop because we have the mean slurry velocity ${\widehat{U}}_s$ as an input parameter here, in contrast to the critical velocity prediction task where this parameter is unknown. We also define the slurry friction factor ($f_{sl}$) as the dimensionless parameter obtained from the frictional pressure drop:

$$f_{sl}=0.5\frac{\widehat{D}\frac{d\widehat{p}}{d\widehat{z}}}{{\widehat{\rho }}_l{\widehat{U}}_s^2}.$$

(5)

Therefore, the dimensionless parameters governing the slurry friction factor are as follows:

$$f_{sl}=f(Re,Fr,s,C_v,\delta ).$$

(6)

It has also been found that the friction factor $f_w$ of the clean water with the same flow parameters is useful for prediction of the slurry flow friction factor. $f_w$ can be derived by the Colebrook-White correlation which gives a Darcy-Weisbach friction factor as a function of the Reynolds number and the roughness of the pipe ($f_w=f_{CW}(Re,{ϵ}_r)$). Therefore, we also add $f_w$ as an extra feature which potentially helps the model performance in prediction. As could be observed, the functional relationship (6) among the dimensionless parameters, plus $f_w$ has six features based on which the target $f_{sl}$ should be predicted.

3. Machine Learning Methodology

We now briefly outline the background to the methods we have used in this study. For more detail on these methods the reader is referred to [23,24] for SVR modeling and [25,26] for ANN.

3.1. ANN Modeling

Artificial neural networks (ANNs) are composed of a large number of interconnected processing elements called neurons or cells, that are tied together with weighted connections. Neural networks are inspired by the system of biological neurons whose connections are provided by synapses [25,27]. Learning process in neural networks occurs in a similar way through training and provision of the true input and output dataset, where the connection weights are iteratively being adjusted to solve the specific problem at hand.

The most widely applied feed forward ANN for supervised regression and classification purposes is multi-layer perceptron (MLP), which consists of an input layer, an output layer, and one or more hidden layer(s) in which each layer has a weight matrix, W and a bias vector, b [26]. Figure 1 illustrates the architecture of an MLP network. Observe that each node in every layer of MLP, including the bias node, is fully connected to all the nodes in the subsequent layer. The number of nodes in the input layer is equal to the number of input parameters. The output layer may also contain more than one nodes, corresponding to the number of predictions the network is responsible for making. However, the number of hidden layers and the number of their nodes are adjustable hyperparamters so that the model satisfies the desired approximation and suitable generalization capability.

Considering the feed-forward process of a single data point, for obtaining the values of all the $n^{\left[l\right]}$ nodes in layer l we first calculate the vector $z^{\left[l\right]}$($n^{\left[l\right]},1$) which is a linear function of the values of nodes in the previous layer, $a^{[l-1]}$($n^{[l-1]},1$), i.e., $z^{\left[l\right]}={\sum }_{j=1}^{n^{[l-1]}}{w}_{ij}^{\left[l\right]}{a}_j^{[l-1]}+b_i^{\left[l\right]}$. Subsequently, a nonlinear activation function, $\psi \left(z\right)$, is applied element-wise to the vector $z^{\left[l\right]}$ to get the final values of all the nodes in layer l contained in the vector $a^{\left[l\right]}=\psi \left(z^{\left[l\right]}\right)$. If we have m data points in the training batch, we can write the feed-forward equations in matrix formation as:

$$\begin{array}{cc}\hfill & Z^{\left[l\right]}=W^{\left[l\right]}{A}^{[l-1]}+b^{\left[l\right]},\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill & A^{\left[l\right]}=\psi \left(Z^{\left[l\right]}\right),\hfill \end{array}$$

(8)

where the matrix $A^{\left[l\right]}$($n^{\left[l\right]},m$) contains the obtained values of the nodes in layer l, and $W^{\left[l\right]}$($n^{\left[l\right]},n^{[l-1]}$) and $b^{\left[l\right]}$($n^{\left[l\right]},1$) are the adjustable parameter (weight) matrix and the bias vector respectively.

The activation function acts as a mathematical gate between the inputs that are fed to a neuron and its output that is going to the next layer. Non-linear activation functions allow the model to create complex mappings between each layer’s input and output, which are vital for learning and modeling complex data [28]. Whereas, using linear activation functions leads to a model as simple as linear regression with a significant under-fitting problem. An additional important aspect of activation functions is that they should be computationally inexpensive. The most common activation functions for MLPs are Logistic sigmoid function, $\psi \left(z\right)=\frac{1}{1+e^{-z}}$, Hyperbolic tangent function, $\psi \left(z\right)=\frac{e^z-e^{-z}}{e^z+e^{-z}}$, ReLU function, $\psi \left(z\right)=max(0,z)$, and Leaky ReLU function $\psi \left(z\right)=max(0.01z,z)$.

The training is an iterative process during which the network tries to “learn” the relationship between the provided input(s) and the corresponding output(s) by altering the weights and biases to achieve a satisfactory prediction within a reasonable error margin. The weights and biases are slightly adjusted during each iteration through the training set until the mentioned task is accomplished. At this stage, the learning process is finished on the training set, and the model is ready to be examined on the unseen data. Iterating once over the entire training set is called one epoch. The optimized number of epochs for the training purpose depend on many hyper-parameters such as the learning rate, the optimization algorithm, complexity of the dataset itself, etc.

The back propagation algorithm is the most popular method for modifying and adjusting the weights and biases in every iteration, in which the difference (error) in the ground truth and the obtained outputs is propagated back to each layer and the weights and biases get adjusted accordingly [29]. The main goal of the back propagation process is that the prespecified loss function is minimized. The most widely used loss function for training a regression problem is the mean squared error:

$$L(y,y^{\prime })=\frac{1}{m}\sum _{i=1}^m{(y_i^{\prime }-y_i)}^2,$$

(9)

where $y^{\prime }$ and y are the predicted and true outputs respectively.

All machine learning models including neural networks are considered to be satisfactory in terms of their prediction(s) in case they perform well on the unseen dataset (test set) which was not used in the training process. In other words, the learning model has a suitable generalization capability if the out of sample error is within an acceptable margin. It is possible that a learning model performs well on the training set, but fails to make accurate predictions on the test set. this issue is referred to as overfitting or variance problem, in which the generalizability of the model is poor. On the other hand, if the model or MLP structure is too simple so that the performance on both the training and test sets are poor, the model has high bias or underfitting problem. To examine a specific trained model, one should perform a cross-validation procedure in which the generalizability of the model is monitored on the data that was not seen by the network during the training. The most suitable model is the one with the lowest cross-validation error. This can be accomplished either via defining another set (development set) or by performing K-fold cross validation [30].

Finding the most suitable model for a supervised learning task is essentially a trial and error process which often involves conducting a grid search over different hyperparameters of the model. For instance, MLP hyperparameters include the learning rate, weight matrices initialization, optimization algorithm, network architecture, regularization parameter, number of training epochs, etc.

3.2. SVR Modeling

Support vector regression (SVR) is the most common application form of support vector machines (SVMs). Suppose our training set contains $\{(x_1,y_1),(x_2,y_2),\dots ,(x_m,y_m)\}$ where $x_i\in {\mathbb{R}}^n$ are input variables (features) and $y_i\in \mathbb{R}$ are the corresponding output (target) value. The modeling goal in $ϵ$-SVR is to propose a function $f\left(x\right)$ that has a maximum of $ϵ$ deviation from the actual target values for all of the input variables in the training set [24,31]. Meanwhile, the obtained function should be as flat as possible to prevent the high variance issue. For this purpose, the loss function is penalized only in case the predicted output deviation from the target is more than $ϵ$. SVR considers the following linear estimation function:

$$f\left(x\right)=\langle w,x\rangle +b,$$

(10)

where $w\in {\mathbb{R}}^n$ and $b\in \mathbb{R}$ denote the weight vector and bias respectively, and $\langle .,.\rangle$ represents the dot product in the feature space ${\mathbb{R}}^n$. A viable solution for increasing the flatness of $f\left(x\right)$ is to minimize the norm of w. Therefore, we transform to the following convex optimization problem:

$$\begin{array}{c}\hfill \begin{array}{cccc}\mathrm{minimize}\hfill & \hfill & & \frac{1}{2}{\parallel w\parallel }^2\hfill \\ \mathrm{subject}\mathrm{to}\hfill & \hfill & & \left\{\begin{array}{c}{\displaystyle y_i-\langle w,x_i\rangle -b\le ϵ}\hfill \\ {\displaystyle \langle w,x_i\rangle +b-y_i\le ϵ}\hfill \end{array}\right.\hfill \end{array}\end{array}$$

(11)

The optimization problem (11) assumes there exists a function f such that the prediction errors for all the data points are within the $ϵ$ margin. However, this assumption might not be always satisfied. One might also consider allowing for some errors by defining a “soft margin” loss function which was primarily adapted to the SVMs by Cortes and Vapnik [32] to prevent overfitting and enhance the generalizability of the proposed function/model. To fulfill this purpose, slack variables $\xi$ and ${\xi }^{*}$ can be introduced to modify the infeasible constraints of the optimization problem (11). The corresponding formulation is as follows:

$$\begin{array}{c}\hfill \begin{array}{cccc}\mathrm{minimize}\hfill & \hfill & & \frac{1}{2}{\parallel w\parallel }^2+C\sum _{i=1}^m({\xi }_i+{\xi }_i^{*})\hfill \\ \mathrm{subject}\mathrm{to}\hfill & \hfill & & \left\{\begin{array}{c}{\displaystyle y_i-\langle w,x_i\rangle -b\le ϵ+{\xi }_i}\hfill \\ {\displaystyle \langle w,x_i\rangle +b-y_i\le ϵ+{\xi }_i^{*}}\hfill \\ {\displaystyle \xi ,{\xi }^{*}\ge 0}\hfill \end{array}\right.\hfill \end{array}\end{array}$$

(12)

where the constant $C>0$ determines the penalizing degree of the output predictions with the errors larger than $ϵ$, which also regulates the flatness of f. Therefore, the corresponding $ϵ$-intensive loss function ${\left|\xi \right|}_{ϵ}$ is described by:

$$\begin{array}{c}\hfill {\left|\xi \right|}_{ϵ}=\left\{\begin{array}{c}{\displaystyle 0\mathrm{for}\left|\xi \right|\le ϵ}\hfill \\ {\displaystyle \left|\xi \right|-ϵotherwise.}\hfill \end{array}\right.\end{array}$$

(13)

It is found that the optimization problem (12) can be solved more simply when converted to its dual formation. For this purpose, Lagrange multipliers can be used as described in Fletcher et al. where the idea is to construct a Lagrange function from the primal objective function and the corresponding constraints, by introducing a dual set of variables. Considering the Lagrangian function and its properties, by performing some mathematical operations (see details in [24]), we arrive at the following dual optimization problem:

$$\begin{array}{c}\hfill \begin{array}{cccc}\mathrm{maximize}\hfill & \hfill & & \left[-\frac{1}{2}\sum _{i,j=1}^m({\alpha }_i-{\alpha }_i^{*})({\alpha }_j-{\alpha }_j^{*})\langle x_i,x_j\rangle -ϵ\sum _{i=1}^m({\alpha }_i+{\alpha }_i^{*})+\sum _{i=1}^m{y}_i({\alpha }_i-{\alpha }_i^{*})\right]\hfill \\ \mathrm{subject}\mathrm{to}\hfill & \hfill & & \sum _{i=1}^m({\alpha }_i-{\alpha }_i^{*})=0\mathrm{and}{\alpha }_i,{\alpha }_i^{*}\in [0,C],\hfill \end{array}\end{array}$$

(14)

where ${\alpha }_i$, ${\alpha }_i^{*}$ are Lagrange multipliers. In practice, only some of the coefficients $({\alpha }_i-{\alpha }_i^{*})$ are non-zero due to the specific character of the quadratic programming problem (14). The corresponding input vectors $x_i$ whose coefficients are non-zero are referred to as support vectors (SVs). The SVs could be considered as the data points that represent the information content of the entire training dataset. The final form of the estimation function using (14):

$$f\left(x\right)=\sum _{i=1}^m({\alpha }_i-{\alpha }_i^{*})\langle x_i,x\rangle +b.$$

(15)

From Equation (15) it can be deduced that $w={\sum }_{i=1}^m({\alpha }_i-{\alpha }_i^{*})x_i$. Figure 2 demonstrates a schematic of the diagram of a non-linear support vector regression using $ϵ$-sensitive loss function. As could be observed, according to (13) if the predicted value of a data point (blue dots) are within the $ϵ$-tube the loss is zero; Otherwise (red dots), the loss is equal to the magnitude of the difference between the predicted value and the radius of the tube $ϵ$. As can be observed in Figure 2 the SVR algorithm tries to situate the tube around the data points with the help of the support vectors (green dots).

The true power of SVMs can be achieved by introducing non-linearity to the original algorithm. This could be accomplished by mapping the training data points to another feature space $\Phi :{\mathbb{R}}^n\to {\mathbb{R}}^k$ with higher dimension $k>n$, in which the dot product operation can be substituted by a kernel function, i.e., $K(x_i,x_j)=\varphi \left(x_i\right)\varphi \left(x_j\right)$. In other words, the kernel function represents the dot product in a higher dimensional feature space. Substituting the kernel function in Equation (15) introduces the mentioned non-linearity to the SVM algorithm:

$$f\left(x\right)=\sum _{i=1}^m({\alpha }_i-{\alpha }_i^{*})K(x_i,x)+b.$$

(16)

It is important to note that kernel functions ought to have some key specific characteristics so that they correspond to a dot product operation in some other feature space. The two most widely used kernel functions in SVM algorithm are as follows:

$$\begin{array}{cccc}\hfill \mathrm{Polynomial}\mathrm{Function}& :\hfill & \hfill & K={(u^{T}v+1)}^P,\hfill \end{array}$$

(17)

$$\begin{array}{cccc}\hfill \mathrm{Gaussian}\mathrm{Radial}\mathrm{Basis}\mathrm{Function}& :\hfill & \hfill & K=exp\left(-\frac{{\parallel u-v\parallel }^2}{2{\sigma }^2}\right),\hfill \end{array}$$

(18)

where u, v are kernel arguments, P is degree of the polynomial, and $\sigma$ is the width of radial basis function (RBF).

Similar to ANNs it is important to come up with an optimized set of hyper-parameters for SVR algorithm so that the proposed hypothesis function f offers an acceptable generalization performance and lacks the high bias or variance problem. The tunable hyper-parameters of SVR are $ϵ$, C, the kernel type, and the corresponding kernel parameters, e.g., P for the polynomial and $\gamma =\frac{1}{2{\sigma }^2}$ for the RBF kernels. Choosing a specific kernel type is usually based on the application domain knowledge and is also required to reflect the distribution of the dataset. C determines the trade-off between the flatness of the hypothesis function and the degree up to which deviations larger than $ϵ$ are tolerated. In other words, it also has a regularization effect, such that the smaller the value of C, the more significantly the objective function in (12) is regularized. The $ϵ$ parameter defines the radius of a “tube” zone in which the loss function is zero: larger $ϵ$ selection leads to a proposed hypothesis function with more flatness and also fewer number of support vectors. Therefore, it can be deduced that both C and $ϵ$ affect the model complexity.

4. Modeling Approach

The purpose of this study is to develop learning models using ANN and SVR algorithms for prediction of the critical velocity and frictional pressure drop of slurry flow in pipe geometry. For critical velocity, we use the four dimensionless features, $\delta$, s, $R{e}_{pw}$, and $C_v$ developed in Section 2.1 as inputs to develop the above mentioned learning algorithms with satisfactory generalizability. However, for the prediction of frictional pressure drop we also need to understand the effect of the slurry flow regime on the friction factor.

Figure 3 shows a schematic of the frictional pressure drop as a function of the mean slurry velocity for different flow regimes. The slurry flow regime is governed by the competition between the turbulent eddies and the particle settling tendency due to gravity. The former tends to suspend the solid particles in the carrier liquid while the latter drives the particles to settle at the bottom of the pipe. The frictional pressure drop of a slurry flow depends on the different existing stresses and forces whose nature and strength strongly depend on the flow regime [6,16].

At low flow rates, the turbulent eddies are not strong enough to suspend the solid phase. As a result, a considerable portion of the pipe is occupied by stationary sedimentation bed above which there is a heterogeneous layer with a recognizable solids concentration gradients. This regime of the slurry flow is also referred to as the bed-load regime. As observed in Figure 3 the frictional pressure drop decreases with the mean velocity in this regime. This is explained by the fact that at low velocities, the stresses and forces are dominated by the solids phase, that are weakened as the slurry velocity increases.

As the mean superficial velocity increases, the turbulent eddies become more capable of suspending the solids until all the static bed layer is eroded, and there is a moving bed layer at the bottom of the pipe whose concentration is close to maximal packing. As the flow rate is further increased, we reach the heterogeneous or fully suspended regime where there is a solid concentration gradient in the direction of gravity. At extremely high flow rates, turbulent eddies become significantly more dominant and the solid phase becomes progressively more homogeneously distributed in the carrier liquid. As shown in Figure 3 the frictional pressure drop increases with the mean velocity through saltation flow, heterogeneous, and homogeneous regimes. Furthermore, the pressure drop increase rate is also increasing at higher velocities, as the liquid phase role becomes more dominant in the suspension stresses.

As noted above, the frictional pressure drop behavior is noticeably affected when the regime changes from the bed-load to saltation flow, i.e., at the critical velocity. Therefore, we can introduce this prior knowledge to our predictive modeling approach. Figure 4a,b shows the work flow chart for developing our predictive models. We develop two separate learning models with satisfactory accuracy and generalization capability for the bed-load and heterogeneous flow regimes according to the work-flow chart illustrated in Figure 4a. For this task, we also need to train the two mentioned models with separate datasets representing the corresponding regimes. For examining the generalizability of the developed predictive model for frictional pressure drop, we primarily check what the flow regime is, based on the developed model for critical velocity. Subsequently, we feed the six dimensionless parameters (see Section 2.2) as features to the corresponding predictive learning model for frictional pressure drop prediction. This procedure is illustrated in Figure 4b for further clarification of our integrated method scheme for prediction of the slurry friction factor. Consequently we have a dataset for critical velocity, and two distinct datasets for frictional pressure drop: in bed-load regime and the rest of the regimes.

We develop the most suitable ANN and SVR predictive models for each of the three datasets via grid search among their corresponding hyperparameters. The chosen hyperparameters of ANN for tuning include the architecture of the network, i.e., the number of hidden layer(s) and neurons in each hidden layer, activation function, number of training epochs, and learning rate, and the ones for SVR include C, $ϵ$, kernel type, and kernel parameter (polynomial degree for polynomial function, and $\gamma =\frac{1}{2{\sigma }^2}$ for the radial basis function). Then we pick the one with the best validation score as our ultimate proposed model. For the purpose of model development we take 80% of each dataset randomly as the training set and the remaining 20% as the test set. We perform 5-fold cross validation on the training set to examine the generalization capacity of the model on the data that it did not get trained on. The best model with specific sets of hyperparameters is chosen based on this validation score.

5. Results and Discussion

As the magnitude of the input features are significantly different, the data should be normalized before being fed to the training algorithms. If the inputs are of different scales, the weights connected to the inputs with larger scales will be updated much faster compared to others, which can considerably hurt the learning process. On the other hand, there are also a variety of practical reasons why normalizing the inputs can make training faster and reduce the chances of getting stuck in local optima. We use the standard normalization as follows:

$$x_{norm,i}=\frac{x_i-u_{train}}{{\sigma }_{train}},$$

(19)

where $x_{norm,i}$ is the normalized input of the $i$ sample, and $u_{train}$ and ${\sigma }_{train}$ are the mean and standard deviation of the data points in training set. The output is also normalized in similar way as in (19).

Table 2. Parameters of the experimental data considered for comparison with critical velocity.

Source	Data Sets	$\widehat{\mathit{D}}$$\left[\mathbf{m}\right]$	${\widehat{\mathit{d}}}_{\mathit{p}}$$\left[\mathbf{mm}\right]$	${\widehat{\mathit{V}}}_{\mathit{c}}$$\left[\mathbf{m}/\mathbf{s}\right]$	${\mathit{C}}_{\mathit{v}}$	s
Kokpinar et al. [10]	42	0.15	1.09–5.34	1.06–3.00	0.011–0.091	1.04–2.6
Graf et al. [33]	12	0.102; 0.152	0.45–0.88	1.55–2.42	0.007–0.07	2.65
Durand [12]	7	0.15	0.44–2.04	2.19–2.71	0.05–0.15	2.6
Avci [34]	15	0.052	0.29–3.2	0.27–1.58	0.05–0.30	1.04–2.68
Yotsukura [35]	11	0.108	0.23–1.15	1.83–2.96	0.05–0.25	2.6
Wicks [36]	2	0.027; 0.14	0.25	0.46–0.79	0.01	2.6
Sinclair [37]	11	0.025	2.205	0.32–0.52	0.03–0.18	1.74

shows parameters of 100 experimental data points collected from the literature, measuring critical velocity, which we use to train and test our proposed models. Additionally, Figure 5a–e demonstrate the estimation of the probability density function and box plot of all the input features along with the output, which provides insightful information about the distribution and statistical parameters of dataset. Each data point is the result of an experimental test by the listed authors, performed in different flow loop facilities. As can be observed, these experiments cover a wide range of particle sizes ${\widehat{d}}_p=0.23\mathrm{mm}-5.34\mathrm{mm}$, pipe diameters $\widehat{D}=0.025\mathrm{m}-0.152\mathrm{m}$, mean solids concentrations $C_v=0.007-0.30$, and also different density ratios $s=1.04-2.68$. Most of the data are taken from the measurements conducted by Kokpinar et al. [10] who used coarse particles, with different materials to also see the effect of s on the critical velocity. They used sand, coarse sand, coal, blue plastic, black plastic, fine tuff, and coarse tuff with specific densities of $s=2.60,2.55,1.74,1.20,1.35,1.31,1.04$ respectively.

We train and obtain a validation score (loss) on the randomly chosen 80 data points (training set) and report the out of sample results on the remaining 20 data points.

Table 3. Optimum hyper-parameters obtained by the SVR algorithm.

Case	C	$\mathit{ϵ}$	Kernel Type	Kernel Parameter	Validation Loss
Critical Velocity	40	0.05	RBF	0.2	0.059
Heterogeneous Friction Factor	800	0.025	RBF	0.1	0.097
Bed-load Friction Factor	10	0.05	RBF	0.6	0.123
Batch-trained Friction Factor	50	0.025	RBF	0.5	0.149

shows the the optimum hyperparameters for SVR algorithm,

Table 4. Parameters of the experimental data considered for comparison for frictional pressure drop.

Source	Regime	$\widehat{\mathit{D}}$$\left[\mathbf{m}\right]$	${\widehat{\mathit{d}}}_{\mathit{p}}$$\left[\mathbf{mm}\right]$	${\widehat{\mathit{U}}}_{\mathit{s}}$$[\mathbf{m}/\mathbf{s}]$	${\mathit{C}}_{\mathit{s}}$	s	$-\frac{\mathit{d}\widehat{\mathit{p}}}{\mathit{d}\widehat{\mathit{z}}}$$[\mathbf{kPa}/\mathbf{m}]$
Gillies et al. [3]	Heterogeneous	0.103	0.09; 0.27	1.49–7.77	0.10–0.40	2.65	0.37–5.32
Schaan et al. [39]	Heterogeneous	0.053	0.085–0.1	0.99–5.02	0.15–0.40	2.44–2.66	0.27–7.20
Matousek [40]	Heterogeneous	0.155	0.37	4.72–8.98	0.12; 0.26	2.65	0.99–3.53
Doron et al. [38]	Heterogeneous	0.051	3.00	0.55–1.63	0.042–0.115	1.24	0.22–0.63
Durand [5]	Bed-load	0.15	0.44; 2.04	1.10–2.13	0.085–0.26	2.65	0.52–2.13
Clift [41]	Bed-load	0.44	0.29–0.68	1.73–3.81	0.11; 0.15	2.65; 2.87	0.36–1.05
Yagi [42]	Bed-load	0.08–0.15	0.25–1.28	1.00–2.81	0.15	2.63–2.67	1.01–4.67
Doron et al. [38]	Bed-load	0.051	3.00	0.24–0.55	0.042–0.115	1.24	0.19–0.41

shows parameters of experimental data points collected from the literature, and

Table 5. Optimum hyper-parameters obtained by the ANN algorithm.

Case	Hidden Layers	Neurons	Activation Function	Epochs	Learning Rate	Validation Loss
Critical Velocity	1	16	Leaky ReLU	120	0.08	0.072
Heterogeneous Friction Factor	2	16	Leaky ReLU	350	0.02	0.090
Bed-load Friction Factor	2	14	Leaky ReLU	500	0.02	0.112
Batch-trained Friction Factor	2	18	Leaky ReLU	700	0.01	0.155

shows the optimum hyperparameters for ANN algorithms. Table 3 and Table 5 show their corresponding validation loss respectively. It should be noted that the validation loss refers to the average mean squared error obtained by the 5-fold cross validation. As observed, the optimum SVR model outperforms ANN in terms of the validation score and hence the generalization capability. Therefore, the SVR model is chosen as the ultimate prediction model for critical velocity.

Table 6. Performance of the chosen models on training and test sets.

Case	Chosen Model	Set	$\mathit{AARE}$	$\mathit{\sigma }$	R
Critical Velocity	SVR	training	0.073	0.153	0.959
		test	0.099	0.207	0.920
Heterogeneous Friction Factor	ANN	training	0.017	0.013	0.997
		test	0.026	0.034	0.992
Bed-load Friction Factor	ANN	training	0.025	0.024	0.999
		test	0.054	0.085	0.997
Batch-trained Friction Factor	SVR	training	0.097	0.096	0.963
		test	0.155	0.178	0.926
Integrated method Friction Factor	SVR-ANN	test	0.084	0.215	0.991

shows the performance of the chosen model on training and test sets, in terms of the average absolute relative error (AARE), the cross-correlation coefficient (R), and the standard deviation of error ($\sigma$). We can compare the proposed model performance against the most widely used predictive correlations in literature brought in Table 1. The out of sample average absolute relative error of these models are 0.099, 0.153, 0.308, 0.322, 0.412, and 0.447 corresponding to the prediction of the proposed SVR model, Kokpinar et al. [10], Oroskar and Turian [9], Durand [11], Yufin [22], and Zandi et al. [13] respectively. It is evident that the prediction error of critical velocity has reduced considerably in the present work. Figure 6 shows the parity plot of the experimentally measured and predicted results of the dimensionless critical velocity for the training and test sets with the AARE of 0.073 and 0.099 respectively.

We have also directly compared the performance of our model with that of Kokpinar et. al. [10] with their own 42 experimental data points. Figure 7 shows the parity plot of the corresponding predictions versus the measured dimensionless critical velocities. The AARE of estimations are 0.142 and 0.062 for Kokpinar et al. [10] and present models respectively. As observed in Figure 7 the present model performs better in particular where $\frac{{\widehat{V}}_c}{\sqrt{\widehat{g}\widehat{D}}}>1.8$.

Figure 8 illustrates the effect of hyperparameter C on the loss function (mean squared error) of the training, validation, and test sets. As was mentioned in Section 3.2, C determines the trade-off between the flatness of the hypothesis function, and the degree up to which deviations larger than $ϵ$ are tolerated in SVR algorithm. In practice it also has regularization effect such that the lower the value of C, the more the objective function is regularized. As seen in Figure 8 there is an optimal C where the loss function is minimized in validation and test sets, lower than which the hypothesis function suffers from high bias (under-fitting) and higher than which it suffers from high variance (over-fitting). Obviously, the values of all hyperparameters including C are chosen based on the validation score.

Table 4 shows parameters of experimental data points collected from the literature, measuring frictional pressure drop for heterogeneous and bed-load regimes, which we use to train and test our proposed models. The total number of experimental data points are 365 and 125 for the heterogeneous and bed-load regimes respectively. As can be observed, the experiments mostly used fine particles except for the Doron et al.’s data [38] and part of Durand’s measurements in bed-load regime [5], where particle sizes of ${\widehat{d}}_p=0.23$ mm and $5.34$ mm were used respectively. Pipe diameters of the range $\widehat{D}=0.051\mathrm{m}-0.155\mathrm{m}$ were used in the experiments with flow velocity range of ${\widehat{U}}_s=0.24\mathrm{m}/\mathrm{s}-7.77\mathrm{m}/\mathrm{s}$, and mean delivered solids concentration of $C_s=0.042-0.40$. Most of the experiments were conducted using sand particles with the density ratio of $s=2.44-2.87$ except for Doron et al.’s work where General Electric “Black Acetal” with the density ratio of $s=1.24$ was used [38]. Figure 9a–g show the kernel density estimation and box plot of all the input features along with the output for both heterogeneous and bed-load regime datasets. An interesting observation is that the distribution of $Fr$ and $f_{sl}$ are considerably different comparing the two regimes. The reason is that according to (4) and (5) both of these dimensionless variables include the term ${\widehat{U}}_s^2$ in their equations, and we know that the mean slurry velocity in the bed-load regime is less than that of the heterogeneous regime. Therefore, $Fr$ is considerably lower while $f_{sl}$ is larger in bed-load regime compared to the heterogeneous regime.

Similar to the critical velocity case, we randomly take 80% of the dataset for the purpose of training and validation, and the rest 20% as the test set for evaluating the out of sample performance. As can be observed from Table 3 and Table 5 the most efficient developed ANN models are outperforming SVR for both heterogeneous and bed-load regimes. Figure 10a,b show the parity plot comparing the measured and predicted slurry friction factor for both regimes. The corresponding out of sample results are illustrated in Table 6.

For a fair comparison against the existing correlations and models from literature, we also need to investigate the integrated method performance in terms of predicting the frictional pressure drop. In other words, we would like to determine the out of sample error where we ignore the prior knowledge of the flow regime, which can be the case in real-life scenarios specifically for industrial applications. To serve this purpose, we feed each data point to the developed SVR algorithm for critical velocity prediction, and compare the predicted result with the mean slurry velocity as a means to identify the regime. For this process the key assumption is $C_v=C_s$ at the critical velocity which is a reasonable assumption to make. After the regime identification, we feed the data point to the corresponding model for predicting the frictional pressure drop. The out of sample results for integrated method prediction is shown in Table 6.

Once again, we can compare the out of sample AARE against that of some recognized correlations and models available in literature for predicting the pressure drop. For slurry friction factor prediction in heterogeneous regime, the AARE of the correlations developed by Zandi and Govatos [13], Durand and Condolios [5], and Turian and Yuan [6] are 0.643, 0.449, and 0.348 respectively, whereas for the bed-load regime, the AARE of the proposed models by Gruesbeck et al. [14], Penberthy et al. [15], and Turian and Yuan [6] are 0.837, 0.769, and 0.529 respectively. It is clear that the prediction performance of the current study with AARE of 0.084 significantly outperforms that of the mentioned models.

Figure 11 illustrates the effect of the epoch number, a key hyperparameter for ANNs, on the loss function of the training, validation, and test sets for the heterogeneous regime ANN model. As can be observed, there is an optimal epoch number for training, after which the validation loss starts to increase. In other words, after around 400 training epochs the model is over-fitting on the training dataset.

To investigate whether the proposed integrated method is indeed required for the purpose of a satisfactory prediction for frictional pressure drop, we have also performed a batch train using all of the 490 frictional pressure drop data points, without any supervised or unsupervised classification based on the flow regime. We trained and tested another learning model under the mentioned condition with the similar procedure as other developed models. Table 3 and Table 5 show that the SVR model performance is more satisfactory compared to ANN in terms of generalization capacity. Figure 12 illustrates the corresponding parity plot for the measured slurry friction factor against the predicted values.

For comprehending and comparing the performance of the batch-trained model with the integrated method, the corresponding parity plots are shown in Figure 13a,b. Figure 13a shows the measured and predicted slurry friction factor for the integrated method. As could be observed, there are four heterogeneous data points whose regime was incorrectly classified as bed-load (blue squares), and three bed-load data points with false classification. As shown, the predicted slurry friction factor for misclassified heterogeneous data points tend to be higher than the measured value, whereas the reverse is true for the misclassified bed-load data points. The reason is that generally, the value of slurry friction factor in bed-load regime is more than that of the heterogeneous regime. However, the out of sample results of the integrated method is more satisfactory compared to the batch-trained model, with the AARE of 0.084 for the former and 0.155 for the latter, as shown in Table 6. Consequently, it can be comprehended that although the integrated method prediction highly relies on the performance of regime classification, i.e., the SVR model for critical velocity prediction, it is considered to be more efficient in practice to predefine a regime classification method, such as the one accomplished in this work, prior to feeding it to the model for satisfactory prediction of the friction factor.

6. Summary

We have developed a robust integrated method using ANN and SVR algorithms for prediction of critical velocity and frictional pressure drop by identifying and implementing existing knowledge of the main slurry regimes. The proposed model clearly outperforms the estimation of existing well-known and widely used correlations and models for the prediction of critical velocity and frictional pressure drop. Furthermore, it overcomes the limitations of previous machine learning models which only targeted the estimation of frictional pressure drop in the heterogeneous regime.

The features have been extracted based on dimensional analysis of geometrical and flow parameters that are involved in the governing equations of a solid liquid slurry flow in pipe. This ensures that we preserve all the data information with the least input dimension, which is one of the main goals in developing machine learning algorithms and other methods for prediction. Indeed this is a relatively simple step that can be taken for any physical/mechanical scenario. Additionally, we have shown that the slurry friction factor estimation noticeably improves with regime classification before feeding the data to the developed model.

One of the limitations of the proposed integrated method is that its accuracy highly relies on the regime’s classification performance. However, the overall prediction accuracy can be improved by ensuring that the data used for training the critical velocity and frictional pressure drop models are provided from the same distribution. Another limitation of this study is the limited number of data points available for efficient training of the proposed models. Additionally, using more complex machine learning methods along with more available data can be considered as a suitable future task.

In general, the message of the paper is that one should not discard old methodologies in assuming that new machine learning algorithms will automatically solve all problems. The challenge in industrial applications where we need to predict important variables, is to integrate new predictive methodologies with the old and with our prior physical knowledge/know how. In this respect our results are promising in showing significant advance in predictive abilities with a small investment in dimensional analysis.