Evolution of High-Viscosity Gas–Liquid Flows as Viewed Through a Detrended Fluctuation Characterization

: We characterize the long-term development of high-viscosity gas–liquid intermittent ﬂows by means of a detrended ﬂuctuation analysis (DFA). To this end, the pressures measured at different locations along an ad hoc experimental ﬂow line are compared. We then analyze the relevant time-series to determine the evolution of the various kinds of intermittent ﬂow patterns associated with the mixtures under consideration. Although no pattern transitions are observed in the presence of high-viscosity mixtures, we show that the dynamical attributes of each kind of intermittence evolves from one point to another within the transport system. The analysis indicates that the loss of a long-range correlation between the pressure responses are due to the discharge processes.


Introduction
Today, the transportation of high-viscosity gas-liquid mixtures constitutes an important flow assurance problem in several oil producing countries. As an example, the production of petroleum in Mexico has declined over the last two decades, while the heavy and extra-heavy crude oil varieties account for nearly 60% of the total reserves. Although diverse technologies are available to increase the productivity (e.g., [1]), their application might not necessarily be possible, or may require to be tailored to suit the specific needs of each particular case. Therefore, understanding how certain flow characteristics evolve as the associated mixtures progress inside the pipelines is crucial to the development of adequate flow enhancement techniques.
Other problems and technical difficulties may arise in the present context. For instance, Matsubara and Naito [2] previously noted that the traditional flow maps were substantially modified with high-viscosity gas-liquid flows. The boundary lines separating intermittent flows from other patterns were drastically changed, such that at low local flow rates (or equivalently, at low local superficial velocities) the pattern belonged exclusively to the intermittent class. Similar observations were made more recently by Hernandez et al. [3] in a longer experimental flow loop.
In addition, after investigating the characteristics of the flow patterns produced with relatively high-viscosity oils and air, Foletti et al. [4] concluded that the agreement with current predictive models was rather poor. This pointed out to the necessity of conducting further work along these lines. As a matter of fact, with lower viscosity oils there was room for improvement, as the results reported by Khaledi et al. [5] clearly indicated.
Experiments carried out with oils (whose viscosities were in the range ≤10 Pa·s) also showed substantial differences between the measured data and the predicted values produced by the existing correlations. In particular, Zhao et al. [6] identified four flow patterns, three of which were intermittent. The need for improved mechanistic models that could handle such viscosities was stressed. Similar patterns were caused at still higher viscosities (∼6 Pa·s) in horizontal pipes ( [7]).
One relevant quantity affecting the properties of the intermittence is the local liquid holdup. Al-Safran et al. [8] reported a new empirical correlation that is valid in the range of 0.1-0.58 Pa·s. In a preceding article, Farsetti et al. [9] studied the frequencies produced by intermittent oil-air flows with moderate-to-high-viscosities, in addition to other flow parameters. According to the reported results, the frequency appeared to depend on the liquid holdup. Interestingly, only the estimates obtained from the ad hoc correlation proposed by Gokcal et al. [10] compared favourably with the data. The same kind of conclusions were drawn by Okezue [11], who studied high-viscosity oil-gas mixtures (from 1.1 to 4 Pa·s) over a wide range of flow conditions. More recently, Baba et al. [12] conducted additional experiments to further characterize the slug frequencies. It is noted that these flow traits have a definitive influence in the design of field operations, as well as of fluid handling facilities and equipment (e.g., [12]). A statistical analysis of the measured data was previously applied to this kind of flows by Losi et al. [13]. Basic flow features such as slug frequencies, lengths and pressure drops, were discussed in terms of the gas superficial velocities and other parameters. Furthermore, the authors derived improved correlations based on the probability density functions thus obtained.
In view of the forgoing arguments, our present aim is to highlight some of the key traits of high-viscosity gas-liquid flows by means of the well established detrended fluctuation analysis (or DFA). This alternative viewpoint has not been fully exploited in this context, and may provide a robust description of the properties of the intermittent patterns emerging naturally as the flow evolves.
Even though different methods were applied in the past to investigate low viscosity two-phase flows, only the work by Zahi et al. [14] concerns the use of DFA (to the best of our knowledge). As a matter of fact, these authors compared various fractal methods with the DFA. Within the superficial velocities of the low viscosity phases involved, it was found that the numerical value of the scaling exponent was lowest for the bubbly flow pattern, while the highest value was attained with the intermittent slug flow pattern. The former result was attributed to the random dynamics of the flow, and the latter to the periodicity of the slugging. More importantly, the DFA was shown to perform better in comparison with other fractal techniques with varying noise levels. Therefore, Zahi et al. [14] concluded that the DFA based approach can effectively underline the flow pattern characteristics. The ability of the method to provide information of the transitions was also stressed.
The possibilities and implications of the detrended fluctuation analysis were clearly illustrated by Peng et al. [15]. Because the DFA is an approach that relies on a simple scalar parameter (the scaling exponent) to represent the correlation properties of a given signal, it may be applied to extract useful dynamical information from the pressure time-series obtained at various locations. Moreover, the DFA permits the detection of long-range correlations embedded in seemingly non-stationary time-series, while avoiding spurious artifacts related with stationary effects. This technique has been applied to such diverse fields of interest as DNA, heart rate dynamics, neuron spiking, human gait, long-time weather records, cloud structure, economical time-series, solid state physics, and even reservoir characterization (e.g., [16][17][18]).

Test Apparatus
The experiments were conducted in the flow loop shown in Figure 1. This setup was designed to promote the long-term evolution of the mixture, since the pipe's length-to-diameter ratio was of order L/d ∼ 10 2 . The test section was constructed with (schedule 80) steel tubes with an internal diameter of 0.0762 m (3 in) and overall length of 54 m. It should be mentioned that a single 3 m section, made of transparent PVC tube, was fitted 3 m upstream from the discharge point for visualization purposes.
Glycerin was selected as the Newtonian, high-viscosity, liquid phase. Its measured dynamic viscosity was η = 1.1 Pa·s (1100 cP) at a room temperature of 25 • C. The liquid phase was supplied to the test section by a progressive cavity pump (Seepex Mod. BN35-24). Regardless of the pressure buildup at its discharge plane, the pump was capable of delivering constant mass flow rates in the interval [0.0 kg/s, 6.1 kg/s]. The outlet of the test section was connected to a separator tank with an internal capacity of 1.5 m 3 . The pressure in its interior was kept constant. On the other hand, a twin-scroll Kaeser Aircenter SK.2 compressor supplied a constant mass of dry air (at room temperature) for the pressures within the interval [0.0 Pa, 1.6 × 10 6 Pa]. The mass flow rate was finely tuned with a regulator and globe valves to produce a choked condition at the injection port. The mixture was produced at the 3-way connection tube shown in Figure 1.
Both mass flow rates were measured at the inlet with Endress-Hauser Coriolis meters. The Promass 83F80 DN80 3" model used with the glycerine had an accuracy of 0.1% across the entire measuring range. Similarly, the Promass 83F50 DN50 2" model used with the air had an accuracy of 0.05% across the entire measuring range. Also, all the pressures were measured with an array of conventional MEAS U5300 transducers with a resolution of 0.1% across their measuring ranges. The two measuring intervals selected for these instruments were [0.0 Pa, 1.03 × 10 5 Pa], and [0.0 Pa, 3.45 × 10 5 Pa]. Table 1 summarizes the reference values of the fluid properties. Table 1. Fluid properties at standard conditions (i.e., @ P = 1.013 × 10 5 Pa, T = 25 • C).

Fluid Viscosity Density Interfacial Tension (Pa·s) (kg/m 3 ) (N/m)
Impedance computed tomography scans (or ICT scans) provided snapshots of the internal flow configuration. The tomograms were produced with the International Tomography Systems (ITS) ERT and ECT instruments. Each tomogram is constituted by snapshots taken at a sampling rate of 30 frames per second in two separate measuring planes. Then, the computerized system reconstructs the flow structure from the collected images. Although the ECT-ERT systems also yields an accurate measurement of the holdup within the measuring planes, the corresponding times series are similar to those of the pressure. These measurements will be discussed in a future paper.

Methodology
The experiments were produced in accordance with the inlet mass flow rates indicated in Table 2. In total, 15 different combination pairs (q g , q l ) were considered, and 10 experiments were conducted for each one of them. Table 2. Experimental matrix. Mass flow rates of the liquid (q l ) and the gas (q g ) phases injected into the test section. Every experiment conformed to the following procedure: first, a steady state, single phase flow of glycerin was produced at a given flow rate. Owing to the elevated viscosity, the upper bound on the Reynolds number (Re l = (ρud/η) l ∼ 10 2 ) implied that the Poiseuille velocity profile developed fully in just under 20 pipe diameters. The steady state condition was then verified by making sure that ∆p/∆t ≈ 0 in all pressure transducers. Once this condition was reached, a specific mass flow rate of air was injected (see Table 2). Because Re l /Re g ∼ 10 −2 -10 −4 , the head loss was essentially determined by the liquid phase; accordingly, the two-phase (or mixture) Reynolds number was simply Re TP ∼ Re l . Even though the two-phase flow developed rather quickly, the system was nevertheless allowed to settle for at least 2 minutes before taking any kind of measurement. At this point the flow properties were measured and collected by the data acquisition system.
In order to remove any external influence that might affect the measurements, only the differential pressures were used to determine the pressure gradients of interest. The measuring ports were located at 0, 18, 22, and 43 m, downstream from the inlet. These ports are labeled with P1, P2, P3 and P4, respectively, in Figure 1.
The experiments were conducted at an average ambient temperature of 25 • C, as measured inside the laboratory. In order to verify that the temperature of the fluids in the flow loop did not vary as the experiments were conducted, two Type K thermocouples were installed at the inlet and outlet sections. It is worth mentioning that no temperature fluctuations were observed within their resolution limit (0.5 • C across the measuring range). Nevertheless, the mixture's viscosity was checked for each test by taking samples (from the previously indicated locations) and by measuring their viscosities with a Brookfield DV2T viscometer.
Following Kline-McClintok's method to determine the experimental uncertainties, it was readily verified that the pressure time-series were measured with an uncertainty of ±10% with a confidence level of 95% ( [19,20]). It is worth mentioning that this uncertainty is well within the normally accepted range in most multiphase flow contexts (e.g., [21][22][23]). Accordingly, at least 3 experiments were conducted for each entry of the experimental matrix (i.e., for every (q g , q l ) combination). Around 900 samples were gathered for the measured variables per test. The universal set contained nearly 750,000 data points comprising the simultaneous measurement of mass flow rates, pressures and temperatures. In general terms, it was observed that the statistical distribution of the data points conformed to Student's t-distribution with a 95% confidence level.

Results and Discussion
The set of electric impedance tomograms shown in Figure 2 portray the internal phase distribution for three different flow regimes (in correspondence with the low, intermediate and high liquid flow rates). Each image constitutes a snapshots of the actual structure by the two phases as they traverse the measuring device.  Table 2. From top to bottom, image (a) shows the flow structure for q g > q l , while (b) shows it for q g ∼ q l , and (c) for q g < q l .
Respectively, image (a) illustrates the classical structure of an elongated bubble pattern resulting when the relatively low liquid flow rate, q l = 1.3 kg/s, is introduced in the pipe. By increasing q l to 3.7 kg/s the bubbles become shorter. This leads, in turn, to an augmented frequency. Similarly, by further rising the liquid flow rate to q l = 6.1 kg/s the number of smaller bubbles and the corresponding frequency increase considerably.

Pressure Response
The present analysis focuses on the pressure signals alone, because the corresponding time-series encode information of every effect related with the flow. This approach proves to be convenient in field applications ( [24,25]). For illustration purposes, Figure 3 shows a typical set of pressure measurements obtained at the U-turn section. The curves are deliberately plotted with different scales, in order to highlight their distinctive characteristics for all the flow rate combinations given in Table 2.
Clearly, to every combination (q g , q l ) corresponds a particular fluctuation pattern. From left to right the fist image shows how, with low liquid flow rates, the resulting fluctuations tend to be relatively large with respect to the average pressure. The converse is true when more liquid is pumped into the system because the mean pressure increases substantially. In general terms, these non-periodic amplitudes become quite important as they account for significant over-pressures inside the pipe. Obvious external manifestations of such dynamical effects may include vibrations, as well as a very noticeable recoil of the pipe during the ejection of the liquid slugs.
In contrast to the preceding behavior, the fluctuation amplitudes seem to be about the same order of magnitude when the gas mass flow rates increase. However, the resulting frequencies appear to increase with q g . One may also notice that the intermittence tends to become more regular. This is easily explained in view of the fact that: (a) many more bubbles are formed, and (b) these tend to be better distributed throughout the entire length of the pipe. This latter effect is purely conditioned by the elevated viscosity of the liquid phase.
Interestingly, as will be shown, these broad conclusions may not be directly extended to other regions of the flow system. Instead of analyzing the time-series in those regions, the data is subjected to the Detrended Fluctuation Analysis. In preparation for that, Figures 4-6 depict the pressure differences near the inlet, at the mid-section and near the outlet of the pipe. No information has been lost to this process, because the data has not been filtered. These ∆P are calculated at the locations specified previously, while the flow rate combinations are the same in all figures.

Spectral Evolution of the Flow
In order to extract useful information from the pressure time-series, it proves convenient to first analyze the signals in the frequency domain. Since the primary interest of the analysis concerns the pressure fluctuations, the mean pressure is subtracted from the time-series. The set of images in Figures 7-9 show the spectrograms for the flow rates given in the experimental matrix ( Table 2). Note that the three sets correspond to the three different locations along the pipe. Again, the columns represent q l and the rows q g (their values are given in each snapshot).
These plots provide information about the frequency bandwidths that are excited at specific times for different inlet flow rates (q g , q l ). Moreover, since the three sets correspond to the three specified locations along the tube, they also convey important information on the evolution of the intermittence. The color scale indicates the amplitude of the excited modes in a given frequency neighbourhood. Thus, for example, in the upper left image of Figure 7 one may observe that modes with frequencies in the interval (0.02 Hz ≤ f ≤ 0.07 Hz) are excited from 60 s to 120 s. This relatively long time spans most of the experiment. Notice that a strong peak, centered at 0.05 Hz, is produced at approximately t = 80 s after the experiment is initiated. In contrast, for (q g , q l ) = (0.005 kg/s, 3.7 kg/s) and (q g , q l ) = 0.005 kg/s, 6.1 kg/s) only a very narrow, low-frequency band is excited towards the end of the experiment (t ≈ 140 s).
One may conclude that, in a broad sense, the periodic fluctuations appear as high intensity (red) spots in the spectrograms. Conversely, noisy flow fields are represented by a rather diffuse, low intensity, distribution of interconnected bands with a color scheme much closer to the blue. Nevertheless, the abundance of spectral realizations induced by most (q g , q l ) combinations is clearly exemplified by the amount of modes excited in distinct bandwidths, and at different times. For the sake of comparison consider the images corresponding to (q g , q l ) = (0.015 kg/s, 3.7 kg/s) and (q g , q l ) = (0.015 kg/s, 3.7 kg/s), i.e., two neighbouring flow rate combinations of the experimental matrix, which clearly illustrate the complexity of the fluctuation patterns.
Perhaps the most interesting feature portrayed by these images is the particular space-wise evolution of the fluctuation patterns. Consider in this case the spectrograms for (q g , q l ) = (0.005 kg/s, 1.3 kg/s) in Figures 7-9. What initially constitutes a single, narrow frequency band extending in time (Figure 7), eventually unfolds into a couple of short frequency bands excited at the same time (seen as two red spots at 60 s ≤ t ≤ 120 s in Figure 8), plus a wide range of less energetic modes, in the mid-section of the tube. Finally, on the downstream section of the tube, two similar frequency bands are excited at two different times (i.e., The two red blobs centered at f ≈ 0.05 Hz now appear at t ∼ 75 s and t ∼ 125 s in Figure 9). Notice that the less energetic modes have already decayed at this final stage. Interestingly, other flow combinations seem to be relatively immune to such changes and, therefore, tend to persist in time.

Detrended Fluctuation Analysis (DFA)
The characteristics of the spectral evolution, and consequently of the fow intermittence, may be better understood in terms of the detrended fluctuation analysis (DFA). The reason is that two-phase flows exhibit complex self-similar fluctuations over a broad range of space and time scales, much in the same way as in many physical, biological, physiological and economic systems, where multiple component feedback interactions play a central role.
Because of the nonlinear mechanisms controlling the underlying interactions, the output signals of this two-phase flow are typically non-stationary. Furthermore, they are characterized by embedded trends in heterogeneous segment patches with different (local) statistical properties. In this case, traditional methods, such as the spectral and the auto-correlation analyses, may not be necessarily well suited to study this kind of processes. The DFA methodology, in contrast, enables a proper quantification of long-range correlations with non-stationary, fluctuating signals.
The procedure relies on the profiling of the pressure, P, through an 'integration' of the time-series where p is the mean value of the signal and p(i) is the value of the pressure at time i. First, the initial data set is divided into equal segments of length n. Then, a linear approximation P n is computed with a least squares fit performed over each separate segment. The resulting fit represents the trend in the given section. Next, the average fluctuation δp(n) of the signal about the trend is determined with With these elements one may proceed to directly plot log δp(n) as a function of log n. The scaling exponent α then corresponds to the slope of the least squares, linear regression fit to the data points Accordingly, log δp(n) ∼ α log n, and the numerical value of the scaling exponent (also known as the fractal scale, autocorrelation exponent, or self-similarity parameter) is given by As in previous cases the overline indicates the mean value of the respective quantity. Power law relationships of the kind represented by Equation (3) suggest the presence of self-similar fluctuations. Furthermore, depending on the actual value of α, the following cases may be considered: This DFA analysis scheme produced the results summarized in Table 3 for the experimental data sets obtained in the laboratory. It is readily seen that the table is divided in three sections, one for each of the pressure differences p 1 − p 2 , p 2 − p 3 and p 3 − p 4 previously discussed. The rows and columns are in full correspondence with the liquid and gas flow rates of the experimental matrix. All flow rate combinations are considered. Table 3. Numerical values of α. These values provide the DFA footprint of the flows generated with the inlet flow rates given in the experimental matrix (see Table 2). The numerical values acquired by the scaling exponent allows an interpretation of certain aspects stemming from the long range development of the flow. Thus, with high-viscosity gas-liquid mixtures we may draw two general observations: First, it is noticed that the flow in the exit section of the test section is much noisier than in the inlet section. This effect is caused by the ejection of the irregular liquid slugs at the outlet of the pipe. Accordingly, in the long-term development of the flow the autocorrelated (α > 0.5) and anti-correlated (α < 0.5) patterns tend to become noisier (i.e., as white, pink or Brownian noises). This means that the slugs' irregularities increase as the flow progresses towards the outlet. In contrast, the U-turn renders the flow less noisy, possibly due to a regularization effect induced by the secondary flows.
Secondly, in relation with the frames corresponding to (q g , q l ) = (0.005 kg/s, 3.7 kg/s) and (q g , q l ) = (0.005 kg/s, 6.1 kg/s) in figs. 8 and 10, one may notice two characteristic time scales: one corresponds to the flow time scale, while the other to a low frequency mode with a much larger time scale (easily seen as a small red spot appearing at the bottom of the images at t = 150 s). Nonetheless, at higher gas flow rates, the flow in the exit section tends to be more regular, because higher modes are excited (i.e., intense red spots).
For concreteness, however, the flow properties may be described in more detail as follows: 1.
The pressure drop exhibits a space-wise evolution. It may be caused by local variations of the phase volume fractions, which take place as the mixture evolves downstream along the pipe. For example, consider the flow produced with (q g , q l ) = (0.005 kg/s, 1.3 kg/s) (upper most squares of columns 1, 4 and 7 of the table). The pressures appear to be autocorrelated except at the downstream section.
In other words, the pressure measured at P1 corresponds to the pressures measured at P2 and P3, the three being in-phase because α > 0.5 (positive autocorrelation). On the downstream section, however, the pressure measured at P4 does not correspond to any of the measurements at P1, P2 or P3. Hence, three possibilities arise: a) the coherence may be lost to the strong effects induced by the ejection of irregular liquid lumps at the outlet, b) secondary flows at the U-turn may have a disruptive effect on the properties of the pressure waves, and c) a combination of the two. The reason for these possibilities is that the U-turn and the outlet are the only two points in the flow system where the two following effects can take place: secondary flows due to the centripetal acceleration undergone by the flow inside the U-turn section, and the sudden depressurization caused by the ejection of the liquid slugs into the separator tank.

2.
Apparently, very few flows exhibit white noise characteristics or non-stationary fluctuations. The former are more related with higher liquid and gas flow rates, while the latter are more related with low q g and high liquid q l .
The non-stationary cases represented by α > 1 are interesting, because they indicate that the mean values of the pressures are varying with time (or equivalently with position). Even though these variations might be slight, the method is still capable of identifying them. This opens up the prospect of designing techniques for industrial applications. Obvious examples would be the development of methods to detect small leakages in pipelines, or to detect slow corrosion processes.

3.
It is worth noticing that pink noise processes are mostly observed at the outlet section of the pipe. In general, the α values of this section appear to show a relative increase with respect to the values of the inlet section. This suggests that the scaling exponent increases in the direction of the flow.
Since pink noise refers to scalability, the reproduction of self-similar patterns and small scale traits of the signal would suggest the existence of an energy distribution process (analogous to the energy cascade in turbulence). However, it is noted that only the cases corresponding to the inlet flow rates (q g , q l ) = (0.005 kg/s, 4.9 kg/s) and (q g , q l ) = (0.01 kg/s, 6.1 kg/s) maintain this behavior. On the other hand, only the flow combination (q g , q l ) = (0.005 kg/s, 6.1 kg/s) seems to correspond to this process in the U-turn. Overall, from the physical point of view, these characteristics would be mostly related with the ejection effects produced at the outlet of the pipe.

4.
Interestingly, not a single combination of inlet gas-liquid flow rates produced fully random processes in this kind of flow system. The question still remains whether such Brownian noise patterns would eventually emerge in a longer pipeline, or not. The same question could be asked regarding the flow rates.
Anti-correlated flows for which α < 0.5 seem to dominate in the U-turn section of the pipe. This is particularly true with high flow regimes, that is, those produced with elevated inlet mass flow rates. Similar anti-correlations are also observed at the inlet and outlet section of the pipe with high flow rates. These cases indicate that the phases of the pressure waves shift by approximately π radians, as they progress from one pressure port to the next one.

Conclusions
Experiments were carried out with high vicosity mixtures of glycerine and air, in a flow loop with a length to diameter ratio of l/d ∼ 10 2 . The liquid phase had a viscosity of 1.1 Pa·s.
A detrended fluctuation analysis, or DFA, was applied to the experimental measurements obtained in this facility. Concretely, the technique was applied to the pressure time-series, under the assumption that relevant information can be extracted from them in a meaningful manner. It is important to note that the mentioned experiments are significantly different to those so far reported in the open literature. Therefore, a longer term development of the flow is possible, and a variety of dynamical effects are reflected in the spectral content of the registered pressure signals.
The DFA showed that each of the 15 mixture combinations formed by different fractions of air and glycerin have a unique head loss pattern, which can be expressed by power law expressions. In principle, this uniqueness can be exploited by pattern recognition algorithms for the identification of flow regimes, which is an important stage to design separation equipment, slug catchers, gas lift operations, wellhead gathering systems, and production management.
It was noticed that the frequencies are relatively low as compared with those measured with low viscosity water-air or nitrogen-water mixtures. The primary factor producing this outcome is the viscosity of the glycerine.
The DFA analysis enabled the proper identification of various energy states of the flow, in accordance with the observed behavior of the flow pattern at different locations. The corresponding states were determined from the variations of the scaling exponent, which was measured at equally spaced measuring ports on the flow loop.
Author Contributions: All authors contributed equally to this work.