1 Introduction

Gas–liquid two-phase flow refers to a mixed flow consisting of two immiscible media: gas phase and liquid phase. It widely exists in production and transportation of petroleum industry. Flow behavior and flow parameters are important research issues in gas–liquid flow (Chen 2011). Many studies have been made to explore gas–liquid flows, which plays an important role in the production assessment of oil and gas wells. Zhu et al. (2012) focused on the mechanism of sustained casing pressure (SCP) in gas wells. Luo et al. (2014) applied separator control to inhibit severe slugging. Yin et al. (2018) analyzed heat transfer for gas–liquid flow in vertical wellbore annuli. However, due to the changeable flow process and diverse flow structures, gas–liquid flow can be considered as a typical complex system. The flow behavior underlying such flow situation is still blurry, and flow parameters are time varying and difficult to measure. The whole flow process presents obvious instability, instantaneity and randomness, which also brings series of difficulties to the parameter measurement. Characterizing gas–liquid flow and measuring flow parameters represent challenges of great importance, which contribute to the recognition of flow regime and the optimal design of industrial equipment (Qi et al. 2012).

Over the past few decades, a mass of work has been done for exploring gas–liquid flow behavior and brought us many achievements. Qi et al. (2013) employed CFD-PBE model to investigate the hydrodynamics of multiphase flow and accurately predict radial profiles of local gas holdup and bubble diameter. Gao et al. (2016a, b, c, d) characterized the slug to churn flow transition in a 20-mm-diameter vertical pipe by using multivariate pseudo Wigner distribution and multivariate multiscale entropy. Li et al. (2018) revealed the nonlinear dynamic characteristics of gas–liquid flow in small channels. Their analysis found obvious differences in the dynamic characteristics between different flow patterns in small channels, and the sensitivity varies with the change of scale. Furthermore, many researchers were dedicated to the measurement of void fraction, which can be used to determine mixture density and estimate gas or petroleum reserves. Some studies aimed to establish some formulas to describe the correlation of flow structures and response characteristics of sensors. Sardeshpande et al. (2015) utilized electrical capacitance tomography (ECT) facility to predict gas void fraction. He et al. (2019) developed a multi-wire capacitance probe to obtain the average void fraction and the local void fraction. Pan et al. (2018) proposed a slip-ratio based equal-diameter double-circle model for void fraction measurement. In addition, some studies applied machine learning to fit the mapping relationship between sensor response and flow parameters. Li et al. (2016) identified flow patterns via Fisher discriminant analysis and then measured the void fraction via support vector machine in small channels. Wang et al. (2019) combined empirical mode decomposition and artificial neural networks for void fraction estimation. Dang et al. (2019) developed a deep learning method including convolutional neural networks and long short-term memory networks for flow parameters measurement. Despite these existing works, it is still an important challenge to accurately measure flow parameters, due to complex and variable flow conditions, such as media, pipes and temperature. Moreover, complicated flow behavior also affects the measurement of flow parameters. It is in urgent need to develop a method which can effectively characterize gas–liquid flow behavior and achieve accurate measurement of flow parameters.

Nowadays, complex network theory plays an important role in characterizing complex systems. Its main idea is to map a real-world system to a network, where the nodes denote the components of the system and the edges allow describing the relationship among these components. It provides an effective framework for better understanding of complex systems and has made great progress in various fields (Newman 2003; Zhang et al. 2008; Xu et al. 2008; Donges et al. 2011; Hong et al. 2012; Wang et al. 2016; Gao et al. 2018). In particular, complex network theory can effectively extract topological characteristics from the time series. A review can be seen in the literature (Gao et al. 2016a, b, c, d). Zhang and Small (2006) built complex networks from pseudo-periodic time series by considering each cycle as a node, which demonstrates that time series with different dynamical properties have different topological structures. Lacasa et al. (2008) and Luque et al. (2009) proposed the visibility graph and horizontal visibility graph, which allow mapping a time series to a complex network fast and simply. Gao et al. (2016a, b, c, d) proposed a limited penetrable visibility graph (LPVG) to analyze nonlinear time series, containing EEG signals and two-phase flow signals. The LPVG has great noise resistance performance, which enables to extract features underlying experimental measurements.

More recently, deep learning methods (LeCun et al. 2015) have been widely used in computer vision, speech recognition and natural language processing and have achieved great success. Since AlexNet (Krizhevsky et al. 2017) showed an excellent performance on ImageNet, convolutional neural network received a lot of attention and achieved considerable development in recent layers. He et al. (2016) employed shortcut connections to add previous feature maps to current outputs, which can solve the issue of the disappearance of the gradients and easily train a deeper neural network. Further, Huang et al. (2017) designed a dense convolutional network which introduces direct connections from any layer to all subsequent layers instead of shortcut connections, achieving a state-of-the-art performance with fewer parameters and less computation. More and more researchers applied deep learning technology to decode time series due to these advantages. Gao et al. (2019) proposed a spatial–temporal convolutional neural network to evaluate the fatigue conditions of drivers. Yao et al. (2020) combined convolutional neural network, recurrent cells and attention module for information fusion on ECG signals. Sors et al. (2018) employed a 14-layer convolutional neural network for sleep stage detection via single-channel raw EEG.

In this paper, we propose a novel complex network-based deep learning method to characterize gas–liquid flow. Based on multichannel measurements acquired by four-sector distributed conductance sensor, we first construct multiple limited penetrable visibility graphs (LPVGs) and obtain degree sequences as the graph representation for flow behavior characterization. Then, we design a dual-input convolutional neural network to fuse the raw signals and the graph representation for the classification of flow structures and measurement of void fraction. Particularly, we implement the model with two parallel branches with the same structure, each corresponding to one input. Each branch consists of a channel-projection convolutional part, a spatial–temporal convolutional part, a dense block and an attention module. The outputs of the two branches are concatenated and fed into several full connected layers for flow structure classification and void fraction measurement.

2 Materials and methods

2.1 Experiment and data acquisition

The vertical upward gas–liquid two-phase flow experiment was carried out based on the multiphase flow loop facility, which comprises the following parts: a vertical testing pipe, our designed four-sector distributed conductance sensor (Gao et al. 2016a, b, c, d), a water tank, a water pump, an air pump and a gas meter as shown in Fig. 1. A vertical testing pipe with 50 mm inner diameter was used, and the tap water and the air were taken as the liquid phase and the gas phase in our experiment. The four-sector distributed conductance sensor consists of eight titanium alloy concave electrodes, which flush mounted on the inside pipe wall. There are two electrode layers with a pitch of 10 mm. The electrode height is 10 mm, and the electrode angle is 45°. The sensor was employed to obtain gas–liquid two-phase flow information, and its installation location was chosen elaborately so that the gas–liquid flow at the sensor can be developed completely. The water pump and the gas meter were used to precisely adjust the inlet flow parameters, i.e., the water velocity and the gas velocity, in our experiment. The water pump and the gas meter were calibrated before experiment.

Fig. 1
figure 1

Schematic of vertical upward gas–liquid two-phase flow loop facility

Full size image

In order to generate different gas–liquid flow conditions, we pumped the two phases (gas and liquid) into the pipe at many different velocities. In the experiment, we studied the gas velocities in the range from 0.014 to 0.212 m/s and the water velocities in the range from 0.002 to 0.153 m/s. Then, we got 18 flow conditions in total and the gas void fraction is in the range from 12.2% to 95.4%. For each flow condition, four-channel signals were recorded by the four-sector distributed conductance sensor and stored by PXI-4472 synchronization data acquisition board card of NI Company. The sample rate was 2 kHz, and the sampling duration was 20 min. Meanwhile, we equipped a high-speed digital camera for flow pattern definition and recording. During the experiment, we can observe six representative gas–liquid flow structures, i.e., uniform bubble flow (UB), bubble flow with small bubbles (B-SB), bubble flow with large bubbles (B-LB), bubble flow with high velocity (B-HV), slug flow wrapped in bubbles (S-WB) and slug flow with large slugs (S-LS). Figure 2 shows the pictures of these flow structures, which are obtained by the high-speed digital camera.

Fig. 2
figure 2

Experimental snapshots of representative flow structures, i.e., uniform bubble flow (UB), bubble flow with small bubbles (B-SB), bubble flow with large bubbles (B-LB), bubble flow with high velocity (B-HV), slug flow wrapped in bubbles (S-WB) and slug flow with large slugs (S-LS)

Full size image

2.2 Methodology

In this section, we first introduce the construction of multiple LPVGs and consider the degree sequences as the graph representation. Then, we present the dual-input convolutional neural network with the raw signals and degree sequences as inputs. Finally, we describe tasks for the classification of flow structures and measurement of void fraction in gas–liquid two-phase flow.

2.2.1 Limited penetrable visibility graph and degree sequence

Based on the above multichannel measurements collected from the four-sector distributed conductance sensor, we map the measurements to complex networks via limited penetrable visibility graph. For the time series per channel, a LPVG whose limited penetrable distance is 1 is constructed. The LPVG is developed from visibility graph, which is an effective method for mapping measurements with noise to complex networks. Figure 3 shows a schematic diagram from the time series to the LPVG.

Fig. 3
figure 3

Example of a the time series and b its corresponding limited penetrable visibility graph (LPVG), where every node corresponds to series data point in the same order. The blue lines represent the edges of visibility graph, while the red lines represent the added edges of LPVG with the limited penetrable distance being 1

Full size image

For a time series \(\left\{ {x_{t} } \right\}_{t = 1}^{L}\), where L is the length of the time series, we treat each time point as a node and the edges among these nodes are assigned based on whether two nodes can see each other. For visibility graph, we first consider the points of the time series as some vertical bars whose heights are the values of each time point and whose positions are the same as the time series as shown in Fig. 3a. Then, we link every two bars and consider they are visible if the visibility line (i.e., the connection line between two vertices) between them is not blocked as shown by the blue lines. For example, the point o and the point p can see each other, so the edge between them can be inferred. While the visibility line between the point o and the point q is blocked by the point p, so there is no connection between o and q. More formally, for two arbitrary data points (\(t_{a}\), \(v_{a}\)) and (\(t_{b}\), \(v_{b}\)), there is a connection between them only if any other data (\(t_{c}\), \(v_{c}\)) located between them satisfy:

$$v_{c} < v_{b} + \left( {v_{a} - v_{b} } \right)\frac{{t_{b} - t_{c} }}{{t_{b} - t_{a} }}$$
(1)

where \(t_{i}\) is the time of data points and \(v_{i}\) is the observation values of data points, \(i = a, b, c\).

Compared with the visibility graph, the limited penetrable visibility graph adds a limited penetrable distance. Particularly, if we set the limited penetrable distance to be M, there is a connection between two nodes only if the number of intermediate nodes that cut off the visibility line is no more than M. As shown in Fig. 3, the red lines are the added edges when we set the limited penetrable distance to be 1. At this time, the point p is the only node that cut off the visibility line between the point o and the point q, so an edge can be inferred between o and q.

After that, we can obtain the multiple LPVGs \(A_{i,j}^{k} \in \left\{ {0,1} \right\}\), where \(k = 1, 2, 3, 4\) represents the channel of the sensor. We then extract the degree sequence (DS) \(\left\{ {D_{i}^{k} = \mathop \sum \nolimits_{j = 1,j \ne i}^{L} A_{i,j}^{k} } \right\}_{i = 1}^{L}\) as the graph representation so as to characterize the dynamics of the gas–liquid flow.

2.2.2 Dual-input convolutional neural network

As shown in Fig. 4, we design a dual-input convolutional neural network to fuse the raw signals and the graph representation from multiple LPVGs, which contains feature extractor and classifier. We implement the feature extractor via two parallel branches with the same structure, each corresponding to one input. Firstly, each branch extracts spatial–temporal representation from input signals with four channels and 400 sample index. We perform a channel-projection convolutional layer with the kernel size of (1, 1) to fuse the information between channels. Then, a temporal convolutional layer with kernel size of (1, 11), a spatial convolutional layer with kernel size of (4, 1) and an average pooling layer are applied in turn. Note that, we insert the exponential linear unit (ELU) after each convolutional layer as the activation function and apply a dropout operation after the average pooling layer. In spatial convolutional layer, we use depthwise convolution, which performs regular convolution within each depth slice. The depthwise convolution has a lower parameter amount and operation cost. In our model, we set the spatial filter of depthwise convolutions as one, so that the output of depthwise convolutions has the same size as input.

Fig. 4
figure 4

The architecture of dual-input convolutional neural network, where indigo rectangle represents a single layer and orange rectangle represents a network module

Full size image

After extracting spatial–temporal representation, we employ a dense block to further abstract advanced features. Compared to traditional methods that cascade layers together, the dense connectivity increases the diversity of input features and facilitates the reuse of redundant features. Figure 5 shows the details of dense block, which adopts the same structure as DenseNet. The dense block includes three same convolutional blocks, which perform six operations sequentially: batch normalization, applying ELU activation, convolution with the kernel size of (1, 1), batch normalization, applying ELU activation and convolution with the kernel size of (1, 5). For a convolutional block i, its input \(X_{i}\) can be described as follows:

$$X_{i} = Cat\left( {\left[ {X_{0} ,Y_{1} , \ldots ,Y_{i - 1} } \right]} \right)$$
(2)

where \(X_{0}\) is the input of the dense block, \(Y_{i - 1}\) is the output feature map of convolutional block i − 1, and \(Cat\left( \cdot \right)\) is the concatenating function. To ensure that these feature maps can be concatenated, we set the stride of convolution as (1, 1) and pad 0 around the margin of the convolution input to keep the size of feature maps consistent. We then apply an average pooling layer to further reduce the dimension of feature maps.

Fig. 5
figure 5

The architecture of dense block. Conv block means convolutional block, which performs six operations sequentially: batch normalization, applying ELU activation, convolution with the kernel size of (1, 1), batch normalization, applying ELU activation and convolution with the kernel size of (1, 5)

Full size image

Based on the output of dense block, we then use attention mechanism to weight the features at different time steps. Figure 6 shows the attention module used in our method, which is based on the similar principle as the one used in Raffe et al. (2016). In the attention module, one fully connected layer along with a hyperbolic tangent function is applied to the feature vector at each time step. The output size of fully connected layer is the same as the input feature vector. Softmax is then employed for each series, so that the values in each series sum to 1. For each feature, the attention module generates an independent group of values which would be used as weights. At last, the weighted average of input is computed as the output of the attention module. Given a module which produces a feature vector \(F_{t} \in {\mathbb{R}}^{d \times 1}\) at each time step t, attention module computes an output vector \(O \in {\mathbb{R}}^{d \times 1}\) as the weighted mean of the feature vector F by

$$O = \mathop \sum \limits_{t = 1}^{T} \beta_{t} *F_{t}$$
(3)

where T is the total number of time steps, \(\beta_{t} \in {\mathbb{R}}^{d \times 1}\) is the weight computed at each time step t based on \(F_{t}\) and * represents element-wise multiply. The weight \(\beta_{t}\) can be formulated as follows:

$$h_{t} = { \tanh }\left( {WF_{t} + b} \right)$$
(4)
$$\beta_{t} = \frac{{{ \exp }\left( {h_{t} } \right)}}{{\mathop \sum \nolimits_{k = 1}^{T} { \exp }\left( {h_{k} } \right)}}$$
(5)

where \(w \in {\mathbb{R}}^{d \times d}\) and \(b \in {\mathbb{R}}^{d \times 1}\) are the trainable neural network parameters. By introducing the attention mechanism as the above, we weight the time step separately for each feature. It can help the model focus on important segments of the signal, thereby improving the network performance.

Fig. 6
figure 6

The structure of attention module. The feature vector \(F_{t}\) at each time step t is fed into a fully connected layer to produce a probability vector \(\beta_{t}\), where \(t = 1,2, \ldots ,T\) and T is the total number of time steps. The output vector O is the weighted average of \(F_{t}\) with weighting given by \(\beta_{t}\)

Full size image

For the classifier, we first concatenate the outputs of the two branches and then apply several full connected layers for flow structure classification and void fraction measurement. Particularly, for the task of flow structure classification, we use a full connected layer with nodes of 1024 and a full connected layer with nodes of 6 (the number of flow structure categories), and the activation functions are set as tanh and softmax, respectively. We use categorical cross-entropy loss function as the objective function for model optimization. For the task of void fraction measurement, we use a full connected layer with nodes of 1024, a full connected layer with nodes of 18 (the number of flow conditions) and a full connected layer with nodes of 1. The activation functions are set as tanh, softmax and sigmoid, respectively. Note that, in the task of void fraction measurement, we establish the connection between gas void fraction classification and gas void fraction measurement. Here, we utilize two output layers for classification and prediction, respectively. The total loss is

$$L_{\text{total}} = \alpha L_{\text{cla}} + \beta L_{\text{pred}}$$
(6)

where \(L_{\text{cla}}\) is categorical cross-entropy loss function corresponding to classification of gas void fraction, \(L_{\text{pred}}\) is mean square error loss function corresponding to gas void fraction prediction, and \(\alpha\) and \(\beta\) are the hyper-parameters that control the trade-off between classification objective and prediction objective.

3 Results and discussion

3.1 Characterization of gas–liquid flow via LPVG

For each flow condition, the multichannel measurements are split into segments with length of 3 s via a sliding window, and the sliding step is 1 s. To simplify the calculation, we downsample the segments to 200 Hz. Then, we infer multiple LPVGs from each segment and obtain their degree sequences. In order to characterize various flow structures, we count the distributions of the degree values of the LPVGs from all the segments under the six representative structures. We normalize the degree distributions and calculate the frequency of occurrence of each degree value. Figure 7 shows the degree distributions of six flow structures in one channel of the sensor and that of other channels have the similar distribution. From Fig. 7, it is obvious that the degree distribution under bubble flow is narrower and the maximum frequency of occurrence is higher, while that under slug flow is wider and the maximum frequency of occurrence is lower. What’s more, the degree value corresponding to the peak under slug flow is larger than that under bubble flow. Particularly, when the gas velocity is low while water velocity is high, the mixture flow appears as uniform bubble flow. In this flow condition, gas phase is not able to gather due to the impact of water phase, resulting that a large number of tiny bubbles flow randomly in the pipe. The gas phase distributes uniformly in the pipe, and the behavior of mixture flow represents strong randomness. Thus, the degree values are concentrated on few values and there is a sharp peak in the degree distribution. In addition, the degree value corresponding to the peak is small under uniform bubble flow due to the random flow of a large number of bubbles in small sizes. With the decrease in water velocity, the effect of water phase on the gas phase reduces and some small bubbles with various sizes appear in the pipe. The mixture flow turns into non-uniform bubble flow where bubble size becomes larger. The degree value corresponding to the peak slightly increases, and the degree distribution is similar to uniform bubble flow. When water velocity reaches a smaller value, gas phase can aggregate into big bubbles and there are also some small bubbles scattered in the pipe as shown by bubble flow with large bubbles (B-LB) in Fig. 2. The heterogeneity of mixture flow enhances, and the degree distribution becomes smooth. The degree value corresponding to the peak becomes larger due to the larger bubble size.

Fig. 7
figure 7

The degree distributions of six flow structures (UB, B-SB, B-LB, B-HV, S-WB and S-LS) in one channel of the sensor, where the abscissa represents the degree value, and the ordinate represents the frequency of occurrence of each degree value. D means the degree value corresponding to the peak as shown by the red line

Full size image

When the gas velocity and water velocity are both high, the mixture flow appears as bubble flow with high velocity (B-HV) as shown in Fig. 2. We can see a large number of bubbles with various sizes in the pipe. The degree distribution becomes wider, corresponding to strong heterogeneity of mixture flow. Additionally, the flow of the larger bubbles results in a larger degree value corresponding to the peak. With the increase in the gas velocity, gas phase gradually gathers into slugs and the mixture flow turns into slug flow same as slug flow wrapped in bubbles (S-WB) in Fig. 2. There are a large number of bubbles around slugs. The spatial flow structure in this flow condition is desultory; namely, the mixture flow presents strong transient. The flow behavior represents the intermittent quasi-periodic flow characteristic and heterogeneity. The degree values evenly distribute across multiple values instead of concentrating on a value, and the maximum frequency of occurrence is small. What’s more, the degree value corresponding to the peak is larger than that of bubble flow due to the flow of large slugs. When gas velocity keeps the large value and water velocity decreases, we can see legible large slugs same as slug flow with large slugs (S-LS) in Fig. 2. In this flow condition, there are large slugs and some small bubbles in the pipe. Degree distribution further widens and the degree value corresponding to the peak further increases, which indicate that the quasi-periodicity and heterogeneity characteristics of mixture flow intensify. Therefore, our method can effectively characterize the complicated flow behavior and reveal the evolution from bubble flow to slug flow.

3.2 Classification of flow structures

In this section, we aim to accurately recognize six typical flow structures: UB, B-SB, B-LB, B-HV, S-WB and S-LS; each of them contains three flow conditions. For each 1200-second flow condition, we choose the top 80% measurements as training set, the next 10% measurements as validation set and the last 10% measurements as testing set. Then, the multichannel measurements are split into segments of length 3 s via a sliding window, and the sliding step was 1 s. We infer multiple LPVGs from each segment and obtain their degree sequences. The proposed dual-input CNN is implemented with the inputs of raw signals and degree sequences. We save the optimal model that has the lowest loss in the validate set and assess the classification performance through the test set. The method can achieve an accuracy of 95.3%. Then, we exploit the detection performance for each special flow structure. Figure 8 shows the confusion matrixes where the value on the diagonal represents the accuracy of each flow structure. We can see that all flow structures can achieve a high accuracy over 90%, of which B-LB achieves the highest accuracy of 98.3%, while B-SB obtains a lowest accuracy of 92.2%. Since it is similar between B-SB and UB as well as B-HV, there are some misclassified samples between them. However, there are no misclassified samples between bubble flow and slug flow, indicating a large difference between bubble flow and slug flow. These results suggest that our proposed method can effectively recognize gas–liquid flow structures.

Fig. 8
figure 8

Confusion matrixes of flow structures where the value on the diagonal represents the accuracy of each flow structure

Full size image

3.3 Measurement of void fraction

Besides the characterization of flow behavior and the classification of flow structure, the measurement of gas void fraction is also an important challenge in gas–liquid flow due to complicated flow conditions. The proposed method can not only effectively recognize gas–liquid flow structures, but also accurately measure gas void fraction. In this section, we evaluate the performance of our method in the task of gas void fraction measurement. Note that, we add the classification task of gas volume fraction into the model to assist the regression task. As described in Sect. 2.2, we set the number of neuron nodes as 18 in the penultimate layer, corresponding to the number of flow conditions, and then apply a softmax layer for classification. Furthermore, we cascade a fully connected layer with one neuron node and a sigmoid activation function for gas void fraction measurement. We choose mean squared error (MSE) and mean absolute percent error (MAPE) as indicators describing the measurement performance:

$${\text{MSE}} = \frac{1}{N}\mathop \sum \limits_{i = 1}^{N} \left( {\hat{y}_{i} - y_{i} } \right)^{2}$$
(7)
$${\text{MAPE}} = \frac{1}{N}\mathop \sum \limits_{i = 1}^{N} \left| {\frac{{\hat{y}_{i} - y_{i} }}{{y_{i} }}} \right| \times 100\%$$
(8)

where N is the number of samples, \(\hat{y}_{i}\) represents the predicted value and \(y_{i}\) represents the true value. MSE here describes the absolute error of measurements, and MAPE represents the relative deviations. Our method can reach a MSE of 0.0038 and a MAPE of 6.3%.

4 Conclusions

Characterizing gas–liquid flow and measuring flow parameters represent challenges of great importance. We first capture multichannel measurements of gas–liquid flow via four-sector distributed conductance sensor. We obtain 18 flow conditions in total and six representative gas–liquid flow structures. Then, we map the multichannel measurements to multiple LPVGs under each flow condition and obtain their degree sequences. Based on the analysis of degree distribution, we characterize various flow structures and find that the degree distribution under bubble flow is narrower, while that under slug flow is wider. What’s more, the degree value corresponding to the peak under slug flow is larger than that under bubble flow. After that, a dual-input convolutional neural network with the inputs of raw signals and degree sequences is designed for the classification of flow structures and measurement of void fraction. We implement the model with two parallel branches with the same structure, each corresponding to one input. Each branch consists of a channel-projection convolutional part, a spatial–temporal convolutional part, a dense block and an attention module. The outputs of the two branches are concatenated and fed into several full connected layers for the classification and measurement. Our method achieves an accuracy of 95.3% for the classification of flow structures and a MSE of 0.0038 and a MAPE of 6.3% for the measurement of gas void fraction. Our method is a promising solution for characterizing gas–liquid flow and measuring flow parameters.