FLOW VELOCITY MEASUREMENT METHODS USING ELECTRICAL CAPACITANCE TOMOGRAPHY

The current methods of calculating the velocity field based on a twin-plane electrical capacitance system for multiphase flows are presented. The described methods allow to calculate the velocity profile of the multiphase flow in cross-section of industrial installation. The theoretical assumptions of the considered methods are also noticed. The main advantages and disadvantages of the authors’ flow velocity measurement methods are discussed.


Introduction
Modern flowmeters enable us to evaluate how much material has passed through the meter in a certain time (for example kilograms per hour).For single flows, the technology to measure the flow is well understood and flowmeters are widespread.However, many industrial processes use multiphase flows pumped along pipes.The measurement of the multiphase flow of gas, oil, and water is not trivial [2,19,20,23] and significant challenges include:  the reduction of multiphase flow measurement uncertainties which depends on the measurement conditions;  making the flow meters less sensitive to changes in the physical properties of the flow components;  the improvement of stability over time;  the reduction of the flow regime dependency for the multiphase flow meters.Electrical Tomography (ET) is a typical completely nonintrusive technique for measuring multiphase flow in pipes [19,18,22].This technique can be used for measuring different kinds of multiphase flows such as gas/solids [11], gas/liquid [17], and liquid/solids [23] flows.
There exist two widely applied types of ET tomography: Electrical Resistance Tomography (ERT) and Electrical Capacitance Tomography (ECT).Currently more papers about flow measurement using ECT [20] and ERT systems [5] have been published.
The flowmeter for measuring the flow rate based on Electrical Capacitance Tomography technique consists of two separate planes of sensors, which measure electrical capacitance between two pairs of electrodes around a pipe.
To calculate the cross sectional distribution of the permittivity of the material (or concentration profiles) inside the pipe from the capacitance measurements an image reconstruction algorithm is used.The images provide instantaneous information on material distribution in the two pipe cross-sections at given locations.With the sensor dimensions we can then calculate the volume flow rate per zone.As result, the integration of the volumetric flow, if the density is known, mass flow can be derived.
Although acquisition and processing were initially seriously limited by computing speed, the commercial ECT systems, which are now available, can produce twin-plane images with real-time image capture and reconstruction at up to 500 complete frames per second.Such systems include major commercial applications such as pneumatic conveying of solids, many flow mixtures within the oil industry, and gas-solid fluidized beds [6,20].
The basic parameter for flow measurement is a velocity.The article presents the review of the authors' methods, which permit the measurement of a flow velocity more accurately and reliably.

Methods for flow velocity calculation in Tomography image Velocimetries -review
The solids motion has a stochastic character and hence typical solid concentration changes within any pixel of two cross-sections can be treated as stochastic stationary process.
Let x n (iT) and y n (iT), n=0,...,N-1, where N is the total number of image pixels, define instantaneous distributions of material for the nth pixel of the ith image obtained from the cross-sections X and Y respectively within the frame rate resolution T (frame per second).Nowadays, there exist different methods to calculate the transit time τ of flow propagation from reconstructed images of two planes.Most of them are based on cross-correlation function or spectral density of instantaneous distributions of material for pixels of the reconstructed images [2].
These methods for calculating velocity are based on two main assumptions: 1) The trajectories of motion between these sensors are parallel to the pipe axis (the trajectories are perpendicular to the sensor plane), 2) The changes within any pixel of a pipe cross-section are treated as a stationary or quasi-stationary process and hence the correlation function R xy (τ) describes the general dependence between any process from the first plane and another from the second plane.As discussed in the literature [21], the particle transition time τ can be found by cross-correlating corresponding processes x n (t) and y n (t).For the reconstructed images obtained within the frame rate T the cross correlation can be calculated as following: where p = ..., -1, 0, 1, ... is the shift time in the number of frames; n is the pixel index, M is the number of images for which the cross-correlation is calculated; x n (iT), y n (iT) are the numerical values, associated with pixel n from i image obtained from sensors of a plane X and plane Y, respectively, T is time period of a single frame.
In the simplest assumption the shift p  at the peak of the correlogram corresponds to a transit time τ of flow structures between two planes, τ x = p  .•T.Knowledge of the transit time τ and distance d between the two sensors enables the solids velocity within nth pixel to be calculated from the following equation (2): In total, N velocities, which have created a velocity profile, can be calculated [v n ], n=0,…,N-1.
The calculated correlogram has a clearly discernible peak if the flow patterns take a form similar a rectangular impulse (such as flow instabilities referred to as "slugs" and "plugs" in the dense pneumatic conveying transport) over the correlation interval.Consequently, the flow regime (stratified, slug, annular or bubble flow and flow direction) has significant influence for finding a correlogram peak.
The peak may be found by the greatest single value, centre of area or polynomial fitting.The polynomial fitting gives the most consistent results though all the other techniques are available in the software.
Although cross-correlation technique has many advantages, there are, however, many disadvantages.In practice, the following factors have a great influence on finding a peak of crosscorrelation function.The correlation is essentially the square of the signal, so within any volume the largest signal changes, and hence the largest flow structures, dominate the result.The time window must be determined in such a way that correlation function calculated of the flow pattern over the window contains an evident peak.To skip the limitations of classical approach based on the correlation function for calculation of a flow velocity, several methods were proposed in last decades by the authors.

The "best-correlated pixels" method [15]
As been noted, the cross-correlation technique for flow velocity measurements is based on the assumption that solids trajectories are parallel and perpendicular to the sensor plane.Therefore, such a method is flow pattern dependent.The method discussed below calculates the flow velocity, that does not require the assumption about the laminar trajectory within a sensor volume.This is important since in many cases trajectories exhibit very complex shapes.
For example, a dense pneumatic conveying is realised by the occurrence of "slugs" and "plugs", which essentially should be analysed as wave propagation phenomena.In this situation, finding a link between actual solids trajectories and the movement of flow instabilities is not straightforward.This constitutes a serious limitation for applying direct cross-correlation techniques.Another limitation factor is that the ECT system can only measure the average solids concentration in a rather large control volume (resulting from both the spatial resolution and the electrode length) compared to the particle size.Consequently, it is not always possible to distinguish between stationary and moving material in the pipe.Transport in the vertical pipe can be in both directions: upward and downward.In vertical conveying there is the added force of gravity acting against the flow of material causing the plugs to compress further, which partially accounts for the concentration being higher than in a horizontal flow.There are also particles 'raining' down from the back of the previous plug, applying stress to the nose of the next plug.
All these factors mean that the flow is usually very dynamic and non-laminar.Classical methods based on correlating corresponding pixels are not suitable for studying such phenomena.There is also the need to improve precision of flow measurements by respecting its dynamic nature in calculations.
The proposed 'best correlated pixels" method uses a correlation technique in a new way.Figure 1 shows covariant matrices of a sequence of frames with chosen pixels.These matrices are presented as images where the pixel colour is equal to a value covariant function between the chosen pixel and all other pixels.As observed in Fig. 1, the correlation order of pixel neighbourhood is large.In the classical method only corresponding pixels are correlated, but the hypothesis can be advanced by proposing that a signal from one pixel on the first plane is sometimes better correlated with a signal from a noncorresponding pixel on the second plane than with a corresponding pixel.So the idea is to calculate cross-correlation between a pixel from the first plane and some pixels from the second plane.The pixels from the second plane are chosen from the corresponding pixel and its neighbours.The best-correlated pair is taken from the correlation calculations.The idea is shown in Figure 2.
Figure 2 shows cross-section X and the corresponding crosssection Y.The arrows denote the correlation between the pixel X (n,m) and the neighbourhood of the pixel Y (n,m) .The bestcorrelated pair is chosen.
Equation 2 with some modification enables the calculation of the correlation function between a pixel from plane X with several pixels from Y plane as shown in equation 3 below: where (n,m) are the coordinates of the pixel; are the value of the pixel (n,m) on planes X and Y; B is the neighbourhood of pixel (n,m) on plane Y; B  Y.
The average value of maximum correlation function for bestcorrelated pixels is about 5% higher than the maximum correlation function obtained for corresponding pixels.The maximum value of the correlation function is usually about 0.7-0.9 for the normalised correlation function.
The different correlation functions often mean different transit time.In the example shown in Figure 3 the transit time between the corresponding pixels is estimated as 45 ms.The transition time between the best correlated pair is estimated as 65 ms.Thus the solids velocity at pixel X (15,15) is calculated from transit time 45 ms the result is 2.9 m/s, but if the second transition time is used (65ms) the velocity at the same pixel is 2.0 m/s.The difference is quite significant.(15,15) on plane X according to co-ordinates on plane Y: per-to-per- (15,15); best neighbour - (11,13).One frame is equal to 0.01 s If there is a difference between velocity values obtained from corresponding and best correlated pixels there is also a difference between mass flow rate calculated using those two methods of obtaining velocity profile.In one of the experiments mass flow rate equal to 725 kg/h was obtained for the classic method and 736 kg/h for the best-correlated pixels method.The solids flow rate obtained by the weighing method, as reported in [14], was 789 kg/h.

Flow velocity measurement based on a virtual channel concept [7]
From the viewpoint of using an ECT, the following parameters can be used to characterise a flow:  a profile velocity,  an average concentration calculated inside each pixel.
Examples of these characteristics are presented in Fig. 4. As can be seen, these flow parameters are not the same in different cross-section regions.As is shown in the figure, the particles have high and low velocities in different regions of the cross-section (see Fig. 4a).The similar case can be observed for an average concentration profile (see Fig. 4b).The investigation carried out allows the conclusion to be made that flow presents itself as a set of independent variables, which present a certain flow structure changing with time.
Thus, flow structure detection allows one to analyse the behaviour of variables of the separated structures as an alternative to the whole flow analysis.Such variables can be observed as certain spaces between two planes of cross-sections called virtual channels, where the motion of particles has similar properties.Fig. 5 shows an example of such a flow structure as a set of channels.Such analysis allows the behaviour of structures of a flowing medium to be observed and changes in their positions to be defined.In this work, we consider a simplified case of the virtual channel concept where each virtual channel is represented by the connection of identical regions of two planes (see Fig. 5b).These regions are represented by an image called an image of regions, and the flow structure is unequivocally defined by this image (see Fig. 5c).

a) virtual channels concept; b) simplified virtual channels; c: cross-section splitting
The splitting of the cross-section can then be performed periodically, and as a result one obtains sequences of images of regions in which the chosen region corresponds to the virtual channels of the flow.
A stochastic model [15] is used for the detection of virtual channels.In a given model a flow X(t) is defined by an ndimensional vector X(t)=[X i (t)] where i=0,..,N-1, N -the number of pixel elements.The vector element X i (t) determines probability change in dielectric permittivity (concentration) inside an i-th pixel, called an elementary process (EP).In notions of this model splitting can be defined as the grouping of the EP processes for both cross-sections into some groups so that their characteristics have similar properties due to chosen homogeneity criterion.In our case a transit time of solids between two planes was selected as criterion.
To calculate the transit time {t i } i=0,...,N-1 a classical method was used based on finding the maximum value of a crosscorrelation function.As mentioned above there exist probabilistic relationships between the separated velocities {v i } in some regions.For the grouping of EP were proposed to create k possible different groups of obtained set {t i } by "k-means" clustering algorithm.Such grouping is aimed to split EP processes corresponding to the homogeneous elements of a set {t i }.
The analysis was performed for tomographic images obtained from a pneumatic conveying flow rig at the Department of Chemical Engineering, UMIST.The details of the system are described by [11].All the researches were carried out for reconstructed tomographic images of a vertical pipe of a dense pneumatic conveying system.A sequence of a number of plugs appeared, on average, every 3 to 4 seconds.At the end of each sequence some granular material (most likely from the tail of the "train") dropped downwards under the gravity.Firstly, the macroscopic transport in the vertical pipe can take place in both directions: upward and downward.Secondly, in the micro scale, particle transport can be downward ("raining" particles), while, at the same time, plugs in the macroscopic sense travel upwards.The results presented below correspond to the reference gas velocity of 2.2 m s-1 (for an empty pipe).
The proposed k-means clustering algorithm [3] has performed a non-hierarchical divisive group analysis on the input data: a set of EP processes and a set of {t i }.The number of group (group elementary processes with similar transit time) is an input to the kmeans clustering algorithm.Groups are described by centroids, which are the point whose parameter values -transit time -are the average of a transit time t i of all elementary processes in the chosen group.It is also known as a representative element.The similarity between EP i and EP j processes are measured by using Euclidean distance squared between time t i and t j .The main goal of the k-means algorithm is to minimise dissimilarity in elementary processes within each group while maximising this value between elementary processes in a different group.Each EP process is assigned to the centroid (and hence group) with the closest Euclidean distance.New centroids of the k groups are computed after all objects have been assigned.The steps of assigning objects to centroids and computing new centroids are repeated until no objects are moved between the groups.
This method of grouping of elementary processes is conceptually simple: 1) Select k initial cluster centroids, c 1 , c 2 , c 3 , ..., c k .2) Assign each instance t i into the group whose centroid is the nearest to t i .3) For each group, recompute its centroid based on elements ti contained in this group.4) Go to (2) until convergence is achieved.
As already mentioned notion "virtual channel" is represented by group of EP processes.Once the virtual channels have been established a series of parameters like the mean concentration, the velocity, mass flow rate.
The mean concentration k C is calculated as a mean value from the concentration in all pixels, which belong to the region of respective channel: IAPGOŚ 1/2017 p-ISSN 2083-0157, e-ISSN 2391-6761 where C i -the material concentration in pixel i, I k -the set of indexes of pixels which belong to the region k, N k -the number of pixels in region k.
In a similar way the velocity inside a channel k V can be calculated: where V i -the velocity in pixel i, which belong to the region k of respective channel The mass flow rate for k channel k F can be calculated in two ways.Firstly, (let us call this a classical method) the velocity and concentration profiles can be used: where 3600 -the multiplier derived from the number of seconds in one hour; S -the area of a single pixel; Q -the bulk density of the material.
Another way to calculate the mass flow rate k F using the concentration k C and the velocity The virtual channel concept has verified for different behaviour of the solid particles in a vertical pipe flow:  determined virtual channels during a plug -(material move upwards along the pipe)  determined virtual channels at a tail of plug (some material dropped downwards under the gravity).Fig. 6 presents the example profiles of velocity, concentration, mass flow rate and image of regions for tail of plug.The selected channels parameters are given in table 1.

Weighted-mean phases method [8]
Although cross-correlation technique has many advantages, there are, however, many disadvantages.In practice, the following factors have a great influence on finding a peak of crosscorrelation function.The correlation is essentially the square of the signal, so within any volume the largest signal changes, and hence the largest flow structures, dominate the result.
The time window must be determined in such a way that correlation function calculated of the flow pattern over the window contains an evident peak.
In practice, the typical concentration changes within a pixel for the two planes as the function of a frame number has no evident peak, and thus using one to determine the transit time is not easy.Hence, classical methods calculating a flow velocity based on the correlation of the corresponding pixels are not suitable for certain flow-patterns and finding the peak of a correlogram is ambiguous.
These limitations seriously restrict the application of crosscorrelation in flow measurement and the range of operating conditions where reliable answers can be obtained.The weightedmean phases method (Mosorov 2006) overcomes these limitations permitting the flow velocity to be unambiguously calculated with high robust.It is based on a spectrum analysis of the pixel values of both cross-sections of the tomographic system.
In the case of the weighted-mean phases method, a transit time τ x is calculated as the weighted-mean value of the time translations  k of the same k-harmonics of x n and y n signals: where A k x and A k x are the amplitudes and φ k x and φ k y are the phases of the kth spectral components S xn (jkω 0 ), S yn (jkω 0 ) of the x n and y n respectively.
The proposed method can easily be applied to the existing flowmeters based on twin plane tomographic system by upgrading the software.The one has no limitation on the application for other types of modality such as gamma or optical tomography.
The verification of the proposed method that velocity calculated by the proposed method is a little more accurate compared with theoretical velocity (around 5%) than one calculated by the correlation method [2].The probable reason that the phase method gives better results than the correlation method is that correlation is essentially the square of the signal and hence the largest signal changes dominate the result; therefore the dispersion of these calculated values is small.Otherwise the A weighted-mean phases method is more sensitive even for smaller signal changes.This is probably also the reason that the dispersion of velocity values by the proposed method is higher in comparison with those derived by the correlation method.The results using the pneumatic conveying flow show that the proposed estimator achieves a bit more accuracy.Finally, the proposed estimator is relatively fast in software implementation because two independent FFT transformations are needed instead of correlation function calculation.

Conclusion
The described methods enable to obtain the velocity vector (three components vector), not only one component as in classical approach.There are based on the concept of "best correlated pixels", where the pair of pixels has the largest values of the crosscorrelation against the other pairs.Contrary to the classical method, this concept does not require the additional assumption concerned with the solids direction within the sensor volume.This is crucial since in many industrially important cases solids trajectories have a complex three-dimensional character.
This method, combined with electrical capacitance tomography, enables determine the particle velocity field, in gas/solids flow.However, it can be combined with other tomographic modalities.The presented results were obtained both for gas/solids flows in vertical and horizontal channels.Future application of this technique for a solids mass flow measurement will help to develop a more accurate solids mass flow meter and allow to analyse, with high accuracy, swirl flow.

Fig. 2 .Fig. 3 .
Fig. 2. The principle of the cross correlation between the pixel Xn,m and the neighbourhood of the pixel Yn,m

Table 1 .
Selected channels parameters