Three-dimensional visualisation of gas-water two-phase flow based on bubble mapping method and size projection algorithm

Electrical impedance tomography (EIT) has been successfully applied on gas-water flow applications, but it is incapable to identify small bubbles or the sharp gas-water interface of a large bubble due to its relatively low spatial resolution. A new visualisation approach, bubble mapping method (BM3D), offers a good 3D visualisation of bubble size and distribution. However, the empirical thresholding value method used in BM3D might meet a challenging from various flow setups and conditions in practice. Recently, the size projection algorithm (SPA) was proposed to determine the closest thresholding value for each frame of tomogram by minimising projection error. In this paper, the performances of BM3D and SPA methods are individually analysed and evaluated. Then a new method based on the combination of BM3D and SPA methods is reported to achieve better visualisation of gas-water flow, where the SPA is employed to determine the optimised thresholding values for BM3D method. Experiments are conducted to evaluate the proposed combination method for typical gas-water pipeline flow regimes, including horizontal stratified, bubble, plug, slug, annular flow regimes and vertical bubble, slug, annular flow regimes. The results are compared with the BM3D method, colour mapping method, and highspeed camera video recorded from a transparent chamber. A brief discussion on the effects of reconstruction algorithms and thresholding value for horizontal and vertical flows visualisation is also given.


Introduction
Gas-water two-phase flow is a common and important flow in many industries, where the flow visualisation is of significance for understanding and predicting flow dynamics, process operation, analysis and design of flow control equipment [1,2]. The on-site inspection with a high-speed camera through a transparent chamber might be the most common and direct visualisation method to reveal the gas distribution in water [3][4][5]. However, this method is subject to the availability of the transparent chamber, transparency of continuous fluid. Moreover, high gas void fraction (over 10% [6]) affects the reliability of observation. As an alternative technique, tomographic methods (including optical, ultrasonic or acoustic, γ/X-ray, microwave, resistive or capacitive tomography), are highly considered as a promising technology because of providing 2D/3D images for various multiphase fluids with a feature of "seeing through". Particularly, electrical impedance tomography (EIT) is a non-intrusive and cost-effective visualisation technique with a high temporal resolution (sub-millisecond [7]), which can produce a stack of cross-section images for revealing the disperse phase distribution in water continuous two-phase flows, where the phase difference in conductivity exists. However, tomograms generated by EIT are normally ambiguous with relatively low spatial resolution (up to 5% [8]), which is incapable of identifying very small bubbles and determining sharp gas-water interfaces. However, according to simplified Maxwell relationship as expressed by Equation (1), EIT system could manage almost full gas volume fraction range (close 100%) for gas-water two-phase flow [9] but not the gas-water boundary.
where, c is the gas concentration, σ 1 and σ mc are the water conductivity and reconstructed local conductivity, respectively.
Commonly, grey-level or colour palette based mapping methods [10] were employed for visualising gas-water distribution by converting the different values of gas concentration to different grey levels or colours based on predefined lookup tables. However, the converted images are limited in revealing flow characteristics sufficiently, and may vary greatly in human vision and machine perception because of https://doi.org/10.1016/j.flowmeasinst.2019.101590 Received 23 April 2019; Accepted 5 July 2019 the differences from predefined lookup tables. Recently, a novel established method, called bubble mapping (BM3D) [11], enabled a good 3D visualisation of bubbles size and distribution in gas-water two-phase flow. The thresholding value used in BM3D is based on empirical knowledge, which has a potential challenge from various flow setups and condition in practice. The SPA method [12] was proposed for imaging large bubbles by determining the optimised thresholding value for processing EIT tomogram of large bubbles, but did not perform well for imaging small bubbles. This paper reports a method for 3D visualisation of gas-water flow that utilises principles of bubble mapping method and apply thresholding values determined by the size projection algorithm, providing an improvement visualisation quality on both small and large bubbles.

Methodology
Considering fully developed gas-water flows, typical flow regimes and associated flow conditions are illustrated in Fig. 1a and b in regard to horizontal and vertical flows respectively [1,13]. Except misty flow regimes, the gas distributions at the cross-section of typical flow regimes can be characterized as three categories of (1) small bubbles (e.g. bubble regime, in tails of churn regime, etc.), or (2) only a large bubble (e.g. stratified, slug, plug or annular regimes), or (3) a large bubble surrounded with few small bubbles (e.g. slug, plug or churn regimes). Further, the category (3) could also be treated as (2) since the contribution from those few small bubbles to the local gas volume fraction and also to the flow regime visualisation can be ignorable according to the general flow principle [1,14]. With the assumption of ignoring small bubbles in the category (3), the size projection algorithm is employed to visualise the distinctive large bubble as stated in flow categories (2) and (3), while the bubble mapping method is employed to visualise the small bubbles in the category (1).

Bubble mapping method
The BM3D method [11] is a new established approach aiming at enhancing the capability of EIT to visualise gas-water pipeline flows. With the input of EIT reconstructed gas concentration tomogram, BM3D could transform a stack of cross-sectional tomograms into a collection of individual gas bubbles in the pipe, which reveals a vivid 3D visualisation of disperse gas phase distribution. The BM3D method is mainly based on a predefined lookup table indexed by bubble size and an enhanced isosurface algorithm, whose major procedures are briefly introduced as follow: (1) Re-meshing: According to the spatial resolution of EIT system and the small bubble size in real situation, the pipe space is re-meshed into coarser cube cells. Thereafter, the EIT reconstructed gas concentration data are re-filled into the coarser cube cells with considering the actual gas velocity and data acquisition speed of EIT system. (2) Bubble identification: After the re-meshing, two critical parameters or thresholding values, T l (0.05, resulted from measurement error) and T g (0.4, starting forming large bubble), are employed for identifying the bubble in each cube cell. If the gas concentration is below T l , the cell is assumed being fully occupied by water. If the gas concentration is between T l and T g , a predefined lookup table transforms the reconstructed value (i.e. gas concentration in each cell) into a gas bubble whose volume fraction occupied in cube cell is equal to its gas concentration. However, if the gas concentration is above T g , the cell is assumed being fully occupied by gas, then an enhanced isosurface algorithm is employed to merge the neighboring cells with high gas concentration into a large bubble and to extract the boundary between gas bubble and water.
A set of existing concentration tomogram data collected at horizontal plug flow regime are used to evaluate the visualisation capability of the BM3D method. The plug flow was generated on a gaswater flow loop facility at the University of Leeds, where the superficial velocities of gas and water phases were 0.38 m/s and 1.02 m/s, respectively. The pipeline of flow loop is made of PVC tubes with an internal diameter of 50 mm, and the concentration data were reconstructed by a commercial EIT system with imaging speed of 312.5 fps. On one hand, four frames of concentration tomograms at different points of plug flow are processed by BM3D method, where there is a large bubble, or a large bubble with small bubbles, or tail of a plug bubble with small bubbles, or small bubbles in the pipe crosssection. The cross-sectional images obtained from BM3D method and conventional colour mapping method are compared in Fig. 2. On the other hand, 500 frames of concentration data were imported into BM3D software and processed for 3D visualisation. The visualisation results are compared with the conventional colour mapping method and the video recorded by a high-speed camera, as shown in Fig. 3.
In Fig. 2, when only a large bubble exists in the cross-section ( Fig. 2a and e), the BM3D result is similar to the image obtained from the conventional colour mapping method. However, when small bubbles are introduced, the BM3D method can approximately reveal the distribution of large and small bubbles in water, while the colour mapping method has the quite limited capability on visualisation of small bubbles. In Fig. 3, a stack of cross-sectional tomograms are transformed and displayed as a series of individual large and small bubbles (Fig. 3c) by the BM3D method. And the BM3D result is compared with the image obtained from the colour mapping method and on-site video taken through a transparent chamber by a high-speed camera. The large plug bubbles are visualised by all three methods, while small bubbles are presented by cloud in on-site video and missing in colour mapping-based image, but are highlighted in bubble mapping method. The results demonstrate the BM3D method is a good bubblebased 3D visualisation method with providing the key information, e.g. the size and shape, of both large and small bubbles, which is capable of improving flow regime visualisation and visual recognition. However, the selection of thresholding value T g (0.4) in BM3D method is based on empirical knowledge [15], which may meet a great challenging on imaging accuracy of large bubble in variation of flow setups and image reconstruction algorithms.

Size projection algorithm
The SPA method [12] is a new threshold-based image segmentation method for accurately extracting large bubbles in EIT tomogram, where the optimised thresholding value is machine-determined by a multi-step iterative process for distinguishing the gas bubble and water phase in EIT tomogram. In the SPA method, a projection error between EIT measured voltages and computed voltages of segmented image is employed and minimised for approaching the optimised thresholding value, which means the minimal projection error should be reached when the segmented image is close to the real gas distribution. Once the optimal thresholding value is reached, the large bubble is extracted by converting the original EIT concentration tomogram to binary concentration tomogram with the following Equation (2).
where, C(k) and c k ( ) are the gas concentration value of k-th pixel in the updated tomogram and the original tomogram, respectively. T sp is the optimised thresholding value determined by the SPA method for the frame of cross-sectional tomogram. The principle of SPA method is illustrated in Fig. 4.
In order to evaluate the accuracy of large bubble extraction by SPA method, a 2D gas-water two-phase model (i.e. a large bubble in water) was simulated in COMSOL software, as shown in Fig. 5a. A typical 16electrode EIT sensor and the adjacent sensing strategy were employed, which could generate 104 independent boundary voltages for solving inverse problem. Then the processes of image reconstruction and segmentation were conducted under a mesh with 1536 triangular pixels. As shown in Fig. 5, the reconstructed tomogram ( Fig. 5b) is obtained from Landweber method [16] and the binary images ( Fig. 5c and d) are obtained from threshold-based image segmentation methods with respect to the empirical thresholding value 0.4 (i.e. in BM3D method) and the optimised thresholding value (i.e. in SPA method), respectively. Two evaluation criteria [17], i.e. relative image error (IE) and correlation coefficient (CC) between the setup model and reconstructed image or segmented image, are used to estimate the imaging accuracy of large bubble.
As shown in Fig. 5, the sharp boundary of the large bubble is not clearly identified in EIT reconstructed tomogram, but explicitly extracted by the BM3D and SPA methods with their corresponding thresholding values. Meanwhile, the bubble extracted from SPA method is more accurate than it extracted from the BM3D method, which demonstrates the SPA performance of imaging large bubble is better than BM3D method. However, SPA method cannot identify the small bubbles in water since it is not sensitive to small bubbles. Therefore, SPA could only be used for imaging large bubble in each frame of EIT tomogram by the determination of the optimised thresholding value.

Combination method for 3D gas-water flow visualisation
As illustrated in Section 2.1 and 2.2, BM3D method can provide a good 3D visualisation of gas-water flow with revealing large and small bubbles distribution, while it meets a theoretical challenge on the accuracy of imaging large bubbles in various flow regimes since the  thresholding value of BM3D is a global fixed value. However, SPA method can provide the closest thresholding value for imaging large bubble in each frame of tomogram even it is not sensitive for small bubbles. According to a large scatter of data [1], large bubble starts being formed from small bubbles in pipeline when the gas fraction reaches a certain thresholding value (0.3 corresponding to the experimental conditions in this paper), that is the large bubble will exist, and in contrast almost only small bubbles exist if gas fraction is below the     thresholding value. The combination of BM3D and SPA is considered to overcome the limits in each method, where SPA is employed to determine the optimised thresholding value for BM3D method when large bubble exists. Considering the procedures of two methods, the concentration tomogram data need to be processed by SPA before being imported into BM3D software. The schematic diagram of combination procedures is depicted in Fig. 6. Since unavoidable noise in measurement, the gas concentration value might contain abnormal data in tomogram, such as negative value, and it should be filtered out and thus producing a meaningful gas concentration region, i.e. [0.0, 1.0]. Then the mean gas fraction α i is calculated as the average of gas concentration value at each pixel of a tomogram, as expressed in Equation (3), which is a decisive factor for the optimised threshold determination by employing SPA method. If the mean gas fraction of i-th frame tomogram is upper than 0.3, the frame of concentration tomogram c i is converted to a binary concentration tomogram C i by Equation (2), otherwise let = C c i i . Finally, the BM3D method is employed to transform the updated concentration tomogram data C into a 3D visualisation of gas-water flow by revealing the distribution of both large and small bubbles.
where, α i is the mean gas fraction of the i-th frame tomogram. c k ( ) i is the gas concentration value of the k-th pixel in the i-th frame tomogram. N and the subscript i are the total pixel number and the frame number of gas concentration tomogram.

Evaluation
Before the proposed combination method is evaluated, it should be clarified firstly that the objective of the study is to enhance the 3D visualisation of gas-water flow by replacing the fixed thresholding values in BM3D with optimised thresholding values determined by SPA method. The accuracy of thresholding values determined by the SPA method was demonstrated in paper [12] and Section 2.2 in this paper. The 3D gas-water flow visualisation results of proposed combination method are compared with the BM3D method, along with corresponding images from colour mapping method and high-speed camera video. For the convenience in following text, the video method, the colour mapping method, the BM3D method, and the proposed combination method (i.e. BM3D with SPA) are abbreviated as VD, CM, BM and BS, respectively. In addition, it is worth pointing out that the BM and BS visualisation results include pipe wall.

Visualisation of gas-water flow in horizontal pipeline
The horizontal gas-water flow experiments were conducted using gas-water facilities at TUV NEL (National Engineering Laboratory at Glasgow, UK). The arrangement of test section is depicted in Fig. 7. ITS V5R system [18] was employed for collecting EIT tomogram data with 312.5 dual-frames per second (dfps). Meanwhile, a high-speed camera was also installed to record the flow structures through a transparent chamber for comparison. The gas phase was nitrogen with 0 mS/cm conductivity and 1.25 kg/m 3 density, and the water phase was salty water with 33.5 mS/cm conductivity and 1049.1 kg/m 3 density. Since the involved facilities can produce gas-water flow with 0%-100% GVF, it was able to generate common flow regimes in horizontal pipeline. The superficial gas and water velocities were controlled to achieve certain flow conditions for generating typical horizontal flow regimes (including stratified, bubble, plug, slug and annular regimes). The flow conditions referenced by testing facilities are listed in Table 1.
The resultant visualisation results of the four methods are shown in Fig. 8. The VD images were generated by connecting several screenshots, and the gas and water phases in the CM images were represented by red and blue colours. For the stratified flow in Fig. 8a, the flow regime is clearly recognised by four images. Only in the BM image, few small bubbles exist at the gas-water interface. For the bubble flow in Fig. 8b, the small bubbles cannot be identified in the VD and CM images, but clearly visualised in the BM and BS images. When it comes to the plug and slug flow in Fig. 8c and d, the large bubbles can be located in all images. However, the small bubbles are missing in the VD and CM images, but clearly visualised in the BM and BS images. It is noted to point out that the proposed combination method only reveals the large bubble visualisation with eliminating the small bubbles at the existence position of a large bubble. For the annular flow in Fig. 8e, the thin water film at the top position is not revealed in the CM, BM and BS images since it is too thin to be identified by EIT system. Comparing the BM and BS visualisation results, the only difference is the elimination of small bubbles at the existence position of large bubble, because the SPA method ignores the existence of few small bubbles for achieving accurate visualisation of large bubbles in proposed combination method.

Visualisation of gas-water flow in vertical pipeline
The vertical gas-water flow experiments were performed on upward gas-water flow loop facilities with 50 mm-diameter pipeline at OLIL (Online Instrumentation laboratory) in the University of Leeds. The arrangement of test section is depicted in Fig. 9, and FICA system [8] is employed for collecting EIT tomogram data with 1000 dual-frame per second (dfps), and a high-speed camera is utilised as well. The gas phase was compressed air with 0 mS/cm conductivity and 1.29 kg/m 3 density, and the water phase was tap water with 0.35 mS/cm conductivity and 1000 kg/m 3 density. The superficial gas and water velocities were controlled to achieve certain flow conditions for generating typical vertical flow regimes (including bubble, slug and annular regimes). The flow conditions referenced by testing facilities are listed in Table 2.
The resultant visualisation results of the four methods are shown in Fig. 10. For the bubble flow in Fig. 10a, except the CM image, all images visualised the bubbles distribution by clearly showing the bubbles size and location. For the slug flow in Fig. 10b, the VD image can approximately estimate the location of slug bubble, while it can Fig. 9. The arrangement of test section at OLIL.

Table 2
The flow conditions for typical flow regimes. hardly estimate the exact bubble size and shape because of too many small bubbles. The CM image reveals the approximate shape and location of slug size, but the sharp bubble boundary and small bubbles cannot be revealed. The BM and BS images illustrate the bubble location, shape, sharp boundary and small bubbles. However, the size of slug bubble in the BS image is larger than in the BM image, and the small bubbles at the position of slug bubble are eliminated. For the annular flow in Fig. 10c, both the VD and CM images show the flow regime, but they cannot determine the size of annular air-core. However, the air-core size is clearly shown in the BM and BS images, where the size of annular air-core in the BS image is larger than the BM image, which more closely reflect the true phenomenon of annular flow regime with a thin layer of water. Comparing the BM and BS visualisation results, the proposed combination method reveals a bigger size of large bubble and eliminates the existence of small bubbles at the position of large bubble when the large bubble exists. This is because the thresholding values determined by SPA in the proposed combination method are smaller than in BM3D, and SPA method ignores the existence of few small bubbles at the large bubble position for accurately visualising large bubble.

Impact of tomographic algorithms on visualisation
The aim of EIT tomographic algorithms in gas-water two phase flow visualisation is to determine the unknown distribution of gas phase based on the measured boundary voltages. Due to the so-called "softfield" effect and ill-posed problem, the expected precision is incapable from current inverse solution. There are many tomographic algorithms developed for EIT, which can be classified into two categories, qualitative non-iterative algorithms and quantitative iterative algorithms [19]. As the fast speed is required for online gas-water flow measurement and visualisation, only qualitative non-iterative algorithms (i.e. SBP and MSBP algorithms) are discussed in this paper.
SBP algorithm was firstly produced by Kotre [20] based on the principle of linear back-projection (LBP). Later, the modified SBP algorithm (MSBP) was proposed by Wang [21] based on an approximation of inverse relation, i.e. + ≈ − x x 1 1/(1 ) at < x 1, which extends the application range further. Actually, the nonlinear approximation should be satisfied by → x 0 in mathematics. But the condition < x 1 makes MSBP algorithm work better to speed up the inverse process of imaging flow in vertical layout since the disperse phase has a more homogeneous distribution in vertical flow. It was further demonstrated that the MSBP algorithm has a better correlation of gas concentration and the conductivity ratio than the SBP algorithm [9]. As shown in Fig. 11, a comparison of SBP and MSBP results on vertical slug flow is conducted. The SBP and MSBP show a similar distribution trend in conductivity tomogram, as shown in Fig. 11a and b. But according to the comparison of gas concentration profiles in Fig. 11c, different thresholding values are needed for extracting the same bubble from SBP and MSBP images reconstructed with same measured boundary voltages, which is the impact of tomographic algorithms on the thresholding value selection for vertical flow visualisation.
However, the disperse phase in horizontal flow layout presents a highly heterogeneous, and the reconstructed conductivity change ratio might be close to 1 or even far beyond 1. Therefore, the condition of MSBP will be no longer satisfied in principle, and the reconstructed conductivity tomogram from MSBP will show few pixels having abnormal value. As shown in Fig. 12, a comparison of SBP and MSBP results on horizontal plug flow is conducted, which clarifies that the MSBP is not applicable for imaging horizontal pipeline flow. Therefore, the SBP algorithm is employed to visualising horizontal gas-water flow, while the MSBP algorithm is employed to visualising the vertical gaswater flow in this work.

Impact of thresholding values on bubble mapping
In the proposed combination method, the size projection algorithm determines the optimised thresholding values for extracting large bubble's boundary when a large bubble exists. As the mean gas fraction is less than 0.3, the SPA method will not be employed, which remains the original concentration tomogram data to form small bubbles using the BM3D method.
For the horizontal pipeline flow, the thresholding values determined by SPA method are statistically analysed, and its histograms of probability distribution are given in Fig. 13. Although the thresholding values used in the proposed method are quite different with the  thresholding values used in BM3D method, the 3D visualisation results of large bubbles from two methods are little difference, as shown in Fig. 8. That is because the employment of SBP algorithm makes the artefacts zone in reconstructed tomogram very narrow, as shown in Fig. 12c. Therefore, it is demonstrated that the thresholding values used in BM3D method for horizontal flow do not have a significant impact on the quality of visualisation. For the vertical pipeline flow, the thresholding values determined by SPA method are also statistically analysed, and its histograms of probability distribution are given in Fig. 14. Most of the thresholding values used in the proposed method are smaller than the thresholding value used in BM3D method, which results in the significant difference of large bubble visualisation in two methods. That is due to the employment of MSBP algorithm making the artefacts zone in reconstructed tomogram relatively wide, as illustrated in Fig. 11b and c. Therefore, it is demonstrated that the thresholding values used in BM3D method for vertical flow have a significant impact on the quality of visualisation.
Further, a simulation is conducted to investigate the accurate thresholding values for extracting the real bubbles with different size in MSBP reconstructed tomograms. A sequence of setups (Fig. 15a) were modeled according to the cross-sectional configuration of a large bubble in vertical flow, where the background phase is tap water (0.35 mS/cm) and the disperse phase is air (0 mS/cm). A typical 16electrode EIT sensor and adjacent sensing strategy were employed to generate the 104 independent boundary voltages for solving the inverse problem. The bubbles were reconstructed by MSBP algorithm and compared with the real setup bubbles based on the COMSOL simulation data, as shown in Fig. 15b~d, which demonstrates that the thresholding values for extracting the reconstructed bubbles with different size are different. Therefore, the fixed thresholding value (i.e. 0.4) in BM3D method tends to underestimate the size of the slug bubble or annular air-core in vertical flow. In addition, the artificial effect of blurry boundary in MSBP reconstructed results is not a narrow zone, as shown in Fig. 15b~d, which will be mistakenly converted into small bubbles at the boundary of a large bubble.

Conclusions
In this paper, a combination method utilises the size projection algorithm (SPA) to enhance the bubble mapping based 3D visualisation of gas-in-water two-phase flows. With the reasonable assumption, the large bubble's boundary is extracted by the SPA method determined optimised thresholding values, while small bubbles are formed from the original concentration tomogram data using the bubble mapping method. The evaluation results demonstrate a better visualisation performance of the proposed combination method. For horizontal gaswater flow, the thresholding value in BM3D method does not affect the visualisation very much. For vertical gas-water flow, the BM3D method tends to underestimate the large bubble size, while the proposed combination method offers a better estimation. In addition, the artificial effect of blurry boundary from conventional tomographic algorithm makes BM3D method mistakenly create small bubbles nearby a large bubble. With the employment of the SPA method, this effect can be fully removed, which are revealed in images from both vertical and horizontal flow layouts.