Improving sand flow rate measurement using the wavelet transform and ultrasonic sensors

In numerous gas wells around the world, sand is produced in addition to the gas. Production of sand creates several issues with the production-well/upstream/downstream facilities. The main issue with sand production is the erosion of various transport equipment such as chokes and pipe. Severe erosion might lead to damage to pressure holding devices. The same might lead to leakage of hydrocarbon products which is an environmental issue. In the worse cases, it might lead to rupture or explosion of such equipment. In addition, sand may stagnate in various vessels (e.g. separators) and cause a reduction in gas processing capacity (Oyeneyin, 2015). This causes increased operational costs for removal of the sand and/or repair of the wellhead facility (Jaimes Plata et al., 2012). For instance, sand production may cause erosion in a subsea choke which may necessitate costly replacement of the choke module in the subsea wellhead.

Methods to prevent the consequences of sand production can be categorized as two main techniques. The first technique is to prevent sand production using various downhole prevention devices. For instance sand screens, gravel packs, slotted liners, frack-and-pack and consolidation methods can be used to prevent sand production (Matanovic et al., 2012; Kaiser et al., 2000; Ben Mahmud et al., 2020). Failure of such downhole devices (i.e. due to sand erosion) can result in a sudden increase in sand flow rate. Therefore a sand measuring device for early detection of sand production is very useful to alert the operator to reduce/stop gas production from that well.

The second technique is sand management in which a small amount of sand production is allowed. In this case, suitable precautions need to be considered to limit the side effects of sand production. For instance, the amount of sand production would need to be less than a certain allowable limit. Also some consideration of sand production needs to be taken into account during the design of the facility. For example, some form of sand removal facility needs to be considered in the incoming separators of the downstream facility. An erosion allowance may be considered to be included in the wall thickness of the pipe. This sand control management is vital for accurate sand measurement to ensure that the sand production is less than the allowable limit for each specific field.

To correctly apply sand prevention and sand management techniques, accurate sand flow rate measurement is necessary. It is reported that the major concern of the operators about sand flow rate measuring devices available in the market is their poor measurement accuracy (Emiliani et al., 2011). This research is conducted to improve sand flow rate measurement accuracy. Such a method can be used for early detection of sand production. Early detection can notify the operator of any problem in sand prevention techniques (e.g. damage to downhole sand prevention devices). Also, improving sand measurement leads to enhancement of the accuracy of sand flow measurement. The same will greatly support operators to apply sand management strategies as in such a method, accurate sand flow measurement is needed. The other problem with existing sand detectors is the difficult and lengthy procedure for calibrating these devices. Therefore, finding a method with simpler calibration procedure is needed by the industry. This is another aspect which is addressed by this paper.

Finally, one difficulty with commonly used sand detectors is the requirement for having fluid velocity as a pre-requisite for sand flow rate estimation. For instance, in several applications such as subsea wellheads, there is no flow meter by which the fluid velocity can be measured. Therefore, having a method for sand flow rate measurement in which the fluid velocity is not needed can really simplify sand measurement in these applications and the paper also presents a method in which fluid velocity measurement is not needed.

This paper focuses on sand flow rate measurement using ultrasonic sensors. Advanced signal analysis of ultrasonic signals using the wavelet transform is then used for improved sand flow rate measurement. Figure 1 shows the research flow chart for the proposed method of determining sand flow rate.

Figure 1:

A background review is now provided in which the ultrasonic sensors and wavelet transform method used in this research are explained. Then, the details of the test setup used for this research is elaborated. Then, the results and discussions are presented about applying a wavelet transform to ultrasonic sensor data for measuring sand flow rate.

Ultrasonic sensors are widely used by the petroleum production industry for various measurement purposes such as corrosion (Feydo et al., 2017), sand rate measurement (Musa et al., 2005), and non-destructive testing (Sinclair and Malkin, 2020). While Acoustic Sand Detectors (ASDs) in the industry use ultrasonic sensors there are a number of concerns about the accuracy of these measurements (Emiliani et al., 2011). This paper explains a new approach to improve the accuracy by applying advanced signal analysis methods to the ultrasonic sensor output data.

Typically ASDs are installed after a bend in the pipe through which the sand is flowing (Gao et al., 2015). When a gas containing some sand particles passes around a bend in the pipe, the sand particles make impact with the pipe wall due to their inertial energy. This causes the creation of ultrasonic waves that are then picked-up and translated by the ASD into electrical signals. One of the commonly used ASDs in the market uses the following relationship to calculate sand flow rate (Ibrahim and Haugsdal, 2008): (1)(1)Q=M−ZS.

In this relationship, Q is sand flow rate, M is the output of the ASD, Z is a “Zero” value representing the offset due to background noise, and S is a “Step” value showing the gain value of an ASD output to the sand flow rate. The Step and Zero parameters are obtained during the calibration process for each application. Calibration is performed using the data collected by various experiments at various fluid velocities and sand flow rates. As a result, a lookup table is obtained for each of the zero and step parameters. These lookup tables show the recommended values of step and zero parameters at various velocities. Such calibrations are a lengthy and difficult activity to practically perform.

Other research has suggested other relationship to calculate the sand flow rate using the raw value obtained from an ASD. For instance, Gao et al. (2015) proposed using the following relationship: (2)(2)Q=M−ZKV2.

In this relationship, Q is sand flow rate, M is the output of the ASD, Z is the Zero value to compensate for background noise, V is the fluid velocity and K is a constant value.

Sampson et al. (2002) suggested the following relationship for calculating sand flow rate using an ASD: (3)(3)Q=M2−Z2CV2.

In this relationship, Q is sand flow rate, M is the output of the ASD, Z is the Zero value to compensate for background noise, V is the fluid velocity and C is a constant value.

Ultrasonic signals are a subgroup of acoustic signals with frequencies above human hearing capabilities (i.e. above 20 kHz). There are two common types of propagating ultrasonic signals which have longitudinal and transverse propagation. In longitudinal propagation, particles are oscillating along the same direction as the wave propagates. In transverse (shear) propagation, the direction of particle oscillation is perpendicular to the direction of wave propagation. In gas and liquids, the ultrasonic waves propagate using longitudinal waves only, while ultrasonic waves propagate in solids using both longitudinal and transverse forms (Boyd and Varley, 2001).

After an ultrasonic sensor is installed on a pipe, a signal processing technique needs to be adopted in order to correlate the output of the ultrasonic sensor with the sand flow rate. Various researchers have used different signal processing techniques to measure sand flow rate. For instance, Gang et al. (2015) used a Short-term Fourier Transform (STFT) to analyse the signal from acoustic sensors. STFT is an extension of the Fourier Transform for non-stationary signals (those signals which change characteristics over time). STFT converts a signal from the time domain to the time-frequency domain.

The Wavelet Transform is a superior method for analysing non-stationary signals and is used in various applications such as image processing, data compression, de-noising, etc. (Shukla, 2013). Similar to the STFT, this method uses the time-frequency domain to analyse non-stationary signals. Since sand flow rate is a non-stationary signal and due to the common use of the wavelet transform methodology, this method is used for signal analysis in this paper.

The Wavelet transform has two main categories i.e. the Continuous Wavelet Transform (CWT) and the Discrete Wavelet Transform (DWT). The CWT of a continuous signal is defined by Polikar (2006): (4)(4)CWTxψ(τ,s)=1|s|∫x(t)ψ*(t−τs)dt,(4)where ψ(t) is the wavelet function, CWT is a function of τ and s parameters and τ is a translation parameter while s is a scale parameter. Translation represents the time information while scale represents the inverse of frequency.

The Wavelet function (ψ(t)) is a rapidly decaying wave shape oscillatory signal with a zero mean. There are various types of wavelet function. Figure 2 represents one of the commonly used wavelets called a Morlet wavelet (Cohen, 2019) while another commonly used is the Ricker wavelet.

Figure 2:

As shown in this figure, the wave shape is non-zero for a limited period and is symmetrical about its centre. In CWT, a signal is represented by a series of wavelet functions. After using a wavelet transform, a signal in the time domain is converted to a signal in the time-frequency domain. Since there are indefinite time-frequency range selections to transform τ and s (equivalent to any arbitrary selection in time-frequency space), calculating CWT needs a lot of computation. In on-line industrial applications such as sand flow rate measurement, it is preferred to have minimal calculations when the calculation is to be performed by the measuring device. For this reason, researchers have tried to reduce the calculations by limiting the selections in the time-domain space. For instance, Burt and Adelson (1983) has proposed Multiresolution analysis which led to the creation of the Discrete Wavelet Transform (DWT). The DWT is used for discrete signals by using a series of high frequency and low frequency band-pass digital filters. Convolution is often used during the process of passing a discrete signal through a digital filter. The convolved relationship for passing a signal through a digital low pass filter used in the DWT is given by Polikar (2006): (5)(5)x[n]*h[n]=∑−∞∞x[k].h[n−k],(5)where h[n-k] is the impulse response of the digital low pass filter (shifted by k samples in time). x[n] is the discrete signal and in the DWT, the discrete signal passes through a digital high pass filter.

The output signals from low pass and high pass filters are subsampled by discarding every other sample. The subsampled output of the high pass filter constitutes the first series of DWT coefficients. Then, the low pass and high pass filters are again applied to the subsampled signals from the low pass filter as follows (Polikar, 2006): (6)yhigh[n]=∑−∞∞z[n].g[2n−k]ylow[n]=∑−∞∞z[n].h[2n−k],(6)(6)where h[n] represents the impulse response of the low pass filter and g[n] represents the impulse response of the high pass filter while z[n] represents the subsampled output of the low pass filter from the previous stage of the DWT. The subsampled output of the high pass filter forms the second series of DWT coefficients. This procedure is continued until the required levels of the DWT have been reached. Figure 3 demonstrates the procedure to calculate DWT coefficients (Polikar, 2006).

Figure 3:

Here we briefly explain the test set-up used for this experiment. The test facility consisted of an open flow loop in which a mixture of air and solid particles were passed through a pipe. First the air from the surrounding environment was pumped through the pipe using a compressor. Using a variable speed drive (VSD), the velocity of the air was able to be adjusted. Using a special sand hopper, sand particles could be injected into the pipe and the rate of sand injection could be remotely adjusted using a suitable computer interface. At the end of the pipe, the air/sand mixture entered into a collection basket where the sand particles were collected and the air filtered in order to pass to the surrounding environment. A temperature/velocity sensor was installed in the pipe in order to measure the temperature and velocity of the air passing along the pipe. The output of this temperature/velocity sensor was sent to a National Instruments USB-8009 data acquisition (DAQ) system, which transferred the velocity and temperature to the main data recording computer. This DAQ performed sampling at a rate of 48 kS/s (thousand samples per second) and had 14-bit resolution recording. Considering the relatively slow variation of temperature/velocity, the speed and resolution of this DAQ was adequate for this application.

An ultrasonic sensor was installed after a bend in the pipe. This ultrasonic transducer was mounted on an aluminium support which was curved to fit snugly with the pipe curvature on the pipe side and flat on the other so that good contact was made with the ultrasonic sensor on the other side. Figure 4 shows the test facility used for this experiment. For comparison, an Acoustic Sand Detector was installed on this pipe and was used for other experiments which are not covered in this paper (Seraj and Evans, 2020).

The ultrasonic contact-type transducer from Olympus Company could measure ultrasonic signal at a centre frequency of 1 MHz and had a diameter of 0.5 inch. In order to have accurate measurements, the sensor requires good contact with the pipe under test, which was the reason for using the aluminium support as explained earlier. Figure 5 represents the block diagram of the flow loop test facility (Seraj and Evans, 2020).

Figure 4:

The output of the ultrasonic sensor was connected to a Femto amplifier which amplifies the signals with an adjustable gain from 10² to 10⁸. The amplified signal is then passed to the data acquisition system (DAQ) which was a National Instruments Model PXI-1033. This DAQ samples the ultrasonic signal at a sample rate of up to 5 MS/s with a resolution of 20-bits. Having such a high sampling rate with good resolution assisted the sampling of rapidly changing ultrasonic signals.

Labview software was installed on the main computer which stored the sampled signal in a structured data format called Hierarchical Waveform Storage (HWS). This information was later read in MATLAB. Then using the signal processing module in MATLAB and other scripting capabilities in this software, DWT was applied to the recorded ultrasonic sensor signals and the results were further analysed to find a relationship between DWT coefficients and the sand flow rate. The results of these analyses are explained in the following sections of this paper. Various components of the flow loop test are shown in Fig. 6.

Figure 5:

Figure 6:

In this research, experiments were conducted with the flow loop running at various air velocities and sand flow rates. The air velocities examined were up to 27 m/sec and the sand flow rate was in the range of 0 to 5 g/sec- typical rates for sand production issues.

The wavelet transform method was used to analyse the data obtained from an ultrasonic sensor. The Wavelet Transform converts a signal from the time domain to the time-frequency domain. Converting the signal to the time-frequency domain helps to identify the frequency components of the signal during various time periods. To better illustrate this issue, Fig. 7 shows the Spectral Analysis output of the wavelet transform of ultrasonic signal collected from the transducer at a fluid velocity of 7 m/sec and sand rate of 35 g/sec. The figure shows the magnitude of signal using different colours where in this case, the magnitude is in units of V2 (Volts to the power of 2). Higher magnitude is shown with a yellowish colour in this figure.

Figure 7:

The dotted cone shape in this figure is the “cone of influence”. The information inside the cone (above the dotted line) is reliable. The information near the cone (near dotted line) and outside the cone (below the dotted line) is not reliable due to “boundary effects” (Addison, 2002). This is due to the fact that the signal is finite (with a limited amount of sampling).

As can be seen from this figure, the signal shows higher magnitude in frequencies around 16 kHz. Therefore this study focuses on evaluating the magnitude of the wavelet transform in a frequency range around 16 kHz and tries to correlate this strength with sand flow rate. Also, it can be seen that the strength of the signal at a frequency of 16 kHz changes at various times. For instance, the magnitude of signal at 16 kHz is more in the time domain around 0.5 sec compared with that at 1.5 sec. Note that in this experiment, the ultrasonic signal is recorded for 2 seconds. This shows that the characteristics of the signal (e.g. frequency content) varies over time. Therefore wavelet transform method is a suitable tool for analysing the ultrasonic non-stationary signals which vary over time.

As mentioned before when using the DWT, the frequency domain can be divided into various intervals. In this application where ultrasonic sensor with bandwidth of 1 MHz is used, the frequency intervals were defined in Table 1.

As can be seen in this table, by moving from one frequency interval to the next interval, the frequency range was divided by two. The strength of the DWT versus sand flow rate in these frequency ranges are shown in Fig. 8.

Figure 8:

The X-axis in Fig. 8 shows the sand flow rate in grams per second. The Y-axis in this figure shows the sum of the squares of DWT coefficients at that frequency interval (“strength of wavelet coefficients” in this paper). It should be noted that the ultrasonic signal at various velocities are averaged for each sand flow rate before applying the DWT method. The effect of velocity is discussed later in this paper.

The strength of the wavelet coefficients is stronger in frequency interval 5 (frequency range from 31.25 to 62.5 kHz) and 6 (frequency range from 15.625 to 31.25 kHz).

In order to combine the response in these two frequency ranges, the magnitude of wavelet coefficients in these frequency ranges were summed. Figure 9 shows the sum of wavelet coefficients’ strength across these two frequency ranges versus the sand flow rate.

Figure 9:

In this figure, the asterisk represents the strength of DWT coefficients versus sand flow rates. The dashed line shows the interpolation between these values. Note that the strength of the wavelet coefficient in the frequency range 15.625–62.5 kHz has a relatively linear relationship with the sand flow rate. This shows that there is a correlation between the wavelet coefficient strength in these frequency ranges and the sand flow rate. Considering this correlation has a nearly linear relationship, a best-fit line shows the strength of wavelet coefficients in the frequency range 15.625–62.5 kHz is proportional to the sand flow rate. Figure 10 shows the best-fit line to these data.

Figure 10:

The best-fit line passes through the origin, which is logical since at no sand flow rate there are no sand vibrations and hence no transducer output (it is considered that outliers are within error limits). Consequently the linear relationship between the strength of the wavelet coefficient across the frequency range 15.625–62.5 kHz is shown below: (7)Q=aPa=2837.(7)

In this relationship, Q is the sand flow rate (in g/s), and P is the sum of squares of wavelet coefficients of ultrasonic signal in the frequency range 15.625–62.5 kHz (in Volt²), while a is the proportional coefficient (i.e. the slope).

This relationship provides a very easy way to estimate the sand flow rate using a very simple (and cheap) ultrasonic sensor with the discrete wavelet transform. In industry, it is preferred to have a simple relationship so that there is a simple implementation of this recording method in the control system.

Parameter a in this relationship must be obtained for each application site which can be regarded as a site calibration. Since there is only one parameter to be estimated during calibration rather than two, this is an easy task as the installation should ideally have very few calibration tests to keep costs at a minimum. In the DWT method presented in this paper, only one parameter must be obtained while in the commercial ASDs, two look-up tables are needed for step and zero values. Therefore, less experiments are needed during calibration of the sand measurement method proposed in this paper.

To examine how much the strength of the wavelet transform coefficients would be across the frequency range 15.625–62.5 kHz at different velocities, MATLAB analysis was used. Figure 11 shows the wavelet coefficient strength across the frequency range of 15.625–62.5 kHz at the different velocities.

Figure 11:

This figure confirms that although there is some variation from the linear relationship best-fit to the data, the trend of the data is along the best-fit line.

Using the linear relationship (5), the test results indicate a sand flow rate maximum deviation of 17.7%, which is considered acceptable for raw data measurement, and it is noted that at higher flow velocities, these values have lower error than at low velocity values. To compare the results with the commonly used ASDs, one of these ASDs was also installed on the flow loop test facility used in this experiment (refer to block diagram of the test setup in Fig. 5). The output of the ASD was collected for various sand flow rates and fluid velocities. Relationship (1) was used to estimate sand flow rates that has more than 25% uncertainty in calculation compared with 17.7% of this method. Therefore, the proposed method in this paper provides more accurate sand flow rate measurement compared with the commercial ASDs commonly used in the industry.

Thus, it is considered that this relationship will remain true for all higher velocity values (greater than or equal to 5 m/sec) tested in this paper. At these higher velocities, the sand particles properly mix with the gas phase and produce a homogeneous mixture. At lower velocities (less than 5 m/sec), a portion of the sand particles still mix in the gas phase, while other portions of sand particles are gathered at the bottom of the pipe and slowly move (ripple) along on the base of the pipe. Therefore at higher velocities, there is more chance that the sand particles hit the pipe wall at the bend where the ultrasonic sensor/ASD is installed. Therefore measuring sand flow rate using an ultrasonic sensor is more accurate at high velocity values (equal to or greater than 5 m/sec).

As shown in Equations (1)–(3), various relationships to calculate sand flow rate from ASD data need the fluid velocity. For instance, in Equation (1), the step and zero values are obtained from a look-up table which specifies these values at various fluid velocities. Therefore to use this relationship, fluid velocity is needed to estimate the zero and step value. Also in relationship (2) and (3), the fluid velocity is directly used in the sand flow rate calculation. Typically flow meters are used to measure fluid velocity. Where a flow meter is not available, there is concern how to measure fluid velocity needed for sand flow rate measurement. The advantage of using the method proposed in this paper – relationship (5) – is that it does not need fluid velocity as input for sand flow rate calculations.

The other concern operators have with commonly used ASDs is the lengthy and demanding procedure for calibration. For instance in relationship (1), the ASD must be calibrated at various velocities and sand flow rates to obtain the look-up tables for zero and step values. Typically it is not practically possible to change the fluid velocity in an operating gas field. For instance, when ASD is used around a wellhead, then to get the high velocity, it is necessary to open the choke valve. But during the early life of the field when the ASD is calibrated, the reservoir pressure is relatively high. So, a wide opening of the choke valve might result in relatively high pressure in the wellhead area which may not be acceptable for the downstream facility. Therefore, it is likely not possible to reach high velocities for ASD calibration, and so the calibration is incomplete.

In the proposed method in this paper, only one parameter (coefficient “a” in relationship (5)) must be calibrated. To obtain this parameter, there is no need to calibrate the sand measuring device over a wide range of fluid velocities and sand flow rates. Instead, the calibration is performed over a few sand flow rates. This greatly simplifies the operator’s calibration procedure.

This paper proposes to use an ultrasonic sensor and an advance signal processing method, the Discrete Wavelet Transform (DWT), to measure sand flow rate. This paper demonstrates that there is a linear relationship between sand flow rate and energy using the DWT across the frequency range 15.625–62.5 kHz. It is shown that this relationship can measure sand flow rate independent of the fluid velocity. It is also noted that this method provides more accurate results at higher velocities (equal to or greater than 5 m/sec).

This linear relationship which can be applied in real time, and provides a simple and inexpensive solution to accurately calculate sand flow rate using off-the-shelf equipment. Compared to the costly commercial products, this solution provides a more cost effective and accurate sand measurement method.

This paper addresses the concerns about the accuracy of commonly used ASDs in the market. The proposed DWT and ultrasonic sensor method can help various operators to measure the sand flow rate more accurately. Such an improvement in sand flow rate measurement leads to better implementation of sand prevention and control techniques.

Also, the proposed method for sand flow rate calculation removes some of the concerns about calibration of the existing ASDs on the market. The proposed method does still require site-specific calibration (which is an accepted part of commercial products), but considering there is a requirement to obtain only one parameter during calibration, the proposed method is much simpler that of the commercial products.

Furthermore, in the proposed method, there is no need to measure fluid velocity as a pre-requisite to estimate sand flow rate, while for most of the ASDs in the industry, fluid velocity is needed. Since in a number of applications such as at subsea wellheads, the flow meter may not be available to measure fluid velocity, using this proposed method will help the operators to accurately measure sand flow rate.