Flow status identification based on multiple slow feature analysis of gas–liquid two-phase flow in horizontal pipes

In the SFA method, slowness is regarded as a measure of how fast a signal fluctuates. The definition of slowness, given a signal x(k), is:

$$\begin{align}\Delta (x) = {\langle \dot x{(k)^2}\rangle _k}\end{align} \tag{ 1 }$$

where $\dot x(k) = x(k) - x(k - 1)$ denotes the difference which approximates the first-order derivative of x, and ${\langle \cdot \rangle _k}$ denotes the time average calculated using N data samples.

Assuming a J-dimensional input time series signal ${\mathbf{x}}(k) = {[{x_1}(k),{x_2}(k), \ldots ,{x_J}(k)]^T}$, the objective of SFA is to find a set of projection directions that make the output signal change as slowly as possible

$$\begin{align}\mathop {\min }\limits_{{g_j}( \cdot )} \Delta ({s_j}): = {\langle \dot s_j^2\rangle _k}\end{align} \tag{ 2 }$$

under the constraints

$$\begin{align}\left\{ \begin{array}{l} \displaystyle {\langle {s_j}\rangle _k} = 0{\text{ }}({\text{zero mean}}) \\ \displaystyle {\langle s_j^2\rangle _k} = 1{\text{ }}({\text{unit variance}}) \\ \displaystyle \forall i \ne j:{\langle {s_i}{s_j}\rangle _k} = 0{\text{ }}({\text{decorrelation and order}}) \\ \end{array} \right..\end{align} \tag{ 3 }$$

Each s_j (k) combines input components into temporally slow features in a linear way and w _j is the direction of projection:

$$\begin{align}{s_j}(k) = {\mathbf{w}}_j^T{\mathbf{x}}(k).\end{align} \tag{ 4 }$$

The mapping from x(k) to s(k) can be written as

$$\begin{align}{\mathbf{s}}(k) = {\mathbf{Wx}}(k)\end{align} \tag{ 5 }$$

where ${\mathbf{W}} = {[{{\mathbf{w}}_1},{{\mathbf{w}}_2}, \ldots ,{{\mathbf{w}}_J}]^T}$ is the coefficient matrix and ${\mathbf{s}}(k) = {[{s_1}(k),{s_2}(k), \ldots ,{s_J}(k)]^T}$ is the matrix of slow features.

The coefficients W can be derived by solving the following generalized eigenvalue problem [20]

$$\begin{align}{\mathbf{AW}} = {\boldsymbol{{\textbf{BW}}\Omega }}\end{align} \tag{ 6 }$$

where ${\mathbf{A}} = {\langle {\boldsymbol{\dot{\textbf{x}}}}{{\boldsymbol{\dot{\textbf{x}}}}^T}\rangle _k}$ denotes the covariance matrix of the derivative of the input x, and ${\mathbf{B}} = {\langle {\mathbf{x}}{{\mathbf{x}}^T}\rangle _k}$ denotes the covariance matrix of x, and $\boldsymbol{\Omega} = {\text{diag}}\left\{ {{\omega _1},{\omega _1}, \ldots ,{\omega _J}} \right\}$ is a diagonal matrix with ${\omega _j} = {\langle \dot s_j^2\rangle _k}$, which are also the optimal values of the objective function in formula (2). The eigenvalues ω in Ω are arranged from small to large and represent different slownesses. The smaller the eigenvalue, the slower the corresponding feature. According to different slownesses, s can then be further divided into two groups [24]

$$\begin{align}\left\{ \begin{array}{l} \displaystyle{\mathbf{s}}\ \,=\, [\mathbf{s}_d^T\quad \mathbf{s}_e^T]^T \\ \displaystyle \mathbf{s}_d = \mathbf{W}_{1:R}\mathbf{x} \\ \displaystyle \mathbf{s}_e = \mathbf{W}_{R + 1;J}\mathbf{x} \\ \end{array}\right.\end{align} \tag{ 7 }$$

where R is the number of the dominant slowest features. In this work, it is set to the number of slow features that are slower than the 0.4-upper quantile of the set {Δx(j)}. s _d measures the dominant slow variations of input signals and captures the general trend of process variations. s _e represents residuals with fast variations that can be seen as short-term noises or fluctuations [25].

The idea of the multiple SFA method is that data relating to different flow statuses have different characteristics, so the gas–liquid two-phase flow process can be divided into different modes according to different flow regimes. With prior process knowledge, multiple SFA-based local models can be established for different flow regimes. The steps are as follows:

Step 1: Obtain process data X(N× J) containing N sampling points and J process variables.

The gas–liquid two-phase flow process is divided into several modes according to different flow regimes. The data from different flow regimes are then divided into X₁(N₁ × J), X₂(N₂ × J), ..., and X _M (N_M × J), each of which, i.e. X _m (N_m × J) (m= 1,2, ..., M), contains the modeling data in the mth flow regime.

Step 2: Establish local models to monitor static variations.

Denoting the first R_m rows of W _m and the rest as W _m,d (R_m × J) and W _m,e ((J−R_m ) × J), respectively, multiple static models are constructed for different flow regimes by performing SFA [23]

$$\begin{align}\left\{ \begin{gathered} {{\mathbf{s}}_{m,d}} = {{\mathbf{W}}_{m,d}}{\mathbf{X}}_m^T \hfill \\ {{\mathbf{s}}_{m,e}} = {{\mathbf{W}}_{m,e}}{\mathbf{X}}_m^T \hfill \\ \end{gathered} \right.\end{align} \tag{ 8 }$$

where s _m,d refers to the derived slowest features that tend to catch the general trend of process variations in the mth flow regime; s _m,e represents the fastest features, which can be seen as short-term noises or fluctuations.

Step 3: Establish local models to monitor dynamics.

Because the slow features vary slowly and reflect the overall steady-state distribution, the derivative of the slow features can reflect the changing trend and dynamic characteristics of the flow statuses. The first-order derivatives of the slow features are considered to detect process dynamics. Multiple dynamic models are constructed for different flow regimes based on ${{\boldsymbol{\dot{\textbf{s}}}}_{m,d}}$ and ${{\boldsymbol{\dot{\textbf{s}}}}_{m,e}}$

$$\begin{align}\left\{ \begin{gathered} {{\boldsymbol{\Omega}}_{m,d}} = {\langle {{{\boldsymbol{\dot{\textbf{s}}}}}_{m,d}}{\boldsymbol{\dot{\textbf{s}}}}_{m,d}^T\rangle _k} \hfill \\ {{\boldsymbol{\Omega }}_{m,e}} = {\langle {{{\boldsymbol{\dot{\textbf{s}}}}}_{m,e}}{\boldsymbol{\dot{\textbf{s}}}}_{m,e}^T\rangle _k} \hfill \\ \end{gathered} \right..\end{align} \tag{ 9 }$$

Step 4: Calculate statistics and control limits.

To measure the static variations of flow status, the T² statistics of the mth model are defined based on s _m,d and s _m,e

$$\begin{align}\left\{ {\begin{array}{*{20}{c}} \displaystyle{T_{m,d}^2 = \textbf{s}_{m,d}^T{\textbf{s}_{m,d}}} \\ \displaystyle{T_{m,e}^2 = \textbf{s}_{m,e}^T{\textbf{s}_{m,e}}} \end{array}} \right..\end{align} \tag{ 10 }$$

Because of the non-Gaussian distribution of the slow features in SFA, the control limits of the T² statistics in the mth status, $T_{m,d\lim }^2$ and $T_{m,e\lim }^2$, are calculated using kernel density estimation (KDE) [21]. The confidence level is chosen as α= 0.99 in this work. The T² statistics reflect whether the current status conforms to the model.

To detect process dynamics, two S² statistics are defined based on ${{\boldsymbol{\dot{\textbf{s}}}}_{m,d}}$ and ${{\boldsymbol{\dot{\textbf{s}}}}_{m,e}}$:

$$\begin{align}\left\{ {\begin{array}{*{20}{c}} \displaystyle{S_{m,d}^2 = {\boldsymbol{\dot{\textbf{s}}}}_{m,d}^T{\boldsymbol{\Omega }}_{m,d}^{ - 1}{{{\boldsymbol{\dot{\textbf{s}}}}}_{m,d}}} \\ \displaystyle{S_{m,e}^2 = {\boldsymbol{\dot{\textbf{s}}}}_{m,e}^T{\boldsymbol{\Omega }}_{m,e}^{ - 1}{{{\boldsymbol{\dot{\textbf{s}}}}}_{m,e}}} \end{array}} \right..\end{align} \tag{ 11 }$$

The control limits in the mth flow regime, $S_{m,d\lim }^2$ and $S_{m,e\lim }^2$, are also calculated using KDE. The S² statistics reflect whether the dynamics of the current status and the corresponding model status are consistent.

In status identification, Bayesian inference is developed using the T² static statistics to obtain an intuitive index for identifying the current flow status, while the S² dynamic statistics are used to detect process dynamics. The online status identification and analysis steps are summarized as follows:

Step 1: Collect the real-time data x_on, and then normalize the data to get ${{\boldsymbol{\bar{\textbf{x}}}}_{{\text{on}}}}$

$$\begin{align}\left\{ \begin{gathered} {{\bar x}_{{\text{on}},j}} = \frac{{{x_{{\text{on}},j}} - {{\bar v}_{m,j}}}}{{{s_{m,j}}}} \hfill\\ {{\bar v}_{m,j}} = \frac{1}{{{N_m}}}\sum\limits_{n = 1}^{{N_m}} {{x_{m,nj}}} \hfill\\ {h_{m,j}} = \sqrt {\frac{1}{{{N_m}}}\sum\limits_{n = 1}^{{N_m}} {({x_{m,nj}} - } {{\bar v}_{m,j}}{)^2}} \\ \end{gathered} \right.\end{align} \tag{ 12 }$$

where x_on,j is the online sampling data for the jth process variable. The mean and the standard deviation of the jth process variable in the mth modeling data matrix are represented by ${\bar v_{m,j}}$ and h_m,j , respectively.

Step 2: The slow features are obtained by projecting ${{\boldsymbol{\bar{\textbf{x}}}}_{{\text{on}}}}$ onto the local model W _m:

$$\begin{align}\left\{ \begin{gathered} {{\mathbf{s}}_{{\text{on}},d}} = {{\mathbf{W}}_{m,d}}{\boldsymbol{\bar{\textbf{x}}}}_{{\text{on}}}^T \hfill \\ {{\mathbf{s}}_{{\text{on}},e}} = {{\mathbf{W}}_{m,e}}{\boldsymbol{\bar{\textbf{x}}}}_{{\text{on}}}^T \hfill \\ \end{gathered} \right..\end{align} \tag{ 13 }$$

Step 3: Calculate the T² online monitoring statistics.

$$\begin{align}\left\{ {\begin{array}{*{20}{c}} \displaystyle{T_{{\text{on}},d}^2 = {\mathbf{s}}_{{\text{on}},d}^T{{\mathbf{s}}_{{\text{on}},d}}} \\ \displaystyle{T_{{\text{on}},e}^2 = {\mathbf{s}}_{{\text{on}},e}^T{{\mathbf{s}}_{{\text{on}},e}}} \end{array}} \right..\end{align} \tag{ 14 }$$

The two statistics are used to measure whether the current status corresponds to the status of the mth model. If either of them (or both) exceeds its control limit, this indicates that the status of gas–liquid two-phase flow has changed.

Step 4: Bayesian inference is developed using the T² statistics.

Here, the two T² statistics and the two control limits are combined to provide an intuitive and comprehensive index of the process status [26]. The event of the mth typical flow regime is represented by L_m . In the Bayesian inference of this work, the within-class conditional probability of x_on with respect to $T_m^2 \in \{ T_{m,d}^2,T_{m,e}^2\}$ is computed as:

$$\begin{align}\left\{ {\begin{array}{*{20}{c}} \displaystyle{{\text{Pr}}\left[ {T_{m,d}^2|{L_m}} \right] = \exp \left\{ { - \frac{{T_{m,d}^2}}{{T_{m,d\lim }^2}}} \right\}} \\[12pt] \displaystyle{{\text{Pr}}\left[ {T_{m,e}^2|{L_m}} \right] = \exp \left\{ { - \frac{{T_{m,e}^2}}{{T_{m,e\lim }^2}}} \right\}} \end{array}} \right..\end{align} \tag{ 15 }$$

Correspondingly, the without-class conditional probability of x_on with respect to $T_m^2 \in \{ T_{m,d}^2,T_{m,e}^2\}$ is computed as

$$\begin{align}\left\{ {\begin{array}{*{20}{c}} \displaystyle{{\text{Pr}}\left[ {T_{m,d}^2|{{\bar L}_m}} \right] = \exp \left\{ { - \frac{{T_{m,d\lim }^2}}{{T_{m,d}^2}}} \right\}} \\[12pt] \displaystyle{{\text{Pr}}\left[ {T_{m,e}^2|{{\bar L}_m}} \right] = \exp \left\{ { - \frac{{T_{m,e\lim }^2}}{{T_{m,e}^2}}} \right\}} \end{array}} \right..\end{align} \tag{ 16 }$$

The posterior probability is then defined as the identification index

$$\begin{align}\left\{ \begin{gathered} {\text{Pr}}\left[ {{L_m}|T_{m,d}^2} \right] = \frac{{{\text{Pr}}\left[ {T_{m,d}^2|{L_m}} \right]P({L_m})}}{{{\text{Pr}}\left[ {T_{m,d}^2|{L_m}} \right]P({L_m}) + {\text{Pr}}\left[ {T_{m,d}^2|{{\bar L}_m}} \right]P({{\bar L}_m})}} \hfill \\ {\text{Pr}}\left[ {{L_m}|T_{m,e}^2} \right] = \frac{{{\text{Pr}}\left[ {T_{m,e}^2|{L_m}} \right]P({L_m})}}{{{\text{Pr}}\left[ {T_{m,e}^2|{L_m}} \right]P({L_m}) + {\text{Pr}}\left[ {T_{m,e}^2|{{\bar L}_m}} \right]P({{\bar L}_m})}} \hfill \\ \end{gathered} \right..\end{align} \tag{ 17 }$$

As described above, s _d represents the main information, which catches the general trend of process variations; s _e represents residuals with fast variations and contains very little information. Therefore, the two probabilities in formula (17) are given different weights

$$\begin{align}{\text{Pr}}\left[ {{L_m}|T_m^2} \right] = \lambda {\text{Pr}}\left[ {{L_m}|T_{m,d}^2} \right] + (1 - \lambda ){\text{Pr}}\left[ {{L_m}|T_{m,e}^2} \right]\end{align} \tag{ 18 }$$

where λ (0.5 < λ < 1) and (1 − λ) are the weights assigned to $Pr[{L_m}|T_{m,d}^2]$ and $Pr[{L_m}|T_{m,e}^2]$, respectively (λ is set as 0.8 in this work). The current flow status is then assigned to the status with the largest identification index $Pr[{L_m}|T_m^2]$

$$\begin{align}L_m^* = \mathop {\arg \max }\limits_{{L_m}} \left\{ {{\text{Pr}}\left[ {{L_1}|T_1^2} \right],{\text{Pr}}\left[ {{L_2}|T_2^2} \right], \ldots ,{\text{Pr}}\left[ {{L_M}|T_M^2} \right]} \right\}.\end{align} \tag{ 19 }$$

Step 5: Calculate the S² online monitoring statistics.

$$\begin{align}\left\{ {\begin{array}{*{20}{c}} \displaystyle{S_{{\text{on}},d}^2 = {\boldsymbol{\dot{\textbf{s}}}}_{{\text{on}},d}^T{\boldsymbol{\Omega }}_{m,d}^{ - 1}{{{\boldsymbol{\dot{\textbf{s}}}}}_{{\text{on}},d}}} \\ \displaystyle{S_{{\text{on}},e}^2 = {\boldsymbol{\dot{\textbf{s}}}}_{{\text{on}},e}^T{\boldsymbol{\Omega }}_{m,e}^{ - 1}{{{\boldsymbol{\dot{\textbf{s}}}}}_{{\text{on}},e}}} \end{array}} \right..\end{align} \tag{ 20 }$$

The two statistics are used to measure consistency with the process dynamics. Once either of them (or both) exceeds its threshold, the dynamic characteristics of the current process do not correspond to the model.

The concurrent modeling and identification strategy for the steady states and process dynamics is shown in figure 1. The four possible results and descriptions are listed in table 1.

Zoom In Zoom Out Reset image size

Table 1. Process status as indicated by monitoring statistics.

Case	$T_{m,d}^2$ $T_{m,e}^2$	$S_{m,d}^2$ $S_{m,e}^2$	Process status
1	√	√	Process is in the flow status corresponding to the model
2	√	×	Process is in the flow status corresponding to the model, but some different dynamic characteristics are detected
3	×	√	Process is not in the flow status corresponding to the model but has some similar dynamic characteristics
4	×	×	Process is not in the flow status corresponding to the model

('√' refers to the situation whereby both statistics are within control limits;'×' refers to the situation whereby one or both statistics are outside control limits)

In gas–liquid two-phase flow processes, the flow status is closely related to the pressure at the pipe inlet, the temperature, the velocity of each phase fluid, the medium distribution, and the phase holdup. The sensors used to measure gas-water two-phase flow are shown in table 2.

A conductance sensor is composed of six ring electrodes, and the potential difference is measured between each pair of electrodes. In the experiments, changes to the flow structure lead to changes in the mixture's conductivity and cause changes in electric potential difference. Therefore, the potential difference can be used to calculate the water holdup [7]. The measurement data from the conductance sensor are preprocessed according to formula (21) and selected as modeling data, which reflect the phase holdup information of the gas–liquid two-phase flow when the conductive phase is continuous.

$$\begin{align}{V_{\text{n}}} = \frac{{{V_{\text{w}}}}}{{{V_{{\text{meas}}}}}}\end{align} \tag{ 21 }$$

where V_w denotes the value when the tube is filled with water, V_meas denotes the actual measurement data, and V_n denotes the normalized voltage that reflects phase holdup.

A capacitive sensor consists of two electrode plates. The dielectric constant of the mixed fluid is determined by the dielectric constants of each phase and the phase distribution. There is a corresponding relation between the measured capacitance and phase holdup [8]. The measurement data of the capacitive sensor are selected as modeling data to reflect the phase holdup information when the non-conductive phase is continuous. They are complementary to the data of the conductance sensor. These data are preprocessed according to formula (22)

$$\begin{align}{\text{RCD}} = \frac{{{V_{{\text{meas}}}} - {V_{\text{I}}}}}{{{V_{\text{w}}} - {V_{\text{I}}}}}\end{align} \tag{ 22 }$$

where RCD denotes the relative change of capacitance and V_meas, V_I, and V_w denote the voltage of the fluid, the voltage when the pipe is filled with non-conductive phase, and the voltage when the pipe is filled with water, respectively.

The transmitting piezoelectric chip (TPC) of a continuous-wave ultrasonic Doppler (CWUD) sensor emits a continuous sinusoidal ultrasonic wave, and the receiving piezoelectric chip (RPC) receives the wave reflected by dispersed phase droplets. According to the Doppler effect, the average Doppler shift and the average velocity of the dispersed phase in the sample volume are obtained by calculating the Doppler shift of the ultrasonic signal between the TPC and RPC [27]. The Doppler shift data of the CWUD sensor are selected as modeling data to reflect the flow velocity information.

A pressure gauge is installed at the inlet of the horizontal pipe to measure the pressure at the pipe's inlet. The pressure data are selected as the modeling data that mainly reflect the flow environment information.

A resistance array consists of 16 electrodes and provides two-dimensional medium distribution information for the two-phase flow. The boundary voltage data of the different electrodes are preprocessed [28] according to formula (23) and selected as modeling data.

$$\begin{align}{V_{Ri}} = \frac{1}{n}\sum\limits_{j = 1}^n {({V_{ij}} - {V_{ij0}})} \end{align} \tag{ 23 }$$

where V_ij represents the jth (j= 1, 2, ..., n) boundary voltage measured in the ith (i= 1, 2, ..., 16) excitation direction, V_{ij 0} represents the jth boundary voltage measured in the ith excitation direction on the condition that the pipe is full of water, and V_Ri represents the characteristic value of boundary voltages in the ith excitation direction.

Since the ambient temperature changes very little and has little influence on the gas–liquid two-phase flow process, the temperature is not selected.

The experimental data are collected from a horizontal pipe of the gas–liquid two phase flow loop. The experimental device and test section are shown in figure 2. In the experiments, the different flow ratios of each phase are adjusted to form various flow statuses within the experimental device. The experimental data are then obtained and multiple SFA-based local models can be established.

Zoom In Zoom Out Reset image size

In the experiments, the liquid phase is tap water (with a density of 998 kg m⁻³ and a dynamic viscosity of 1.01 × 10⁻³ Pa•s) and the gas phase is dry air (with a density 1.2 kg m⁻³ and a dynamic viscosity of 1.81 × 10⁻⁵ Pa•s). The horizontal pipe made of stainless steel is 16.6 m in total length with an inner diameter of 50 mm. Standard single-phase flow meters are used to measure each phase in front of the mixer inlet. The test section where the combined sensors are installed is located 12 m downstream from the inlet to allow the flow status to fully develop. A camera is used to directly observe and record the flow status. Besides, the pressure gauge and temperature gauge are installed to monitor the flow environment. By controlling the flow rate of the gas and liquid before they reach the mixer, six typical flow regimes and transition statuses are formed, namely, bubble flow, plug flow, slug flow, wave flow, stratified flow, and annular flow [29].

Data acquisition is carried out for the six typical flow regimes and transition statuses between different flow regimes. The collection time of each data set is 10 s. The acquisition frequencies of the conductance sensor, capacitance sensor and pressure gauge are 2 kHz, and the acquisition frequency of the CWUD sensor is 50 kHz. Two thousand eight hundred frames of data can be measured by the resistance array sensor in 10 s. In this way, a set of two-dimensional data can be obtained every 10 s. Each set includes three voltage data columns from the conductance sensor, two voltage data columns from the capacitance sensor, one Doppler shift data column from the CWUD sensor, one voltage data column from the pressure gauge, and 16 data columns from the resistance array sensor. The data have to be resampled to realize time matching because of inconsistent sampling frequencies. After resampling and preprocessing, each data set is collected as a two-dimensional data array Y(2800 × 23). For each typical flow regime, three sets of data are used for modeling, and two sets of data are reserved as test data.

Based on the modeling data of the gas–liquid two-phase flow, six typical flow regimes are modeled according to the method presented in section 2. The coefficient matrices W _m,d , W _m,e and the diagonal matrices Ω _m,d , Ω _m,e are obtained as model matrices. Besides, the control limits for each model are calculated. The first four slow features extracted for each typical flow regime are presented in figure 3.

Zoom In Zoom Out Reset image size

Figure 3 suggests that the slow features of typical flow regimes are different and represent the inherent essential properties, which provides a theoretical basis for the status identification of the two-phase flow process based on multiple slow feature analysis. Five hundred sampling points of each flow regime are used for verification, and the results are shown in figure 4. Four statistics stay below the control limits for most of the sampling points, verifying the validity and accuracy of the models.

Zoom In Zoom Out Reset image size

During online status identification, the data should be collected continuously in real time. Moreover, in the two-phase flow in an actual industrial process, the most common flow regimes are bubble flow, plug flow and slug flow. In the experiment, the amounts of data for each flow regime are also different. Considering the flow regimes of two-phase flow in actual industrial processes and the amounts of data for each flow regime obtained in the experiment, the test sample is set up as follows: the first ∼400 sampling points are in bubble flow, sampling points 401 ∼ 900 are in plug flow, sampling points 901 ∼ 1300 are in slug flow, sampling points 1301 ∼ 1600 are in wave flow, sampling points 1601 ∼ 1900 are in stratified flow, and sampling points 1901 ∼ 2200 are in annular flow. Bayesian inference is developed to obtain an intuitive identification index. For the PCA and ICA methods, according to minimum squared prediction error (SPE) [30], the SPE statistics are used to calculate Bayesian inference-based indexes. The identification results from the multiple PCA method, the multiple ICA method and the proposed method are presented in figures 5–7, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In figures 5–7, flow statuses one, two, three, four, five and six on the vertical axis of the figures respectively denote bubble flow, plug flow, slug flow, wave flow, stratified flow and annular flow. The multiple SFA method has the highest accuracy among the three methods. In general, the flow status identification results using the proposed method are satisfied, with an overall false identification rate of 1.64%. The false identification rates for different flow regimes of the three methods are presented in table 3.

For the multiple PCA method, the false identification rate is high for all the flow statuses except for slug flow and annular flow. Taking the first 400 sampling points, for example, numerous sampling points are misidentified as plug flow and slug flow (as in figure 5). The detailed monitoring statistics are shown in figure 8, which represent the monitoring results using the slug flow model.

Zoom In Zoom Out Reset image size

In figure 8, the T² statistics are below the control limit for almost all of the first 1300 sampling points, indicating that there is no change of flow status. Therefore, the T² statistics barely detect the change of status before the one thousand three hundredth sampling point. Meanwhile, SPE statistics cannot distinguish bubble flow from slug flow. After Bayesian inference, the identification result from the multiple PCA method is disappointing.

For the multiple ICA method, the false identification rate is high for plug flow (as shown in figure 6). To analyze the identification results more meticulously, monitoring statistics using the models of plug flow and slug flow are shown in figures 9 and 10, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In figure 9, none of the two statistics in the first 400 sampling points exceed the control limit, and some of them are even smaller than those of the four hundred and first to nine hundredth sampling points. This implies that the change of status from bubble flow to plug flow cannot be detected. In figure 10, it can be seen that the monitoring results from ICA are better than those of PCA (as shown in figure 8). However, two statistics are below the control limits for some sampling points in the four hundred and first to nine hundredth sampling points, which means that some sampling points in the plug flow are misidentified as being in slug flow.

The process of multiphase flow tends to be dynamic, that is, process data are associated with time. However, PCA and ICA assume that process data are statically independent, causing the dynamics of the process to be ignored. This means that it is difficult to accurately describe the actual variations of the data, and therefore, the identification effect is not ideal.

For the multiple SFA method, some sampling points in the bubble flow are misidentified as plug flow (as in figure 7). Monitoring statistics are used to analyze the flow process in the first 900 sampling points more meticulously and monitor the process dynamics. The monitoring results using the multiple SFA-based models for bubble flow and plug flow are shown in figures 11 and 12, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the plot of figure 11(a), the T² statistics detect the change of flow status after the four hundred and first sampling point. In the plot of figure 11(b), from the four hundred and first to the one thousand three hundredth sampling points, a few of the $S_d^2$ and $S_e^2$ statistics are below the control limits. This indicates that although the status has changed, there are still some similar dynamic characteristics. In figures 12(a) and (b), the monitoring results of the SFA are better than those of ICA (as in figure 9). The $T_d^2$ statistics that capture the main variations are around the control limit before the four hundredth sampling point, and they are higher than those of the four hundred and first to eight hundredth sampling points; some of them exceed the control limit. The S² statistics reveal some similarity of dynamic characteristics between bubble flow, plug flow, and slug flow. Because there are small bubbles in plug flow and slug flow, the flow process exhibits similar properties at certain periods. The monitoring results show that the proposed method can not only identify typical flow regimes but also accurately detects similar dynamic characteristics in different flow statuses.

In gas–liquid two-phase flow processes, the transition statuses between different flow regimes are common and important, have an important influence on the flow process, and should be paid attention to. The experimental data of the transition status from bubble flow to slug flow (as shown in figure 13) are used for testing. The first stage contains 500 sampling points in bubble flow. The second stage contains three data sets representing the transition of the flow status from bubble flow to slug flow, each with 200 sampling points. The third stage contains 500 sampling points in slug flow. In the experiment, the change from bubble flow to slug flow is the result of an increase in the gas flow rate.

Zoom In Zoom Out Reset image size

Bayesian inference is developed to identify the flow status. Because the flow process fluctuates complexly and violently in the transition status, a time sliding window L is used to smooth the Bayesian inference-based index. The transition status and typical flow regimes can then be identified by the following criteria:

$$\begin{align}\left\{ \begin{gathered} \Pr (i) < \beta \;\;\;\;\;({\text{transition}}\;{\text{status}}) \hfill \\ \Pr (i) \geqslant \beta \;\;\;\;\;({\text{typical}}\;{\text{flow}}\;{\text{regime}}) \hfill \\ \end{gathered} \right.\end{align} \tag{ 24 }$$

where Pr(i) denotes the Bayesian inference-based index of bubble flow or slug flow, and β (0 < β< 1) is the threshold (β is set to 0.5 in this work).

The identification results using multiple SFA combined with Bayesian inference are shown in figure 14.

Zoom In Zoom Out Reset image size

In the plot of figure 14, flow statuses one, two and three respectively denote bubble flow, transition, and slug flow. It can be seen that the Bayesian inference-based index of bubble flow has a decreasing tendency, while that of slug flow has an increasing tendency. The transition status is from the four hundred and eighty-fifth to the one thousand and fourteenth sampling points.

To monitor the steady states and process dynamics of the changing process accurately and meticulously, statistics are used for analysis. Two other methods are compared with the proposed method. The distributions of statistics using the three methods are shown in figures 15–17, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the plot of figure 15, the T² statistics cause false alarms at the thirtieth and the two hundred and eighty-sixth sampling points. From the five hundred and first to the seven hundredth sampling points, a few of the T² and SPE statistics exceed the control limits, while there are still many statistics below the control limits. This indicates that the flow status of the process is changing. From the seven hundred and first to the thousandth sampling points, numerous T² and SPE statistics are outside the control limits, while there are a few statistics below the control limits, indicating that the current process is in transition status and has both bubble flow characteristics and slug flow characteristics. After the one thousand and first sampling point, both indexes exceed the control limits for most of the sampling points, indicating a complete change in flow status. In practice, the fluid is mixed with numerous small bubbles in slug status. However, the principal component information extracted by PCA only shows that there is a change in status, with no more information provided.

In the plot of figure 16, the S² and SPE statistics indicate the process is in bubble flow for the first 500 sampling points. As opposed to the results from multiple PCA, both statistics seldom exceed the control limits during the five hundred and first to seven hundredth sampling points, in which the initial change has been ignored. After the seven hundred and first sampling point, the distribution of statistics is similar to that using multiple PCA method. The monitoring results of the above two methods lead to the same conclusion, which is that they can identify the process status, but they cannot distinguish the dynamic potential changing trend or the similar dynamic characteristics among different flow statuses.

In the plot of figure 17, the steady states monitored by $T_d^2$ and $T_e^2$ are illustrated in figure 17(a), and the process dynamics described by $S_d^2$ and $S_e^2$ are plotted in figure 17(b). There is an increasing tendency in the statistics. Two of the T² statistics stay below the control limits for the first 500 sampling points, illustrating that the flow status is bubble flow, while the process appears to operate around the control limit of $S_d^2$. The $S_d^2$ statistics exceed the control limit at around the fifty-first, the two hundred and thirty-ninth and the three hundred and ninety-ninth sampling points, indicating that the dynamic variation has been noticed before the transition status. From the five hundred and first to the seven hundredth sampling points, a few $T_d^2$ statistics exceeding the control limits detect the initial change, which is earlier than PCA and ICA. However, the S² statistics remain similar to the distribution of the first 500 sampling points, which indicates that although the status is changing, the dynamic characteristics of the current process have not changed much and the process is in the initial stage of transition status. This similarity can be seen in the first two photos in figure 13; small bubbles gradually coalesce into larger ones. From the seven hundred and first to the thousandth sampling points, the T² and S² statistics exceed the control limits for some sampling points, while there are still some statistics under the control limits, illustrating that the process is in transition status. After the one thousand and first sampling point, the two T² statistics exceed the control limits for most of the sampling points, indicating that the flow status has changed completely to slug flow. Some of the S² statistics are still below the control limits because of the similar dynamic characteristics of bubble flow and slug flow.

From the above analysis, PCA and ICA can only tell whether the change of flow status has been detected but cannot provide dynamic information or detect similar characteristics in different flow statuses. The flow status changes over time, so the dynamic characteristics of the time-varying flow process should be fully considered. Moreover, although there are great differences among the flow statuses of a gas–liquid two-phase flow, there are still some similarities in finite time that should not be ignored. Unlike PCA or ICA, the proposed method provides more meaningful information by concurrently monitoring both the steady states and the process dynamics. SFA extracts information of two characteristics, which are the slow features and the first-order derivative of the slow features. The information extracted by SFA is more comprehensive and the physical significance is clearer. Therefore, the identification and analysis made by SFA are more in line with the actual situation and have better interpretability.