## Abstract

Electromagnetic flow meters (EMFMs) are the gold standard in measuring flow velocity in process industry. The flow meters can measure the mean flow velocity of conductive liquids and slurries. A drawback of this approach is that the velocity field cannot be determined. Asymmetric axial flows, often encountered in multiphase flows, pipe elbows and T-junctions, are problematic and can lead to serious systematic errors. Recently, electromagnetic flow tomography (EMFT) has been proposed for measuring velocity fields using several coils and a set of electrodes attached to the surface of the pipe. In this work, a velocity field reconstruction method for EMFT is proposed. The method uses a previously developed finite-element-based computational forward model for computing boundary voltages and a Bayesian framework for inverse problems. In the approach, the *v*_{z}-component of the velocity field along the longitudinal axis of the pipe is estimated on the pipe cross section. Different asymmetric velocity fields encountered near pipe elbows, solids-in-water flows in inclined pipes and in stratified or multiphase flows are tested. The results suggest that the proposed reconstruction method could be used to estimate velocity fields in complicated pipe flows in which the conventional EMFMs have limited accuracy.

This article is part of the themed issue ‘Supersensing through industrial process tomography’.

## 1. Introduction

Electromagnetic flow meters (EMFMs) are widely used in industry for measuring the volumetric flow rates of conductive continuous phase liquids and slurries. The conventional EMFMs consist of a circular pipe conducting the flow, a system producing a uniform magnetic field and a pair of point electrodes located on the opposite sides of the circular pipe diameter perpendicular to the magnetic field. Although the conventional EMFM has proven accurate for axisymmetric single phase flows, the lack of axial symmetry can lead to serious systematic errors in flow rate measurements. Moreover, the spatial changes in velocity field can induce further errors. Asymmetric axial flows are often encountered near valves and in pipe elbows and T-junctions, for example, and in multiphase flows.

In order to tackle the flow rate measurement of highly asymmetric single-phase flows, different multi-electrode approaches to electromagnetic flow measurement have been described [1,2]. In the multi-electrode approach, the voltages caused by the moving fluid are measured using multiple electrode pairs typically from eight to 16 electrodes in total. This method has been shown to handle the velocity field changes relatively well. However, the approach only measures the mean flow velocity and does not provide information on the local axial velocity field. This is a problem, especially with multiphase flows in which the actual volumetric flow rate of each phase is obtained by integrating the product of a local velocity field and a local phase volume fraction over the flow cross section. The velocity field information can be extremely important for characterizing multiphase flows such as oil-in-water flows in the oil industry.

For the reconstruction of velocity fields, different approaches have been developed. A modified filtered backprojection method was used in [3]. In [4], a method for reconstructing a one-dimensional stratified flow profile in seven vertical pixels by using precalculated weighting values was proposed. The method uses uniform magnetic fields produced by one pair of Helmholtz coils and seven pairs of electrodes to measure the induced boundary voltages. The method was compared against the flow results given by a Pitot tube with an error less than 8% [5]. An analytical method based on the discrete Fourier transformation assuming that the asymmetric velocity profile is polynomial was proposed in [6]. The analytical method uses both uniform and non-uniform magnetic fields. In [7], the method was extended to reconstruct velocity profiles containing both asymmetric and axisymmetric components given by polynomial and power-law representations, respectively.

Recently, a finite-element (FE)-based computational forward model was developed for simulating boundary potentials in electromagnetic flow tomography (EMFT) [8]. In EMFT, multiple electrode pairs are used to measure induced boundary voltages using at least two magnetic field excitations produced by pairs of coils. The developed forward model was compared against an independent model implemented using Comsol Multiphysics and several different parameters, such as mesh size, velocity profile and uniformity of the magnetic field, that affect the measured boundary voltages were studied. In [9], a reconstruction method for EMFT was proposed by using an inverse problem solution technique and preliminary results were shown.

In this work, a reconstruction method for estimating two-dimensional (2D) velocity field by using the developed forward model and a Bayesian framework for inverse problems is proposed. In the reconstruction method, only the *v*_{z}-component of the velocity field along the longitudinal axis of the pipe is estimated and the other components are assumed to be zero, that is *v*_{x}=*v*_{y}=0. In the reconstructions, 16 electrodes and two pairs of coils producing spatially non-uniform magnetic *B*-field inside the pipe are used resulting in a total of 32 voltage measurements. The proposed reconstruction method is tested using numerical simulations. Different velocity fields often encountered near pipe elbows, solids-in-water flows in inclined pipes, and in stratified or multiphase flows are tested.

The rest of the paper is organized as follows. In §2, the model for the boundary potentials is derived and the computational approximation of the model using the finite-element method (FEM) is described. The method for reconstructing the velocity field is derived in §3. Numerical simulation results testing the proposed method are shown in §4. Finally, the conclusion is given in §5.

## 2. Computing boundary potentials

In this section, the theoretical background for the computational model describing the potential distribution is described. First, the model is derived from the Maxwell's equations and then variational formulation of the model is presented. Finally, the variational formulation is discretized using the FEM. To take into account the external non-uniform magnetic fields an analytical solution describing the *B*-fields generated by circular coils are reviewed in the last subsection.

### (a) Forward model for potential distribution

Maxwell's equations for the time harmonic electric field intensity ** E**=

**(**

*E***) and for the magnetic field**

*r***=**

*H***(**

*H***) can be written as 2.1**

*r**a*and 2.1

*b*where is the position vector,

*Ω*is the spatial domain of interest,

**=**

*J***(**

*J***) is the current density,**

*r**μ*=

*μ*(

**) is the magnetic permeability,**

*r**ϵ*=

*ϵ*(

**) is the permittivity and**

*r**ω*is the angular frequency. The current density can be written as 2.2where

*σ*=

*σ*(

**) is the conductivity,**

*r*

*J*_{a}is the applied current,

**=**

*v***(**

*v***) is the velocity field and**

*r***=**

*B***(**

*B***) is the external magnetic field. The term**

*r**σ*

**is the so-called ohmic current and**

*E**σ*(

**×**

*v***) is due to the Lorentz force caused by the external**

*B**B*-field. By substituting equation (2.2) into equation (2.1b) and assuming that

*ω*≈0 and

*J*_{a}≈0 , i.e. there are no internal current sources, we can write equations (2.1) as 2.3

*a*and 2.3

*b*Further, since ∇×

**=0, the electric field is given by**

*E***=−∇**

*E**u*, where

*u*=

*u*(

**) is the scalar electric potential. Substituting**

*r**σ*

**=−**

*E**σ*∇

*u*into equation (2.3b) and taking the divergence, we obtain 2.4since the divergence of the curl operator is always zero. When the conductivity

*σ*is known or constant, and

**and**

*v***are known, the potential distribution**

*B**u*inside the pipe segment

*Ω*can be solved from equation (2.4) with appropriate boundary conditions. This is the form that is used in traditional electromagnetic flow measurements [10]. The boundary conditions are chosen as 2.5

*a*2.5

*b*and 2.5

*c*where is the outward unit normal, ∂

*Ω*=∂

*Ω*

_{ends}∪∂

*Ω*

_{w}is the surface of the domain

*Ω*, ∂

*Ω*

_{ends}are the ends of the pipe and ∂

*Ω*

_{w}is the surface (wall) of the pipe. The boundary conditions describe (i) no current flows through the surface of the domain, (ii) zero velocity at the pipe wall (no-slip condition), and (iii) vanishing magnetic field at the ends of the pipe which can be justified when the computational domain of the pipe is long enough. Equation (2.4) can be solved analytically in simple geometries with an axisymmetric velocity field perpendicular to the constant magnetic field. More complicated geometries and velocity fields typically require using numerical methods. In the following, the FEM is used for a numerical approximation of equation (2.4).

### (b) Finite-element approximation

Let us now derive the FE approximation of the equation (2.4) with the boundary conditions (2.5). By multiplying both sides with a test function *ϕ*=*ϕ*(** r**), integrating over the domain

*Ω*and using Green's formula we can write [8,9] 2.6where d

*S*is the infinitesimal surface element of the surface ∂

*Ω*. Using the zero boundary conditions equations (2.5), the variational formulation is 2.7The FE approximation of the variational formation equation (2.7) is obtained by approximating the solution as a linear combination of the basis function , choosing the test functions

*ϕ*(

**)≈**

*r**ϕ*

_{h}(

**)=**

*r**φ*

_{i}(

**),**

*r**i*=1,…,

*N*and approximating the velocity field as 2.8where

*N*is the number of nodes in the mesh. The FE approximation can be written in a matrix form as 2.9where and are the vectors of the potential coefficients and the velocity coefficients at the nodes of the FEM mesh, respectively. The matrix is 2.10and the components of the block matrix are 2.11 2.12and 2.13where

*i*,

*j*,

*k*=1,2,…,

*N*. After fixing one of the potentials, say

*u*

_{p}, the matrix equation for the potential coefficients can be solved.

Because *A* is symmetric positive definite matrix, the Cholesky factorization can be computed and the boundary potentials *U* can be solved as
2.14where is the measurement operator interpolating and averaging the potential *u* on the electrodes, *N*_{M}=(*N*_{e}×*N*_{p}) is the number of measurements, where *N*_{e} is the number of electrodes and *N*_{p} is the number of magnetic field excitations.

The goal of EMFT is to estimate the 2D velocity field on the measurement plane. Therefore, only the *z*-components of the velocity field coefficients are estimated and the unknown is written in a form *v*:=[*v*_{z,1},*v*_{z,2},…,*v*_{z,N2D}], where *N*_{2D} is the number of nodes in the 2D measurement plane. The corresponding measurement model can be written as
2.15where is the projection operator extending the *z*-component of the 2D velocity field along the *z*-axis in *Ω*. Note that the matrix *H* needs to be computed only once for the same measurement geometry and magnetic fields. Finally, the linear additive measurement model can be written as
2.16where *e* is the additive measurement noise.

### (c) Computation of the B-field of a coil

To compute the non-uniform magnetic *B*-field in the pipe generated by external coils, we use analytic expressions for the *B*-fields in Cartesian coordinates. The expressions are derived for a single circular current loop but a coil can be modelled by approximating the total magnetic field of the coil as the sum of the *B*-fields of individual loops [8]. The *B*-field components of a single current loop located in the *xy*-plane and centred at the origin with a radius *a* carrying a current *I* can be calculated as follows [11]:
2.17*a*
2.17*b*
2.17*c*where *C*=*μ*_{0}*I*/*π*, *μ*_{0}=4*π*⋅10^{7} Tm A^{−1} is the permeability of free space and *E*(*k*) and *K*(*k*) are complete elliptic integrals with modulus *k* defined by
2.18with the following substitutions *α*^{2}=*a*^{2}+*r*^{2}−2*aρ*, *β*^{2}=*a*^{2}+*r*^{2}+2*aρ*, *r*^{2}=*x*^{2}+*y*^{2}+*z*^{2} and *ρ*^{2}=*x*^{2}+*y*^{2}. Other geometric configurations of the coil can be computed through geometric transformations.

## 3. Reconstructing velocity distribution

The reconstruction of a velocity field in EMFT can be approached using an inverse problem solution technique. More precisely, a statistical (Bayesian) framework for inverse problems is used. In the Bayesian framework for inverse problems, all unknowns are treated as random variables and degree of uncertainty is encoded into the models using probability distributions [12]. The solution of the Bayesian inverse problem is the posterior distribution *π*(*v*|*U*_{meas}) which gives the probability of the unknown *v* for a given measurement realization *U*_{meas}. The posterior distribution is constructed using statistical models for the unknown *v* and the measurement process given by the prior *π*(*v*) and the likelihood *π*(*U*_{meas}|*v*), respectively. Once the posterior distribution is constructed different point and uncertainty estimates, such as conditional mean, maximum a posterior estimate or credibility intervals, can be computed. For more details about the Bayesian framework for inverse problems, see for instance [12,13].

### (a) Posterior distribution

The posterior distribution is
3.1In the case of the additive measurement model, equation (2.16), and assuming that *v* and *e* are mutually independent, the posterior can be written as
3.2where *π*_{e} is the probability distribution of the noise and *U*_{meas} are the measured boundary voltages on the electrodes. If the noise *e* and the unknowns *v* are modelled as Gaussian random variables, and , the posterior is
3.3where and are the Cholesky factors of the inverse covariances of the noise and the parameters, respectively. The maximum a posterior (MAP) estimate can be computed from the posterior as
3.4This problem is equivalent to the minimization problem
3.5
3.6which is a generalized regularized linear least-squares problem. It can formally be solved as [12]
3.7Therefore, online estimation of *v* requires only a few matrix–vector multiplications once the approximate inverse of has been precomputed before any measurements take place.

### (b) Prior model

Prior model for the velocity field is constructed by assuming that it consists of two mutually independent parts [14]
3.8where *v*_{in}(** r**) is spatially varying with zero mean modelling the spatial changes of the velocity field and

*v*

_{mean}is spatially constant with non-zero mean modelling the mean of the velocity field. Since velocity fields are typically spatially smooth this information can be encoded into the covariance matrix of the spatially varying part by setting [15–17] 3.9

where *σ*^{2}_{in,v} is the variance of the spatially varying part and *l* describes the correlation length affecting the smoothness of the velocity field. The spatially constant part is modelled as *v*_{mean}=*q***1**, , where **1** is a vector of ones. Since *v*_{in}(** r**) and

*v*

_{mean}are modelled as mutually independent the prior distribution can be written as where . The parameters

*v*

_{*}and can be chosen based on

*a priori*knowledge about the possible range of the mean velocity of the flow, and

*σ*

^{2}

_{in,v}and

*l*can be chosen based on the expected range of changes of the velocity field around the mean and smoothness of it.

## 4. Results

The proposed velocity field reconstruction method was tested using numerical simulations. First, the simulation geometry, magnetic fields and parameters of the noise and prior models are described. Then, different flow test cases are considered.

### (a) Simulation set-up and parameters

In the simulations, a pipe segment shown in figure 1 was considered. The diameter and length of the pipe were 5 cm and 15 cm, respectively. Four coils with radius of 4 cm and 16 electrodes were located at the middle of the pipe on its surface. The FE mesh for the computation of the potential distribution contained approximately 60 000 nodal points and 333 000 linear tetrahedral elements. The simulated measurement data were generated using a very dense mesh with approximately 270 000 nodal points and 1 400 000 elements. The convergence of the FE approximation as a function of number of nodal points has previously been verified [8]. Therefore, the discretization error can assumed to be small. For the estimation of the velocity field a 2D mesh with *N*_{2D}=581 nodal points and 1096 linear triangular elements was used. Different meshes were used for computing the simulated data and the potential distribution when reconstructing the velocity field to avoid the so-called ‘inverse crime’ [12]. This term refers to obtaining unrealistically good results when using the same meshes since all the possible discretization and modelling errors are neglected in the reconstruction process.

A single-phase water flow was considered with the conductivity *σ*=0.2 mS cm^{−1}. The coils were used to produce two orthogonal magnetic fields (*N*_{p}=2) sequentially by exciting two opposite coils at a time and the induced potential differences between adjacent electrodes (*N*_{e}=16) were measured resulting to *N*_{M}=32 measurements in total. The current amplitudes for the coils were [2,0,2,0] A and [0,−2,0,−2] A corresponding to the first and second excitation, respectively. The coils were modelled as 35 horizontal layers of 23 vertical windings of copper wire having a thickness of 0.68 mm. The outer radius and width of the coil were 5.56 cm and 2.4 cm, respectively. The resulting *B*-field was computed as described in §2c. The *B*-field in the pipe produced by one pair of coils in shown in figure 2. Figure 3 shows the *B*-field on the measurement plane defined by the electrode ring. As it can be seen the *B*-field is non-uniform having the maximum amplitude of 16 mT near the active coils. Moreover, the magnetic field is very weak at the ends of the pipe indicating that the zero boundary condition equation (2.5c) is a feasible choice.

Once the potential differences were simulated, independent and identically distributed Gaussian zero mean random noise (*e*_{*}=0, ) with standard deviation () was added to the data giving a noisy realization of the simulated measurements. In the reconstructions, larger standard deviation () was assumed. The simulated noisy measurements were used to reconstruct the *z*-component of the 2D velocity field at the pipe cross section.

The mean and standard deviation of the prior model were chosen as *v*_{*}=2 m s^{−1} and *σ*_{mean,v}=0.5 m s^{−1}. These parameters indicate that the mean of the flow is within an interval of *v*_{*}±3*σ*_{mean,v} m s^{−1}=[0.5,3.5] m s^{−1} with prior probability of 99.7%. Further, the standard deviation and the correlation length of the spatially varying part of the prior model were chosen as *σ*_{in,v}=1 m s^{−1} and *l*=5 cm, which is equal to the diameter of the pipe. One row of the covariance matrix *Γ*_{v} as a function of distance between the points is shown in figure 4 using three different correlation lengths. Using too small a correlation length compared with the diameter of the pipe will allow rapid spatial changes in the velocity field. On the other hand, using an excessively large correlation length enforces too much smoothness, leading ultimately to a flat velocity field. Both of these extreme cases are unphysical and, therefore, the correlation length should be roughly chosen as the diameter of the pipe. The chosen parameters correspond to a smooth velocity field with maximum and minimum values of ±3*σ*_{in,v} around the mean flow. In all of the simulation cases, the prior model was the same.

### (b) Velocity field simulations

Five different flow test cases given in table 1 were considered. In Case 1, the true velocity field was parabolic and it was computed as
4.1where *r* in the distance from the centre of the 2D cross section of the pipe and *R* is the radius of the pipe. In Case 2, axisymmetric annular flow was considered
4.2where *v*_{jump}=0.5 m s^{−1} and *R*_{1}=0.015 m. In Case 3, asymmetric velocity field with both positive and negative velocities was considered. The true velocity field was generated as
4.3where
4.4and *r*_{0,1}=(0,0.02) m, *r*_{0,2}=(0,−0.02) m, *σ*_{1}=1 m and *σ*_{2}=0.5 m. In Case 4, the asymmetric stratified velocity field was considered
4.5In Case 5, the asymmetric velocity field given by a fourth-order (quartic) polynomial [6] was considered
4.6where *v*_{f} is a multiplicative factor defining the maximum velocity and
4.7
4.8
4.9
4.10The 2D velocity fields were then extended to 3D using the projection operator . In Case 6, the velocity field after a pipe elbow shown in figure 5 was computed using a computational fluid dynamics (CFD) simulation. The CFD simulation was computed using the Reynolds-averaged Navier–Stokes equations with *k*−*ω* turbulence model and OpenFOAM CFD software [18]. The computed steady-state velocity field (fully developed flow) is shown figure 6. The velocity field at the 2D measurement plane was used as a true velocity field. In all of the test cases, *v*_{max} was chosen such that the mean velocity was approximately 1 m s^{−1}.

### (c) Velocity field reconstructions

The 3D potential distributions inside the pipe were computed using equation (2.9) with the known velocity and magnetic fields for different cases. The potential distribution on the measurement plane for Case 5 is shown in figure 7. As it can be seen the potential distribution is non-symmetric due to asymmetric velocity field. Further, the noisy voltage measurements between adjacent electrodes were computed from the potential distribution using equations (2.15) and (2.16). These noisy measurements for different cases are shown in figure 8. The maximum value of the measured voltage is from 0.1 to 0.2 mV depending on the velocity field for 5 cm pipe diameter using coils of 4 cm radius and magnetic field having 16 mT maximum amplitude.

The velocity reconstructions were computed using equation (3.7) and the results are shown in the right column of figure 9. The averages of the velocity estimates were compared against the mean velocities given by a simulated conventional EMFM. Therefore, a voltage measurement in the conventional EMFM was simulated using a uniform magnetic field and one pair of electrodes at the opposite sides of the pipe orthogonal to the magnetic field. The mean velocity in this case was computed as [10]
4.11where *U* is the measured voltage difference, *D* is the diameter of the pipe separating the electrodes and *B* is the magnitude of the uniform magnetic field. Equation (4.11) can be derived from equation (2.4) with certain assumptions including that the velocity field is axisymmetric. The mean velocities and their relative errors against the true mean are given in table 2.

The velocity reconstructions in figure 9 show good agreement with the true velocity fields. In case 2, the velocity field estimate is blurred since all the axisymmetric velocity fields produce almost the same boundary data using a nearly uniform magnetic field. To better distinguish between different axisymmetric velocity fields using more non-uniform magnetic field could be beneficial. The mean velocities computed from the estimates are very close to the true mean having smaller error compared with the mean velocity given by the simulated EMFM due to asymmetric velocity field.

In all of the test cases, the true mean velocity is approximately 1 m s^{−1}. However, the velocity field in different test cases is quite different, indicating that measuring only the mean flow is not enough to characterize highly asymmetric pipe flows. Therefore, this suggests that the proposed tomographic reconstruction approach could be used to estimate 2D velocity fields in complicated pipe flows.

## 5. Conclusion

In this work, a method for reconstructing a 2D velocity field in EMFT was proposed. In EMFT, multiple electrode pairs are used to measure induced boundary voltages using several magnetic coils. The method uses FE-based computational forward model to compute boundary voltages and a Bayesian approach to inverse problem to reconstruct the velocity field. The proposed linear reconstruction method requires only a few matrix–vector multiplications of moderate size once the matrices corresponding the FEM forward model, prior and noise models have been precomputed before any measurements take place. Therefore, the method allows for online estimation of the velocity field. The method was tested using numerical simulations in various test cases with complicated asymmetric velocity fields. The results indicate that the reconstruction method could be used to estimate velocity fields of complicated pipe flows. Further studies will be carried out in order to test the proposed approach with real measurement data with various types of velocity fields.

## Data accessibility

A Matlab package containing the data and scripts to reproduce the figures is available in the Dryad Digital Repository at doi:10.5061/dryad.q2102.

## Authors' contributions

All the authors contributed to developing the velocity reconstruction method. The idea was developed by M.V. Numerical implementation and computation of the results were done by O.L. and K.K. O.L. wrote the manuscript with comments from K.K. and M.V. All the authors have read and approved the manuscript.

## Competing interests

The authors declare that they have no competing interests.

## Funding

This work has been supported by the European Regional Development Fund (ERDF project no. A70738) and by the Academy of Finland (Finnish Centre of Excellence in Inverse Problems Research, project no. 250215).

## Acknowledgements

The authors would like to thank the anonymous referees for useful comments and suggestions.

## Footnotes

One contribution of 10 to a theme issue ‘Super-sensing through industrial process tomography’.

- Accepted February 25, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.