Note: Descriptions are shown in the official language in which they were submitted.
CA 02658810 2013-10-02
WET GAS MEASUREMENT
TECHNICAL FIELD
This description relates to flowmeters.
BACKGROUND
Flowmeters provide information about materials being transferred through a
conduit.
For example, mass flowmeters provide a direct indication of the mass of
material being
transferred through a conduit. Similarly, density flowmeters, or
densitometers, provide an
indication of the density of material flowing through a conduit. Mass
flowmeters also may
provide an indication of the density of the material.
Coriolis-type mass flowmeters are based on the well-known Coriolis effect, in
which
material flowing through a rotating conduit becomes a radially traveling mass
that is
affected by a Coriolis force and therefore experiences an acceleration. Many
Coriolis-type
mass flowmeters induce a Coriolis force by sinusoidally oscillating a conduit
about a pivot
axis orthogonal to the length of the conduit. In such mass flowmeters, the
Coriolis reaction
force experienced by the traveling fluid mass is transferred to the conduit
itself and is
manifested as a deflection or offset of the conduit in the direction of the
Coriolis force
vector in the plane of rotation.
Energy is supplied to the conduit by a driving mechanism that applies a
periodic
force to oscillate the conduit. One type of driving mechanism is an
electromechanical driver
that imparts a force proportional to an applied voltage. In an oscillating
flovvmeter, the
applied voltage is periodic, and is generally sinusoidal. The period of the
input voltage is
chosen so that the motion of the conduit matches a resonant mode of vibration
of the
conduit. This reduces the energy needed to sustain oscillation. An oscillating
flowmeter
may use a feedback loop in which a sensor signal that carries
1
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
instantaneous frequency and phase information related to oscillation of the
conduit is
amplified and fed back to the conduit using the electromechanical driver.
SUMMARY
In one general aspect, a multi-phase process fluid is passed through a
vibratable
flowtube. Motion is induced in the vibratable flowtube. A first apparent
property of the
multi-phase process fluid based on the motion of the vibratable flowtube is
determined,
and an apparent intermediate value associated with the multi-phase process
fluid is
determined based on the first apparent property. A corrected intermediate
value is
determined based on a mapping between the apparent intermediate value and the
corrected intermediate value. A phase-specific property of a phase of the
multi-phase
process fluid is determined based on the corrected intermediate value.
Implementations may include one or more of the following features. The
mapping may be a neural network configured to determine an error in the
intermediate
value resulting from the presence of the multi-flow process fluid. The
apparent
intermediate value may be determined to be within a first defined region of
values prior to
determining the corrected intermediate value, and the corrected intermediate
value may
be determined to be within a second defined region of values prior to
determining the
phase-specific property of a phase of the multi-phase process fluid.
The multi-phase process fluid may be a wet gas. The multi-phase process fluid
may include a first phase and a second phase, the first phase may include a
non-gas fluid,
and the second phase may include a gas. The multi-phase process fluid may
include a
first phase including a first non-gas fluid, and a second phase including a
second non-gas
fluid, and a third phase including a gas.
Determining the first apparent property of the multi-phase process fluid may
include determining a second apparent property of the multi-phase process
fluid. The
first apparent property of the multiphase process fluid may be a mass flow
rate and the
second apparent property may be a density.
One or more measurements corresponding to an additional property of the
process
fluid may be received. The additional property of the multi-phase fluid may
include one
or more of a temperature of the multi-phase fluid, a pressure associated with
the multi-
phase fluid, and a watercut of the multi-phase fluid, and determining an
apparent
intermediate value associated with the multi-phase process fluid based on the
first
2
CA 02658810 2014-11-13
95353-4
apparent property may include determining the intermediate value based on the
first apparent
property and the additional property.
Determining an apparent intermediate value associated with the multi-phase
process
fluid based on the first apparent property may include determining a volume
fraction
associated with an amount of non-gas fluid in the multi-phase process fluid
and a volumetric
flow rate of the multi-phase fluid. Determining an apparent intermediate value
associated with
the multi-phase process fluid based on the first apparent property may include
determining a
first Froude number corresponding to a non-gas phase of the multi-phase fluid
and a second
Froude number corresponding to a gas phase of the multi-phase fluid.
Determining a phase-specific property of a phase of the multi-phase process
fluid
based on the corrected intermediate value may include determining a mass flow
rate of a non-
gas phase of the multi-phase fluid.
Implementations of any of the techniques described above may include a method
or
process, a system, a flowmeter, or instructions stored on a storage device of
flowmeter
transmitter.
In an aspect, there is provided a method comprising: passing a multi-phase
process
fluid through a vibratable flowtube; inducing motion in the vibratable
flowtube; determining a
first apparent property of the multi-phase process fluid based on the motion
of the vibratable
flowtube; determining an apparent intermediate value associated with the multi-
phase process
fluid based on the first apparent property; mapping the apparent intermediate
value to a
corrected intermediate value; and determining a phase-specific property of a
phase of the
multi-phase process fluid based on the corrected intermediate value.
In another aspect, there is provided a flowmeter comprising: a vibratable
tlowtube, the
flowtube containing a multi-phase fluid; a driver connected to the flowtube
and configured to
impart motion to the flowtube such that the flowtube vibrates; a sensor
connected to the
flowtube and configured to sense the motion of the flowtube and generate a
sensor signal; and
a controller to receive the sensor signal and configured to: determine a first
apparent property
of the multi-phase process fluid based on the motion of the vibratable
flowtube; determine an
apparent intermediate value associated with the multi-phase process fluid
based on the first
apparent property; map the apparent intermediate value to a corrected
intermediate value; and
3
CA 02658810 2014-11-13
95353-4
determine a phase-specific property of a phase of the multi-phase process
fluid based on the
corrected intermediate value.
In a further aspect, there is provided a flowmeter transmitter comprising: at
least one
processing device; and a storage device, the storage device storing
instructions for causing the
at least one processing device to: determine a first apparent property of the
multi-phase
process fluid based on the motion of the vibratable tlowtube; determine an
apparent
intermediate value associated with the multi-phase process fluid based on the
first apparent
property; map the apparent intermediate value to a corrected intermediate
value; and
determine a phase-specific property of a phase of the multi-phase process
fluid based on the
corrected intermediate value.
The details of particular implementations are set forth in the accompanying
drawings
and description below. Other features will be apparent from the following
description,
including the drawings, and the claims.
DESCRIPTION OF DRAWINGS
Fig. 1 is a block diagram of a digital mass flowmeter.
Figs. 2A and 2B are perspective and side views of mechanical components of a
mass
flowmeter.
Figs. 3A-3C are schematic representations of three modes of motion of the
flowmeter
of Fig. 1.
Fig. 4 is a block diagram of an analog control and measurement circuit.
Fig. 5 is a block diagram of a digital mass flowmeter,
Fig. 6 is a flow chart showing operation of the meter of Fig. 5.
Figs. 7A and 7B are graphs of sensor data.
3a
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
Figs. 8A and 8B are graphs of sensor voltage relative to time.
Fig. 9 is a flow chart of a curve fitting procedure.
Fig. 10 is a flow chart of a procedure for generating phase differences.
Figs. 11A-11D, 12A-12D, and 13A-13D illustrate drive and sensor voltages at
system startup.
Fig. 14 is a flow chart of a procedure for measuring frequency, amplitude, and
phase of sensor data using a synchronous modulation technique.
Figs. 15A and 15B are block diagrams of a mass flowmeter.
Fig. 16 is a flow chart of a procedure implemented by the meter of Figs. 15A
and
15B.
Fig. 17 illustrates log-amplitude control of a transfer function.
Fig. 18 is a root locus diagram.
Figs. 19A-19D are graphs of analog-to-digital converter performance relative
to
temperature.
Figs. 20A-20C are graphs of phase measurements.
Figs. 21A and 21B are graphs of phase measurements.
Fig. 22 is a flow chart of a zero offset compensation procedure.
Figs. 23A-23C, 24A, and 24B are graphs of phase measurements.
Fig. 25 is a graph of sensor voltage.
Fig. 26 is a flow chart of a procedure for compensating for dynamic effects.
Figs. 27A-35E are graphs illustrating application of the procedure of Fig. 29.
Figs. 36A-36L are graphs illustrating phase measurement.
Fig. 37A is a graph of sensor voltages.
Figs. 37B and 37C are graphs of phase and frequency measurements
corresponding to the sensor voltages of Fig. 37A.
Figs. 37D and 37E are graphs of correction parameters for the phase and
frequency measurements of Figs. 37B and 37C.
Figs. 38A-38H are graphs of raw measurements.
Figs. 39A-39H are graphs of corrected measurements.
Figs. 40A-40H are graphs illustrating correction for aeration.
Fig. 41 is a block diagram illustrating the effect of aeration in a conduit.
Fig. 42 is a flow chart of a setpoint control procedure.
Figs. 43A-43C are graphs illustrating application of the procedure of Fig. 41.
Fig. 44 is a graph comparing the performance of digital and analog flowmeters.
4
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
Fig. 45 is a flow chart showing operation of a self-validating meter.
Fig. 46 is a block diagram of a two-wire digital mass flowmeter.
Figs. 47A-47C are graphs showing response of the digital mass flowmeter under
wet and empty conditions.
Fig. 48A is a chart showing results for batching from empty trials.
Fig. 48B is a diagram showing an experimental flow rig.
Fig. 49 is a graph showing mass-flow errors against drop in apparent density.
Fig. 50 is a graph showing residual mass-flow errors after applying
corrections.
Fig. 51 is a graph showing on-line response of the self-validating digital
mass
flowmeter to the onset of two-phase flow.
Fig. 52 is a block diagram of a digital controller implementing a neural
network
processor that may be used with the digital mass flowmeter.
Fig. 53 is a flow diagram showing the technique for implementing the neural
network to predict the mass-flow error and generate an error correction factor
to correct
the mass-flow measurement signal when two-phase flow is detected.
Fig. 54 is a 3D graph showing damping changes under two-phase flow conditions.
Fig. 55 is a flow diagram illustrating the experimental flow rig.
Fig. 56 is a 3D graph showing true mass flow error under two-phase flow
conditions.
Fig. 57 is a 3D graph showing corrected mass flow error under two-phase flow
conditions.
Fig. 58 is a graph comparing the uncorrected mass flow measurement signal with
the neural network corrected mass flow measurement signal.
Fig. 59 is a flow chart of a procedure for compensating for error under multi-
phase flow conditions.
Fig. 60 is a block diagram of a digital controller implementing a neural
network
processor that may be used with the digital mass flowmeter for multiple-phase
fluid
flows.
Fig. 61 is a flow diagram showing the technique for implementing the neural
network to predict the mass-flow error and generate an error correction factor
to correct
the mass-flow measurement signal when multiple-phase flows are expected and/or
detected.
Fig. 62 is a graphical view of a test matrix for wellheads tested based on
actual
testing at various well pressures and gas velocities.
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
Fig. 63 is a graphical view of raw density errors at various liquid void
fraction
percentages and of wells at various velocities and pressures.
Fig. 64 is a graphical view of raw mass flow errors at various liquid void
fraction
percentages and of wells at various velocities and pressures.
Fig. 65 is a graphical view of raw liquid void fraction errors of wells at
various
velocities and pressures.
Fig. 66 is a graphical view of raw volumetric flow errors for wells at various
velocities and pressures.
Fig. 67 is a graphical view of corrected liquid void fractions of wells at
various
velocities and pressures.
Fig. 68 is a graphical view of corrected mixture volumetric flow of wells at
various velocities and pressures.
Fig. 69 is a graphical view of corrected gas mass flow of wells at various
velocities and pressures.
Fig. 70 is a graphical view of corrected gas cumulative probability of the
digital
flowmeter tested.
Fig. 71 is a graphical view of corrected liquid mass flow error of wells at
various
velocities and pressures.
Fig. 72 is a graphical view of corrected gas cumulative probability of the
digital
flowmeter tested.
DETAILED DESCRIPTION
Techniques are provided for accounting for the effects of multi-phase flow in,
for example, a digital flowmeter. The multi-phase flow may be, for example, a
two-phase
flow or a three-phase flow. In general, a two-phase flow is a fluid that
includes two
phases or components. For example, a two-phase flow may include a phase that
includes
a non-gas fluid (such as a liquid) and a phase that includes a gas. A three-
phase flow is a
fluid that includes three phases. For example, a three-phase flow may be a
fluid with a
gas phase and two non-gas liquids. For example, a three-phase flow may include
natural
gas, oil, and water. A two-phase flow may include, for example, natural gas
and oil.
6
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
Although the digital flowmeter continues to operate in the presence of a multi-
phase fluid, any propertires (e.g., the mass flow rate and density of the
multi-phase fluid)
determined by the digital flowmeter may be inaccurate because the
determination of these
properties using conventional techniques is generally based on an assumption
that the
fluid flowing through the flowmeter is single-phase. Thus, even though the
fluid is not a
single phase fluid, the flowmeter may continue to operate and generate
apparent values of
properties such as the mass flow rate and density of the multi-phase fluid. As
described
below with respect to Figs. 59-72, in one implementation parameters such as
mass flow
rate and density of each of the phases of the multi-phase flow may be
determined from
the apparent mass flow rate and the apparent density of the multi-phase fluid.
In
particular, and as discussed in more detail below, in one implementation, one
or more
intermediate values, such as the liquid volume fraction and the volumetric
flowrate or gas
and non-gas Froude numbers, are determined from the apparent mass flow rate
and
apparent density of the multi-phase fluid and the intermediate value(s) may be
corrected
to account for the presence of multiple phases in the fluid using a neural
network or other
mapping. The mass flow rate and density of each phase of the multi-phase fluid
may be
determined from the corrected intermediate value(s). Using the intermediate
value(s)
rather than the mass flow rate and density of the multi-phase fluid may help
to improve
the accuracy of the mass flow rate and density of each phase of the multi-
phase fluid.
Before the techniques are described starting with reference to Fig. 59,
digital
flowmeters are discussed with reference to Figs. 1 ¨ 39. Various techniques
for
accounting for the effects of multi-phase flow in, for example, a digital
flowmeter are
discussed starting with Fig. 40.
Referring to Fig. 1, a digital mass flowmeter 100 includes a digital
controller 105,
one or more motion sensors 110, one or more drivers 115, a conduit 120 (also
referred to
as a flowtube), and a temperature sensor 125. The digital controller 105 may
be
implemented using one or more of, for example, a processor, a field-
programmable gate
array, an ASIC, other programmable logic or gate arrays, or programmable logic
with a
processor core. The digital controller generates a measurement of mass flow
through the
conduit 120 based at least on signals received from the motion sensors 110.
The digital
controller also controls the drivers 115 to induce motion in the conduit 120.
This motion
is sensed by the motion sensors 110.
Mass flow through the conduit 120 is related to the motion induced in the
conduit in response to a driving force supplied by the drivers 115. In
particular, mass
7
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
flow is related to the phase and frequency of the motion, as well as to the
temperature of
the conduit. The digital mass flowmeter also may provide a measurement of the
density
of material flowing through the conduit. The density is related to the
frequency of the
motion and the temperature of the conduit. Many of the described techniques
are
applicable to a densitometer that provides a measure of density rather than a
measure of
mass flow.
The temperature in the conduit, which is measured using the temperature
sensor 125, affects certain properties of the conduit, such as its stiffness
and dimensions.
The digital controller compensates for these temperature effects. The
temperature of the
digital controller 105 affects, for example, the operating frequency of the
digital
controller. In general, the effects of controller temperature are sufficiently
small to be
considered negligible. However, in some instances, the digital controller may
measure
the controller temperature using a solid state device and may compensate for
effects of
the controller temperature.
A. Mechanical Design
In one implementation, as illustrated in Figs. 2A and 2B, the conduit 120 is
designed to be inserted in a pipeline (not shown) having a small section
removed or
reserved to make room for the conduit. The conduit 120 includes mounting
flanges 12 for
connection to the pipeline, and a central manifold block 16 supporting two
parallel planar
loops 18 and 20 that are oriented perpendicularly to the pipeline. An
electromagnetic
driver 46 and a sensor 48 are attached between each end of loops 18 and 20.
Each of the
two drivers 46 corresponds to a driver 115 of Fig. 1, while each of the two
sensors 48
corresponds to a sensor 110 of Fig. 1.
The drivers 46 on opposite ends of the loops are energized with current of
equal
magnitude but opposite sign (i.e., currents that are 180 out-of-phase) to
cause straight
sections 26 of the loops 18, 20 to rotate about their co-planar perpendicular
bisector 56,
which intersects the tube at point P (Fig. 2B). Repeatedly reversing (e.g.,
controlling
sinusoidally) the energizing current supplied to the drivers causes each
straight section 26
to undergo oscillatory motion that sweeps out a bow tie shape in the
horizontal plane
about line 56-56, the axis of symmetry of the loop. The entire lateral
excursion of the
loops at the lower rounded turns 38 and 40 is small, on the order of 1/16 of
an inch for a
8
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
two foot long straight section 26 of a pipe having a one inch diameter. The
frequency of
oscillation is typically about 80 to 90 Hertz.
B. Conduit Motion
The motion of the straight sections of loops 18 and 20 are shown in three
modes in
Figs. 3A, 3B and 3C. In the drive mode shown in Fig. 3B, the loops are driven
1800 out-
of-phase about their respective points P so that the two loops rotate
synchronously but in
the opposite sense. Consequently, respective ends such as A and C periodically
come
together and go apart.
The drive motion shown in Fig. 3B induces the Coriolis mode motion shown in
Fig. 3A, which is in opposite directions between the loops and moves the
straight sections
26 slightly toward (or away) from each other. The Coriolis effect is directly
related to
mvW, where m is the mass of material in a cross section of a loop, v is the
velocity at
which the mass is moving (the volumetric flow rate), W is the angular velocity
of the loop
(W = W sin cot), and my is the mass flow rate. The Coriolis effect is
greatest when the two
straight sections are driven sinusoidally and have a sinusoidally varying
angular velocity.
Under these conditions, the Coriolis effect is 90 out-of-phase with the drive
signal.
Fig. 3C shows an undesirable common mode motion that deflects the loops in the
same direction. This type of motion might be produced by an axial vibration in
the
pipeline in the example of Figs. 2A and 2B because the loops are perpendicular
to the
pipeline.
The type of oscillation shown in Fig. 3B is called the antisymmetrical mode,
and
the Coriolis mode of Fig. 3A is called the symmetrical mode. The natural
frequency of
oscillation in the antisymmetrical mode is a function of the torsional
resilience of the legs.
Ordinarily the resonant frequency of the antisymmetrical mode for conduits of
the shape
shown in Figs. 2A and 2B is higher than the resonant frequency of the
symmetrical mode.
To reduce the noise sensitivity of the mass flow measurement, it is desirable
to maximize
the Coriolis force for a given mass flow rate. As noted above, the loops are
driven at their
resonant frequency, and the Coriolis force is directly related to the
frequency at which the
loops are oscillating (i.e., the angular velocity of the loops). Accordingly,
the loops are
driven in the antisymmetrical mode, which tends to have the higher resonant
frequency.
Other implementations may include different conduit designs. For example, a
single loop or a straight tube section may be employed as the conduit.
9
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
C. Electronic Design
The digital controller 105 determines the mass flow rate by processing signals
produced by the sensors 48 (i.e., the motion sensors 110) located at opposite
ends of the
loops. The signal produced by each sensor includes a component corresponding
to the
relative velocity at which the loops are driven by a driver positioned next to
the sensor
and a component corresponding to the relative velocity of the loops due to
Coriolis forces
induced in the loops. The loops are driven in the antisymmetrical mode, so
that the
components of the sensor signals corresponding to drive velocity are equal in
magnitude
but opposite in sign. The resulting Coriolis force is in the symmetrical mode
so that the
components of the sensor signals corresponding to Coriolis velocity are equal
in
magnitude and sign. Thus, differencing the signals cancels out the Coriolis
velocity
components and results in a difference that is proportional to the drive
velocity.
Similarly, summing the signals cancels out the drive velocity components and
results in a
sum that is proportional to the Coriolis velocity, which, in turn, is
proportional to the
Coriolis force. This sum then may be used to determine the mass flow rate.
1. Analog Control System
The digital mass flowmeter 100 provides considerable advantages over
traditional,
analog mass flowmeters. For use in later discussion, Fig. 4 illustrates an
analog control
system 400 of a traditional mass flowmeter. The sensors 48 each produce a
voltage
signal, with signal VA0 being produced by sensor 48a and signal Vim being
produced by
sensor 48b. VA0 and VB0 correspond to the velocity of the loops relative to
each other at
the positions of the sensors. Prior to processing, signals VA() and VB0 are
amplified at
respective input amplifiers 405 and 410 to produce signals VA! and Vm. To
correct for
imbalances in the amplifiers and the sensors, input amplifier 410 has a
variable gain that
is controlled by a balance signal coming from a feedback loop that contains a
synchronous demodulator 415 and an integrator 420.
At the output of amplifier 405, signal VAi is of the form:
VA1 = VD sin cot + Vc cos wt,
and, at the output of amplifier 410, signal Vm is of the form:
¨VD sin cot + Vc cos cot,
VB1
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
where VD and Vc are, respectively, the drive voltage and the Coriolis voltage,
and co is
the drive mode angular frequency.
Voltages VA' and VB1 are differenced by operational amplifier 425 to produce:
VDRV = VA1 VB1 = 2VD sin cot,
where VDRV corresponds to the drive motion and is used to power the drivers.
In addition
to powering the drivers, VDRV is supplied to a positive going zero crossing
detector 430
that produces an output square wave FDRv having a frequency corresponding to
that of
VDRV (CD = 27CFDRV). FDRV is used as the input to a digital phase locked loop
circuit 435.
FDRv also is supplied to a processor 440.
Voltages VAI and Vgi are summed by operational amplifier 445 to produce:
VCOR = VA1 VB1 = 2Vc cos cot,
where VcoR is related to the induced Coriolis motion.
Vc0R is supplied to a synchronous demodulator 450 that produces an output
voltage Vm that is directly proportional to mass by rejecting the components
of VcoR that
do not have the same frequency as, and are not in phase with, a gating signal
Q. The
phase locked loop circuit 435 produces Q, which is a quadrature reference
signal that has
the same frequency (co) as VDRV and is 90 out of phase with VDRy (i.e., in
phase with
VcoR). Accordingly, synchronous demodulator 450 rejects frequencies other than
o) so
that Vm corresponds to the amplitude of VcOR at co. This amplitude is directly
proportional to the mass in the conduit.
Vm is supplied to a voltage-to-frequency converter 455 that produces a square
wave signal Fm having a frequency that corresponds to the amplitude of Vm. The
processor 440 then divides Fm by FDRv to produce a measurement of the mass
flow rate.
Digital phase locked loop circuit 435 also produces a reference signal I that
is in
phase with VDRV and is used to gate the synchronous demodulator 415 in the
feedback
loop controlling amplifier 410. When the gains of the input amplifiers 405 and
410
multiplied by the drive components of the corresponding input signals are
equal, the
summing operation at operational amplifier 445 produces zero drive component
(i.e., no
signal in phase with VDRV) in the signal VcOR. When the gains of the input
amplifiers 405
and 410 are not equal, a drive component exists in VcOR. This drive component
is
extracted by synchronous demodulator 415 and integrated by integrator 420 to
generate
an error voltage that corrects the gain of input amplifier 410. When the gain
is too high
or too low, the synchronous demodulator 415 produces an output voltage that
causes the
11
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
integrator to change the error voltage that modifies the gain. When the gain
reaches the
desired value, the output of the synchronous modulator goes to zero and the
error voltage
stops changing to maintain the gain at the desired value.
2. Digital Control System
Fig. 5 provides a block diagram of an implementation 500 of the digital mass
flowmeter 100 that includes the conduit 120, drivers 46, and sensors 48 of
Figs. 2A and
2B, along with a digital controller 505. Analog signals from the sensors 48
are converted
to digital signals by analog-to-digital ("A/D") converters 510 and supplied to
the
controller 505. The A/D converters may be implemented as separate converters,
or as
separate channels of a single converter.
Digital-to-analog ("D/A") converters 515 convert digital control signals from
the
controller 505 to analog signals for driving the drivers 46. The use of a
separate drive
signal for each driver has a number of advantages. For example, the system may
easily
switch between symmetrical and antisymmetrical drive modes for diagnostic
purposes. In
other implementations, the signals produced by converters 515 may be amplified
by
amplifiers prior to being supplied to the drivers 46. In still other
implementations, a
single D/A converter may be used to produce a drive signal applied to both
drivers, with
the drive signal being inverted prior to being provided to one of the drivers
to drive the
conduit 120 in the antisymmetrical mode.
High precision resistors 520 and amplifiers 525 are used to measure the
current
supplied to each driver 46. A/D converters 530 convert the measured current to
digital
signals and supply the digital signals to controller 505. The controller 505
uses the
measured currents in generating the driving signals.
Temperature sensors 535 and pressure sensors 540 measure, respectively, the
temperature and the pressure at the inlet 545 and the outlet 550 of the
conduit. A/D
converters 555 convert the measured values to digital signals and supply the
digital
signals to the controller 505. The controller 505 uses the measured values in
a number of
ways. For example, the difference between the pressure measurements may be
used to
determine a back pressure in the conduit. Since the stiffness of the conduit
varies with
the back pressure, the controller may account for conduit stiffness based on
the
determined back pressure.
An additional temperature sensor 560 measures the temperature of the crystal
oscillator 565 used by the A/D converters. An A/D converter 570 converts this
12
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
temperature measurement to a digital signal for use by the controller 505. The
input/output relationship of the AID converters varies with the operating
frequency of the
converters, and the operating frequency varies with the temperature of the
crystal
oscillator. Accordingly, the controller uses the temperature measurement to
adjust the
data provided by the AID converters, or in system calibration.
In the implementation of Fig. 5, the digital controller 505 processes the
digitized
sensor signals produced by the A/D converters 510 according to the procedure
600
illustrated in Fig. 6 to generate the mass flow measurement and the drive
signal supplied
to the drivers 46. Initially, the controller collects data from the sensors
(step 605). Using
this data, the controller determines the frequency of the sensor signals (step
610),
eliminates zero offset from the sensor signals (step 615), and determines the
amplitude
(step 620) and phase (step 625) of the sensor signals. The controller uses
these calculated
values to generate the drive signal (step 630) and to generate the mass flow
and other
measurements (step 635). After generating the drive signals and measurements,
the
controller collects a new set of data and repeats the procedure. The steps of
the procedure
600 may be performed serially or in parallel, and may be performed in varying
order.
Because of the relationships between frequency, zero offset, amplitude, and
phase,
an estimate of one may be used in calculating another. This leads to repeated
calculations
to improve accuracy. For example, an initial frequency determination used in
determining the zero offset in the sensor signals may be revised using offset-
eliminated
sensor signals. In addition, where appropriate, values generated for a cycle
may be used
as starting estimates for a following cycle.
a. Data Collection
For ease of discussion, the digitized signals from the two sensors will be
referred
to as signals SVi and SV2, with signal SVi coming from sensor 48a and signal
SV2
coming from sensor 48b. Although new data is generated constantly, it is
assumed that
calculations are based upon data corresponding to one complete cycle of both
sensors.
With sufficient data buffering, this condition will be true so long as the
average time to
process data is less than the time taken to collect the data. Tasks to be
carried out for a
cycle include deciding that the cycle has been completed, calculating the
frequency of the
cycle (or the frequencies of SVI and SV2), calculating the amplitudes of SVi
and SV2,
and calculating the phase difference between SVI and SV2. In some
implementations,
these calculations are repeated for each cycle using the end point of the
previous cycle as
13
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
the start for the next. In other implementations, the cycles overlap by 1800
or other
amounts (e.g., 90 ) so that a cycle is subsumed within the cycles that precede
and follow
it.
Figs. 7A and 7B illustrate two vectors of sampled data from signals SVi and
SV2,
which are named, respectively, svl_in and sv2_in. The first sampling point of
each
vector is known, and corresponds to a zero crossing of the sine wave
represented by the
vector. For svl_in, the first sampling point is the zero crossing from a
negative value to a
positive value, while for sv2_in the first sampling point is the zero crossing
from a
positive value to a negative value.
An actual starting point for a cycle (i.e., the actual zero crossing) will
rarely
coincide exactly with a sampling point. For this reason, the initial sampling
points
(start_sample_SV1 and start_sample_SV2) are the sampling points occurring just
before
the start of the cycle. To account for the difference between the first
sampling point and
the actual start of the cycle, the approach also uses the position
(start_offset SV1 or
start offset SV2) between the starting sample and the next sample at which the
cycle
actually begins.
Since there is a phase offset between signals SVI and 5V2, svl_in and sv2_in
may
start at different sampling points. If both the sample rate and the phase
difference are
high, there may be a difference of several samples between the start of svl_in
and the
start of sv2_in. This difference provides a crude estimate of the phase
offset, and may be
used as a check on the calculated phase offset, which is discussed below. For
example,
when sampling at 55 kHz, one sample corresponds to approximately 0.5 degrees
of phase
shift, and one cycle corresponds to about 800 sample points.
When the controller employs functions such as the sum (A+B) and difference (A-
B), with B weighted to have the same amplitude as A, additional variables
(e.g.,
start_sample_sum and start_offset_sum) track the start of the period for each
function.
The sum and difference functions have a phase offset halfway between SVi and
SV2.
In one implementation, the data structure employed to store the data from the
sensors is a circular list for each sensor, with a capacity of at least twice
the maximum
number of samples in a cycle. With this data structure, processing may be
carried out on
data for a current cycle while interrupts or other techniques are used to add
data for a
following cycle to the lists.
14
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
Processing is performed on data corresponding to a full cycle to avoid errors
when
employing approaches based on sine-waves. Accordingly, the first task in
assembling
data for a cycle is to determine where the cycle begins and ends. When
nonoverlapping
cycles are employed, the beginning of the cycle may be identified as the end
of the
previous cycle. When overlapping cycles are employed, and the cycles overlap
by 1800,
the beginning of the cycle may be identified as the midpoint of the previous
cycle, or as
the endpoint of the cycle preceding the previous cycle.
The end of the cycle may be first estimated based on the parameters of the
previous cycle and under the assumption that the parameters will not change by
more than
a predetermined amount from cycle to cycle. For example, five percent may be
used as
the maximum permitted change from the last cycle's value, which is reasonable
since, at
sampling rates of 55 kHz, repeated increases or decreases of five percent in
amplitude or
frequency over consecutive cycles would result in changes of close to 5,000
percent in
one second.
By designating five percent as the maximum permissible increase in amplitude
and frequency, and allowing for a maximum phase change of 50 in consecutive
cycles, a
conservative estimate for the upper limit on the end of the cycle for signal
SVI may be
determined as:
365* sample_rate
end_sample SV1 start¨ sample + __
360 est_freq * 0.95
where start_sample_SV1 is the first sample of svl_in, sample_rate is the
sampling rate,
and est_freq is the frequency from the previous cycle. The upper limit on the
end of the
cycle for signal SV2 (end_sample_SV2) may be determined similarly.
After the end of a cycle is identified, simple checks may be made as to
whether
the cycle is worth processing. A cycle may not be worth processing when, for
example,
the conduit has stalled or the sensor waveforms are severely distorted.
Processing only
suitable cycles provides considerable reductions in computation.
One way to determine cycle suitability is to examine certain points of a cycle
to
confirm expected behavior. As noted above, the amplitudes and frequency of the
last
cycle give useful starting estimates of the corresponding values for the
current cycle.
Using these values, the points corresponding to 30 , 150 , 210 and 330 of
the cycle may
be examined. If the amplitude and frequency were to match exactly the
amplitude and
frequency for the previous cycle, these points should have values
corresponding to
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
est_amp/2, est_amp/2,
-est_amp/2, and -est_amp/2, respectively, where est_amp is the estimated
amplitude of a
signal (i.e., the amplitude from the previous cycle). Allowing for a five
percent change in
both amplitude and frequency, inequalities may be generated for each quarter
cycle. For
the 30 point, the inequality is
(
te
svlin start_sample_SV1 + ____ 30 * sample _ra > 0.475* est_amp SV
_
360 est_freq *1.05
The inequalities for the other points have the same form, with the degree
offset term
(x/360) and the sign of the est_amp SV1 term having appropriate values. These
inequalities can be used to check that the conduit has vibrated in a
reasonable manner.
Measurement processing takes place on the vectors svl_in(startend) and
sv2_in(start:end) where:
start = min (start_sample_SV1, start_sample_SV2), and
end = max (end_sample_SV1, end_sample SV2).
The difference between the start and end points for a signal is indicative of
the frequency
of the signal.
b. Frequency Determination
The frequency of a discretely-sampled pure sine wave may be calculated by
detecting the transition between periods (i.e., by detecting positive or
negative zero-
crossings) and counting the number of samples in each period. Using this
method,
sampling, for example, an 82.2 Hz sine wave at 55 kHz will provide an estimate
of
frequency with a maximum error of 0.15 percent. Greater accuracy may be
achieved by
estimating the fractional part of a sample at which the zero-crossing actually
occurred
using, for example, start_offset_SV1 and start offset_SV2. Random noise and
zero
offset may reduce the accuracy of this approach.
As illustrated in Figs. 8A and 8B, a more refined method of frequency
determination uses quadratic interpolation of the square of the sine wave.
With this
method, the square of the sine wave is calculated, a quadratic function is
fitted to match
the minimum point of the squared sine wave, and the zeros of the quadratic
function are
used to determine the frequency. If
svt = A sin xt +6+ cyst,
16
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
where svt is the sensor voltage at time t, A is the amplitude of oscillation,
xt is the radian
angle at time t (i.e., xt = 27-cft), 6 is the zero offset, ct is a random
variable with
distribution N(0,1), and a is the variance of the noise, then the squared
function is given
by:
svt2 = A2 sin2 xe + 2A(8+ act ) sin xt + Mae, +82 +0_2g.
When xt is close to 27c, sin xt and sin2xt can be approximated as xot = xt-27t
and x0t2,
respectively. Accordingly, for values of xt close to 27c, at can be
approximated as:
a 2t A2x,,2,t + 2A(S+ o-et)xot + 2cSo-e, +82 +a26
=-'t% (A2X02t 2A8x ot + 82 )+ 0-et (2AX0, + + act ).
This is a pure quadratic (with a non-zero minimum, assuming 6=0) plus noise,
with the
amplitude of the noise being dependent upon both a and 6. Linear interpolation
also
could be used.
Error sources associated with this curve fitting technique are random noise,
zero
offset, and deviation from a true quadratic. Curve fitting is highly sensitive
to the level of
random noise. Zero offset in the sensor voltage increases the amplitude of
noise in the
sine-squared function, and illustrates the importance of zero offset
elimination (discussed
below). Moving away from the minimum, the square of even a pure sine wave is
not
entirely quadratic. The most significant extra term is of fourth order. By
contrast, the
most significant extra term for linear interpolation is of third order.
Degrees of freedom associated with this curve fitting technique are related to
how
many, and which, data points are used. The minimum is three, but more may be
used (at
greater computational expense) by using least-squares fitting. Such a fit is
less
susceptible to random noise. Fig. 8A illustrates that a quadratic
approximation is good up
to some 200 away from the minimum point. Using data points further away from
the
minimum will reduce the influence of random noise, but will increase the
errors due to
the non-quadratic terms (i.e., fourth order and higher) in the sine-squared
function.
Fig. 9 illustrates a procedure 900 for performing the curve fitting technique.
As a
first step, the controller initializes variables (step 905). These variables
include
end_point, the best estimate of the zero crossing point; ep_int, the integer
value nearest to
end_point; s[0..i], the set of all sample points; z[lci, the square of the
sample point closest
to end_point; z[0..n-1], a set of squared sample points used to calculate
end_point; n, the
17
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
number of sample points used to calculate end_point (n = 2k + 1); step_length,
the
number of samples in s between consecutive values in z; and iteration_count, a
count of
the iterations that the controller has performed.
The controller then generates a first estimate of end_point (step 910). The
controller generates this estimate by calculating an estimated zero-crossing
point based on
the estimated frequency from the previous cycle and searching around the
estimated
crossing point (forwards and backwards) to find the nearest true crossing
point (i.e., the
occurrence of consecutive samples with different signs). The controller then
sets
end_point equal to the sample point having the smaller magnitude of the
samples
surrounding the true crossing point.
Next, the controller sets n, the number of points for curve fitting (step
915). The
controller sets n equal to 5 for a sample rate of 11 kHz, and to 21 for a
sample rate of 44
kHz. The controller then sets iteration_count to 0 (step 920) and increments
iteration_count (step 925) to begin the iterative portion of the procedure.
As a first step in the iterative portion of the procedure, the controller
selects
step_length (step 930) based on the value of iteration_count. The controller
sets
step_length equal to 6, 3, or 1 depending on whether iteration_count equals,
respectively,
1, 2 or 3.
Next, the controller determines ep_int as the integer portion of the sum of
end_point and 0.5 (step 935) and fills the z array (step 940). For example,
when n equals
5, z[0] = s[ep int-2*step Jength]2, z[1] = s[ep int-step_length]2, z[2] =
s[ep_int]2,
z[3] = s[ep_int+step_length]2, and z[4] = s[ep_int+2*step Jength]2.
Next, the controller uses a filter, such as a Savitzky-Golay filter, to
calculate
smoothed values of z[k-1], z[k] and z[k+1] (step 945). Savitzky-Golay
smoothing filters
are discussed by Press et al. in Numerical Recipes in C, pp. 650-655 (2nd ed.,
Cambridge
University Press, 1995), which is incorporated by reference. The controller
then fits a
quadratic to z[k-1], z[k] and z[k+1] (step 950), and calculates the minimum
value of the
quadratic (z*) and the corresponding position (x*) (step 955).
If x* is between the points corresponding to k-1 and k+1 (step 960), then the
controller sets end_point equal to x* (step 965). Thereafter, if
iteration_count is less than
3 (step 970), the controller increments iteration_count (step 925) and repeats
the iterative
portion of the procedure.
If x* is not between the points corresponding to k-1 and k+1 (step 960), or if
iteration_count equals 3 (step 970), the controller exits the iterative
portion of the
18
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
procedure. The controller then calculates the frequency based on the
difference between
end_point and the starting point for the cycle, which is known (step 975).
In essence, the procedure 900 causes the controller to make three attempts to
home in on end_point, using smaller step lengths in each attempt. If the
resulting
minimum for any attempt falls outside of the points used to fit the curve
(i.e., there has
been extrapolation rather than interpolation), this indicates that either the
previous or new
estimate is poor, and that a reduction in step size is unwarranted.
The procedure 900 may be applied to at least three different sine waves
produced
by the sensors. These include signals SV1 and SV2 and the weighted sum of the
two.
Moreover, assuming that zero offset is eliminated, the frequency estimates
produced for
these signals are independent. This is clearly true for signals SV1 and SV2,
as the errors
on each are independent. It is also true, however, for the weighted sum, as
long as the
mass flow and the corresponding phase difference between signals SV1 and SV2
are large
enough for the calculation of frequency to be based on different samples in
each case.
When this is true, the random errors in the frequency estimates also should be
independent.
The three independent estimates of frequency can be combined to provide an
improved estimate. This combined estimate is simply the mean of the three
frequency
estimates.
c. Zero Offset Compensation
An important error source in a Coriolis transmitter is zero offset in each of
the
sensor voltages. Zero offset is introduced into a sensor voltage signal by
drift in the pre-
amplification circuitry and the analog-to-digital converter. The zero offset
effect may be
worsened by slight differences in the pre-amplification gains for positive and
negative
voltages due to the use of differential circuitry. Each error source varies
between
transmitters, and will vary with transmitter temperature and more generally
over time
with component wear.
An example of the zero offset compensation technique employed by the
controller
is discussed in detail below. In general, the controller uses the frequency
estimate and an
integration technique to determine the zero offset in each of the sensor
signals. The
controller then eliminates the zero offset from those signals. After
eliminating zero offset
from signals SV1 and SV2, the controller may recalculate the frequency of
those signals to
provide an improved estimate of the frequency.
19
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
d. Amplitude Determination
The amplitude of oscillation has a variety of potential uses. These include
regulating conduit oscillation via feedback, balancing contributions of sensor
voltages
when synthesizing driver waveforms, calculating sums and differences for phase
measurement, and calculating an amplitude rate of change for measurement
correction
purposes.
In one implementation, the controller uses the estimated amplitudes of signals
SVi
and SV2 to calculate the sum and difference of signals SVi and SV2, and the
product of
the sum and difference. Prior to determining the sum and difference, the
controller
compensates one of the signals to account for differences between the gains of
the two
sensors. For example, the controller may compensate the data for signal SV2
based on the
ratio of the amplitude of signal SVI to the amplitude of signal SV2 so that
both signals
have the same amplitude.
The controller may produce an additional estimate of the frequency based on
the
calculated sum. This estimate may be averaged with previous frequency
estimates to
produce a refined estimate of the frequency of the signals, or may replace the
previous
estimates.
The controller may calculate the amplitude according to a Fourier-based
technique
to eliminate the effects of higher harmonics. A sensor voltage x(t) over a
period T (as
identified using zero crossing techniques) can be represented by an offset and
a series of
harmonic terms as:
x(t) = a0/2 + al cos(cot) + a2cos(2cot) + a3cos(3oot) + +
bisin(cot) + b2sin(2uot) +
With this representation, a non-zero offset a0 will result in non-zero cosine
terms an.
Though the amplitude of interest is the amplitude of the fundamental component
(i.e., the
amplitude at frequency co), monitoring the amplitudes of higher harmonic
components
(i.e., at frequencies kw, where k is greater than 1) may be of value for
diagnostic
purposes. The values of an and bn may be calculated as:
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
an = ¨2 s X(t)COS modt,
T
and
bn = ¨2 f T X(t) sin ncodt.
T
The amplitude, A,õ of each harmonic is given by:
= 2n + .
The integrals are calculated using Simpson's method with quadratic correction
(described
below). The chief computational expense of the method is calculating the pure
sine and
cosine functions.
e. Phase Determination
The controller may use a number of approaches to calculate the phase
difference
between signals 5V1 and SV2. For example, the controller may determine the
phase
offset of each harmonic, relative to the starting time at 1=0, as:
wn =tan
bn
The phase offset is interpreted in the context of a single waveform as being
the difference
between the start of the cycle (i.e., the zero-crossing point) and the point
of zero phase for
the component of SV(t) of frequency w. Since the phase offset is an average
over the
entire waveform, it may be used as the phase offset from the midpoint of the
cycle.
Ideally, with no zero offset and constant amplitude of oscillation, the phase
offset should
be zero every cycle. The controller may determine the phase difference by
comparing the
phase offset of each sensor voltage over the same time period.
The amplitude and phase may be generated using a Fourier method that
eliminates the effects of higher harmonics. This method has the advantage that
it does not
assume that both ends of the conduits are oscillating at the same frequency.
As a first
step in the method, a frequency estimate is produced using the zero crossings
to measure
the time between the start and end of the cycle. If linear variation in
frequency is
assumed, this estimate equals the time-averaged frequency over the period.
Using the
estimated, and assumed time-invariant, frequency co of the cycle, the
controller
calculates:
21
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
27r
Ii
-= 2c F SV(t) sin (cot) dt, and
71- o
2co
12 = ¨0 SV(t) cos (cot) dt,
7r
where SV(t) is the sensor voltage waveform (i.e., SVi(t) or SV2(t)). The
controller then
determines the estimates of the amplitude and phase:
Amp = + 122, and
Phase = tan'
The controller then calculates a phase difference, assuming that the average
phase and frequency of each sensor signal is representative of the entire
waveform. Since
these frequencies are different for SVI and SV2, the corresponding phases are
scaled to
the average frequency. In addition, the phases are shifted to the same
starting point (i.e.,
the midpoint of the cycle on SV1). After scaling, they are subtracted to
provide the phase
difference:
av_freq
scaled_phase_SV, = phase_SV, ___________
freq_SV,
(midpoint_S V2 - midpoint_SV, )11 freq_SV2
scaled shift SV2 = ______________________________________ , and
-
360
scaled_phase_SV2 = (phase SV2 + scale_shift_SV2 ) av_freq
freq_SV2
where h is the sample length and the midpoints are defined in terms of
samples:
midpoint SY = (startpoint _SVõ + endpoint SVx )
¨ ),
2
In general, phase and amplitude are not calculated over the same time-frame
for
the two sensors. When the flow rate is zero, the two cycle mid-points are
coincident.
However, they diverge at high flow rates so that the calculations are based on
sample sets
that are not coincident in time. This leads to increased phase noise in
conditions of
changing mass flow. At full flow rate, a phase shift of 4 (out of 360 ) means
that only
99% of the samples in the SVi and SV2 data sets are coincident. Far greater
phase shifts
may be observed under aerated conditions, which may lead to even lower rates
of overlap.
22
CA 02658810 2009-02-25
WO 2008/027908 PCT/US2007/077018
Fig. 10 illustrates a modified approach 1000 that addresses this issue. First,
the
controller finds the frequencies (f1, f2) and the mid-points (mi, m2) of the
SVi and SV2
data sets (d1, d2) (step 1005). Assuming linear shift in frequency from the
last cycle, the
controller calculates the frequency of SV2 at the midpoint of SVi (f2mi) and
the frequency
of SVI at the midpoint of SV2 (fh,2) (step 1010).
The controller then calculates the starting and ending points of new data sets
(dim2
and d2,,i) with mid-points m2 and m1 respectively, and assuming frequencies of
firn2 and
f2mi (step 1015). These end points do not necessarily coincide with zero
crossing points.
However, this is not a requirement for Fourier-based calculations.
The controller then carries out the Fourier calculations of phase and
amplitude on
the sets d1 and d2rni, and the phase difference calculations outlined above
(step 1020).
Since the mid-points of d1 and d2.1 are identical, scale-shift_SV2 is always
zero and can
be ignored. The controller repeats these calculations for the data sets d2 and
dim2 (step
1025). The controller then generates averages of the calculated amplitude and
phase
difference for use in measurement generation (step 1030). When there is
sufficient
separation between the mid points m1 and m2, the controller also may use the
two sets of
results to provide local estimates of the rates of change of phase and
amplitude.
The controller also may use a difference-amplitude method that involves
calculating a difference between SVi and SV2, squaring the calculated
difference, and
integrating the result. According to another approach, the controller
synthesizes a sine
wave, multiplies the sine wave by the difference between signals SVI and SV2,
and
integrates the result. The controller also may integrate the product of
signals SVI and
V2, which is a sine wave having a frequency 2f (where f is the average
frequency of
signals SVI and SV2), or may square the product and integrate the result. The
controller
also may synthesize a cosine wave comparable to the product sine wave and
multiply the
synthesized cosine wave by the product sine wave to produce a sine wave of
frequency 4f
that the controller then integrates. The controller also may use multiple ones
of these
approaches to produce separate phase measurements, and then may calculate a
mean
value of the separate measurements as the final phase measurement.
The difference-amplitude method starts with:
SVI (t) = Al sin 270 + and SV2(t) = A2 sin Inft ¨
2 2 j
23
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
where (p is the phase difference between the sensors. Basic trigonometric
identities may
be used to define the sum (Sum) and difference (Diff) between the signals as:
Sum a- ST/i(t)+ A1SV2(t) = 2A1 cos ¨(1) = sin 27rft, and
A2 2
Diff SV1(t) ¨ ¨AiSV2(t) = 2A1 sing.cos 270.
A2 2
These functions have amplitudes of 2Aicos((p2) and 2Aisin(cp2), respectively.
The
controller calculates data sets for Sum and Diff from the data for SVI and
SV2, and then
uses one or more of the methods described above to calculate the amplitude of
the signals
represented by those data sets. The controller then uses the calculated
amplitudes to
calculate the phase difference, (p.
As an alternative, the phase difference may be calculated using the function
Prod, defined as:
Prod a SumxDiff = 4A? cos ¨ sing. cos2aft = sin 27th
2 2
= A?simp sin 47rft ,
which is a function with amplitude A2siny and frequency 2f. Prod can be
generated
sample by sample, and (p may be calculated from the amplitude of the resulting
sine wave.
The calculation of phase is particularly dependent upon the accuracy of
previous calculations (i.e., the calculation of the frequencies and amplitudes
of SVI and
SV2). The controller may use multiple methods to provide separate (if not
entirely
independent) estimates of the phase, which may be combined to give an improved
estimate.
f. Drive Signal Generation
The controller generates the drive signal by applying a gain to the difference
between signals SVI and SV2. The controller may apply either a positive gain
(resulting
in positive feedback) or a negative gain (resulting in negative feedback).
In general, the Q of the conduit is high enough that the conduit will resonate
only at certain discrete frequencies. For example, the lowest resonant
frequency for some
conduits is between 65 Hz and 95 Hz, depending on the density of the process
fluid, and
irrespective of the drive frequency. As such, it is desirable to drive the
conduit at the
24
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
resonant frequency to minimize cycle-to-cycle energy loss. Feeding back the
sensor
voltage to the drivers permits the drive frequency to migrate to the resonant
frequency.
As an alternative to using feedback to generate the drive signal, pure sine
waves having phases and frequencies determined as described above may be
synthesized
and sent to the drivers. This approach offers the advantage of eliminating
undesirable
high frequency components, such as harmonics of the resonant frequency. This
approach
also permits compensation for time delays introduced by the analog-to-digital
converters,
processing, and digital-to-analog converters to ensure that the phase of the
drive signal
corresponds to the mid-point of the phases of the sensor signals. This
compensation may
be provided by determining the time delay of the system components and
introducing a
phase shift corresponding to the time delay.
Another approach to driving the conduit is to use square wave pulses. This is
another synthesis method, with fixed (positive and negative) direct current
sources being
switched on and off at timed intervals to provide the required energy. The
switching is
synchronized with the sensor voltage phase. Advantageously, this approach does
not
require digital-to-analog converters.
In general, the amplitude of vibration of the conduit should rapidly achieve a
desired value at startup, so as to quickly provide the measurement function,
but should do
so without significant overshoot, which may damage the meter. The desired
rapid startup
may be achieved by setting a very high gain so that the presence of random
noise and the
high Q of the conduit are sufficient to initiate motion of the conduit. In one
implementation, high gain and positive feedback are used to initiate motion of
the
conduit. Once stable operation is attained, the system switches to a synthesis
approach
for generating the drive signals.
Referring to Figs. 11A-13D, synthesis methods also may be used to initiate
conduit motion when high gain is unable to do so. For example, if the DC
voltage offset
of the sensor voltages is significantly larger than random noise, the
application of a high
gain will not induce oscillatory motion. This condition is shown in Figs. 11A-
11D, in
which a high gain is applied at approximately 0.3 seconds. As shown in Figs.
11A and
11B, application of the high gain causes one of the drive signals to assume a
large
positive value (Fig. 11A) and the other to assume a large negative value (Fig.
11B). The
magnitudes of the drive signals vary with noise in the sensor signals (Figs.
11C and 11D).
However, the amplified noise is insufficient to vary the sign of the drive
signals so as to
induce oscillation.
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
Figs. 12A-12D illustrate that imposition of a square wave over several cycles
can reliably cause a rapid startup of oscillation. Oscillation of a conduit
having a two
inch diameter may be established in approximately two seconds. The
establishment of
conduit oscillation is indicated by the reduction in the amplitude of the
drive signals, as
shown in Figs. 12A and 12B. Figs. 13A-13D illustrate that oscillation of a one
inch
conduit may be established in approximately half a second.
A square wave also may be used during operation to correct conduit
oscillation problems. For example, in some circumstances, flow meter conduits
have
been known to begin oscillating at harmonics of the resonant frequency of the
conduit,
) such as frequencies on the order of 1.5 kHz. When such high frequency
oscillations are
detected, a square wave having a more desirable frequency may be used to
return the
conduit oscillation to the resonant frequency.
g. Measurement Generation
5 The controller digitally generates the mass flow measurement in a
manner
similar to the approach used by the analog controller. The controller also may
generate
other measurements, such as density.
In one implementation, the controller calculates the mass flow based on the
phase difference in degrees between the two sensor signals (phase_diff), the
frequency of
) oscillation of the conduit (freq), and the process temperature (temp):
= temp - Tc,
noneu_mf = tan en * phase_diff/180), and
massflow = 16 (MF1*Tz2 + MF2*T, + MF3) * noneu_mf/freq,
where T, is a calibration temperature, MF i-MF3 are calibration constants
calculated
5 during a calibration procedure, and noneu_mf is the mass flow in non-
engineering units.
The controller calculates the density based on the frequency of oscillation of
the conduit and the process temperature:
= temp - Tc,
c2 = freq2, and
density = (Di*Tz2 + D2*Tz + D3)/c2 + a4*Tz2,
where D1-D4 are calibration constants generated during a calibration
procedure.
26
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
D. Integration Techniques
Many integration techniques are available, with different techniques requiring
different levels of computational effort and providing different levels of
accuracy. In the
described implementation, variants of Simpson's method are used. The basic
technique
may be expressed as:
rtn+2
ydt4Y Yn+2
3 n
where tk is the time at sample k, yk is the corresponding function value, and
h is the step
length. This rule can be applied repeatedly to any data vector with an odd
number of data
) points (i.e., three or more points), and is equivalent to fitting and
integrating a cubic
spline to the data points. If the number of data points happens to be even,
then the so-
called 3/8ths rule can be applied at one end of the interval:
r , 3h
yat ¨8 kyn 3y+1 +3Y n+2 Y n+3) =
5
As stated earlier, each cycle begins and ends at some offset (e.g.,
start_offset_SV1) from a sampling point. The accuracy of the integration
techniques are
improved considerably by taking these offsets into account. For example, in an
integration of a half cycle sine wave, the areas corresponding to partial
samples must be
) included in the calculations to avoid a consistent underestimate in the
result.
Two types of function are integrated in the described calculations: either
sine
or sine-squared functions. Both are easily approximated close to zero where
the end
points occur. At the end points, the sine wave is approximately linear and the
sine-
squared function is approximately quadratic.
5 In view of these two types of functions, three different
integration methods
have been evaluated. These are Simpson's method with no end correction,
Simpson's
method with linear end correction, and Simpson's method with quadratic
correction.
The integration methods were tested by generating and sampling pure sine and
sine-squared functions, without simulating any analog-to-digital truncation
error.
) Integrals were calculated and the results were compared to the true
amplitudes of the
signals. The only source of error in these calculations was due to the
integration
techniques. The results obtained are illustrated in tables A and B.
27
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
Table A - Integration of a sine function
% error (based on 1000 simulations) Av.Errors (%) S.D. Error (%)
Max.Error (%)
Simpson Only -3.73e-3 1.33e-3 6.17e-3
Simpson + linear correction 3.16e-8 4.89e-8 1.56e-7
Simpson + quadratic correction 2.00e-4 2.19e-2 5.18e-1
Table B - Integration of a sine-squared function
% error (based on 1000 simulations Av.Errors (%) S.D. Error (%)
Max.Error (%)
Simpson Only -2.21e-6 1.10e-6 4.39e-3
Simpson + linear correction 2.21e-6 6.93e-7 2.52e-6
Simpson + quadratic correction 2.15e-11 6.83e-11 1.88e-10
For sine functions, Simpson's method with linear correction was unbiased with
the smallest standard deviation, while Simpson's method without correction was
biased to
a negative error and Simpson's method with quadratic correction had a
relatively high
standard deviation. For sine-squared functions, the errors were generally
reduced, with
the quadratic correction providing the best result. Based on these
evaluations, linear
correction is used when integrating sine functions and quadratic correction is
used when
integrating sine-squared functions.
E. Synchronous Modulation Technique
Fig. 14 illustrates an alternative procedure 1400 for processing the sensor
signals. Procedure 1400 is based on synchronous modulation, such as is
described by
Denys et al., in "Measurement of Voltage Phase for the French Future Defence
Plan
Against Losses of Synchronism", IEEE Transactions on Power Delivery, 7(1), 62-
69,
1992 and by Begovic et al. in "Frequency Tracking in Power Networks in the
Presence of
Harmonics", IEEE Transactions on Power Delivery, 8(2), 480-486, 1993, both of
which
are incorporated by reference.
First, the controller generates an initial estimate of the nominal operating
frequency of the system (step 1405). The controller then attempts to measure
the
deviation of the frequency of a signal x[k] (e.g., SVi) from this nominal
frequency:
x[k] = A sin [(coo Au) kh + (130] + e (k),
28
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
where A is the amplitude of the sine wave portion of the signal, coo is the
nominal
frequency (e.g., 88 Hz), Aw is the deviation from the nominal frequency, h is
the
sampling interval, szto is the phase shift, and 8(k) corresponds to the added
noise and
harmonics.
To generate this measurement, the controller synthesizes two signals that
oscillate at the nominal frequency (step 1410). The signals are phase shifted
by 0 and
ir/2 and have amplitude of unity. The controller multiplies each of these
signals by the
original signal to produce signals yi and y2 (step 1415):
Yi = x[k]cos( okh) = ¨A sin [(20 + A co)kh + 01+ ¨A sin (Acokh + 0), and
2 2
A
Y2
where
2 2
where the first terms of yl and y2 are high frequency (e.g., 176 Hz)
components and the
second terms are low frequency (e.g., 0 Hz) components. The controller then
eliminates
the high frequency components using a low pass filter (step 1420):
5
/ A
y = ¨sin (Acokh + (I))+ 61[k], and
2
/ A
yj = ¨ cos (Acokh + (1)+ e 2[k],
2
where ci[k] and 62[k] represent the filtered noise from the original signals.
The
controller combines these signals to produce u[k] (step 1425):
u[k] = [k] + jy/2 [k])(y; [k ¨1] + jy/2[k ¨1])
= ui [k] + ju2[k]
A2
=-- ____________________ cos (Acoh) + j¨A2 sin (Acoh),
4 4
5 which carries the essential information about the frequency deviation. As
shown, ui [k]
represents the real component of u[k], while u2[k] represents the imaginary
component.
The controller uses the real and imaginary components of u[k] to calculate the
frequency deviation, AS (step 1430):
Af --, ¨1 arctan u2 [k]
al [k]
29
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
The controller then adds the frequency deviation to the nominal frequency
(step 1435) to
give the actual frequency:
f = Af + fo.
The controller also uses the real and imaginary components of u[k] to
determine the amplitude of the original signal. In particular, the controller
determines the
amplitude as (step 1440):
A2 = 4.1q[k]+4[k].
0
Next, the controller determines the phase difference between the two sensor
signals (step 1445). Assuming that any noise (61[k] and 62[k]) remaining after
application of the low pass filter described below will be negligible, noise
free versions of
y1 '[k] and y21[k] (yi*[k] and y2*[k]) may be expressed as:
* A .
5 y1[k] = ¨sm(Acokh + 0), and
2
y[k] = ¨A cos(Acokh ¨(1)).
2
Multiplying these signals together gives:
A 2
* * .
V = Y1 Y2 = km(20) + sin(26,wkh)1
8
Filtering this signal by a low pass filter having a cutoff frequency near 0 Hz
removes the
unwanted component and leaves:
, A2 .
v =
8
from which the phase difference can be calculated as:
1 .
= ¨arcsm _________________________________ .
2 A2
This procedure relies on the accuracy with which the operating frequency is
5 initially estimated, as the procedure measures only the deviation from
this frequency. If a
good estimate is given, a very narrow filter can be used, which makes the
procedure very
accurate. For typical flowmeters, the operating frequencies are around 95 Hz
(empty) and
82 Hz (full). A first approximation of half range (88 Hz) is used, which
allows a
low-pass filter cut-off of 13 Hz. Care must be taken in selecting the cut-off
frequency as
) a very small cut-off frequency can attenuate the amplitude of the sine
wave.
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
The accuracy of measurement also depends on the filtering characteristics
employed. The attenuation of the filter in the dead-band determines the amount
of
harmonics rejection, while a smaller cutoff frequency improves the noise
rejection.
F. Meter with PI Control
Figs. 15A and 15B illustrate a meter 1500 having a controller 1505 that uses
another technique to generate the signals supplied to the drivers. Analog-to-
digital
converters 1510 digitize signals from the sensors 48 and provide the digitized
signals to
the controller 1505. The controller 1505 uses the digitized signals to
calculate gains for
0 each driver, with the gains being suitable for generating desired
oscillations in the
conduit. The gains may be either positive or negative. The controller 1505
then supplies
the gains to multiplying digital-to-analog converters 1515. In other
implementations, two
or more multiplying digital-to-analog converters arranged in series may be
used to
implement a single, more sensitive multiplying digital-to-analog converter.
5 The controller 1505 also generates drive signals using the digitized
sensor
signals. The controller 1505 provides these drive signals to digital-to-analog
converters
1520 that convert the signals to analog signals that are supplied to the
multiplying digital-
to-analog converters 1515.
The multiplying digital-to-analog converters 1515 multiply the analog signals
0 by the gains from the controller 1505 to produce signals for driving the
conduit.
Amplifiers 1525 then amplify these signals and supply them to the drivers 46.
Similar
results could be obtained by having the controller 1505 perform the
multiplication
performed by the multiplying digital-to-analog converter, at which point the
multiplying
digital-to-analog converter could be replaced by a standard digital-to-analog
converter.
5 Fig. 15B illustrates the control approach in more detail. Within the
controller
1505, the digitized sensor signals are provided to an amplitude detector 1550,
which
determines a measure, a(t), of the amplitude of motion of the conduit using,
for example,
the technique described above. A summer 1555 then uses the amplitude a(t) and
a desired
amplitude ao to calculate an error e(t) as:
0 e (t) = ao ¨a (t).
The error e(t) is used by a proportional-integral ("PI") control block 1560 to
generate a
gain Ko(t). This gain is multiplied by the difference of the sensor signals to
generate the
drive signal. The PI control block permits high speed response to changing
conditions.
31
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
The amplitude detector 1550, summer 1555, and PI control block 1560 may be
implemented as software processed by the controller 1505, or as separate
circuitry.
1. Control Procedure
The meter 1500 operates according to the procedure 1600 illustrated in Fig.
16. Initially, the controller receives digitized data from the sensors (step
1605).
Thereafter, the procedure 1600 includes three parallel branches: a measurement
branch
1610, a drive signal generation branch 1615, and a gain generation branch
1620.
In the measurement branch 1610, the digitized sensor data is used to generate
0 measurements
of amplitude, frequency, and phase, as described above (step 1625). These
measurements then are used to calculate the mass flow rate (step 1630) and
other process
variables. In general, the controller 1505 implements the measurement branch
1610.
In the drive signal generation branch 1615, the digitized signals from the two
sensors are differenced to generate the signal (step 1635) that is multiplied
by the gain to
5 produce the
drive signal. As described above, this differencing operation is performed by
the controller 1505. In general, the differencing operation produces a
weighted difference
that accounts for amplitude differences between the sensor signals.
In the gain generation branch 1620, the gain is calculated using the
proportional-integral control block. As noted above, the amplitude, a(t), of
motion of the
0 conduit is determined (step 1640) and subtracted from the desired
amplitude ao (step
1645) to calculate the error e(t). Though illustrated as a separate step,
generation of the
amplitude, a(t), may correspond to generation of the amplitude in the
measurement
generation step 1625. Finally, the PI control block uses the error e(t) to
calculate the gain
(step 1650).
5 The
calculated gain is multiplied by the difference signal to generate the drive
signal supplied to the drivers (step 1655). As described above, this
multiplication
operation is performed by the multiplying D/A converter or may be performed by
the
controller.
2. PI Control Block
The objective of the PI control block is to sustain in the conduit pure
sinusoidal oscillations having an amplitude ao. The behavior of the conduit
may be
modeled as a simple mass-spring system that may be expressed as:
32
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
+ 2cconi+w,z2x = 0,
where x is a function of time and the displacement of the mass from
equilibrium, wr, is
the natural frequency, and C is a damping factor, which is assumed to be small
(e.g.,
0.001). The solution to this force equation as a function of an output y(t)
and an input i(t)
is analogous to an electrical network in which the transfer function between a
supplied
current, i(s), and a sensed output voltage, y(s), is:
y (s) ks
2
(S) s2 + 2cons + con
To achieve the desired oscillation in the conduit, a positive-feedback loop
having the gain Ko(t) is automatically adjusted by a 'slow' outer loop to
give:
0 .5e+ (go), ¨kx0(0)i+0)2õx = O.
The system is assumed to have a "two-time-scales" property, which means that
variations
in K0(t) are slow enough that solutions to the equation for x provided above
can be
obtained by assuming constant damping.
5 A two-term PI control block that gives zero steady-state error may
be
expressed as:
Ko(t)= K pe(t)+ K of e(t)dt,
where the error, e(t) (i.e., ao-a(0), is the input to the PI control block,
and Kp and K, are
constants. In one implementation, with ao = 10, controller constants of
Kp=0.02 and
K,=0.0005 provide a response in which oscillations build up quickly. However,
this PI
control block is nonlinear, which may result in design and operational
difficulties.
A linear model of the behavior of the oscillation amplitude may be derived by
assuming that x(t) equals AEJO, which results in:
= Aej't + jwei't, and
5
_ [A _ (02Alejo)t 2iwikeicot
Substituting these expressions into the expression for oscillation of the
loop, and
separating into real and imaginary terms, gives:
jw{2A.+ (gwn ¨ kKc,).A}= 0, and
A + (gcor, ¨1d(0)A.+(con2 ¨(o2)A = O.
A(t) also may be expressed as:
33
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
A
=
2
A solution of this equation is:
log A(t) = (¨ + id(cl
2
Transforming variables by defining a(t) as being equal to log A(t), the
equation for A(t)
can be written as:
da 1d( o (t)
¨dt = --Cwn
2
where 1(0 is now explicitly dependent on time. Taking Laplace transforms
results in:
0
a(s) = Co)n Ic-Ko (s)/ 2
which can be interpreted in terms of transfer-functions as in Fig. 17. This
figure is of
particular significance for the design of controllers, as it is linear for all
and a, with the
5 only assumption being the two-time-scales property. The performance of
the closed-loop
is robust with respect to this assumption so that fast responses which are
attainable in
practice can be readily designed.
From Fig. 17, the term Ccon, is a "load disturbance" that needs to be
eliminated
by the controller (i.e., kK.0/2 must equal ccon for a(t) to be constant). For
zero
D steady-state error this implies that the outer-loop controller must have
an integrator (or
very large gain). As such, an appropriate PI controller, C(s), may be assumed
to be
Kp(1+1/sT/), where Ti is a constant. The proportional term is required for
stability. The
term ccon, however, does not affect stability or controller design, which is
based instead
on the open-loop transfer function:
a(s) kK p (1 + sTi) kK p I 2(s +11 Tz)
5 C(s)G(s) = __
e(s) 2s21, s2
34
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
The root locus for varying Kp is shown in Fig. 18. For small Kp, there are
slow underdamped roots. As Kp increases, the roots become real at the point P
for which
the controller gain is Kp=8/(kT1). Note in particular that the theory does not
place any
restriction upon the choice of T,. Hence, the response can, in principle, be
made critically
damped and as fast as desired by appropriate choices of Kp and Ti.
Although the poles are purely real at point P, this does not mean there is no
overshoot in the closed-loop step response. This is most easily seen by
inspecting the
transfer function between the desired value, ao, and the error e:
e(s) s2 2
o(s) s2 +0.5kKp(s+11T,) p2 (s)
0 where p2 is a second-order polynomial. With a step input, a0(s)=u/s, the
response can be
written as apt(t), where p(t) is the inverse transform of 1/p2(s) and equals
al exp(-21/4,10+a2exp(-k2t). The signal p(t) increases and then decays to zero
so that e(t),
which is proportional to p', must change sign, implying overshoot in a(t). The
set-point ao
may be prefiltered to give a pseudo set-point ao*:
5 a; (s) = __ 1a0(s),
1+ s
where T, is the known controller parameter. With this prefilter, real
controller poles
should provide overshoot-free step responses. This feature is useful as there
may be
physical constraints on overshoot (e.g., mechanical interference or
overstressing of
components).
The root locus of Fig. 18 assumes that the only dynamics are from the
inner-loop's gain/log-amplitude transfer function (Fig. 16) and the outer-
loop's PI
controller C(s) (i.e., that the log-amplitude a = log A is measured
instantaneously).
However, A is the amplitude of an oscillation which might be growing or
decaying and
hence cannot in general be measured without taking into account the underlying
sinusoid.
5 There are several possible methods for measuring A, in addition to those
discussed above.
Some are more suitable for use in quasi-steady conditions. For example, a
phase-locked
loop in which a sinusoidal signal s(t) = sin(wõt+00) locks onto the measured
waveform
y(t) = A(Osin(o)nt+01:01) may be employed. Thus, a measure of the amplitude a
= log A is
given by dividing these signals (with appropriate safeguards and filters).
This method is
) perhaps satisfactory near the steady-state but not for start-up
conditions before there is a
lock.
CA 02658810 2009-02-25
WO 2008/027908 PCT/US2007/077018
Another approach uses a peak-follower that includes a zero-crossing detector
together with a peak-following algorithm implemented in the controller. Zero-
crossing
methods, however, can be susceptible to noise. In addition, results from a
peak-follower
are available only every half-cycle and thereby dictate the sample interval
for controller
updates.
Finally, an AM detector may be employed. Given a sine wave y(t) = Asincont,
an estimate of A may be obtained from A10.57iF {abs(y)}, where F{} is a
suitable
low-pass filter with unity DC gain. The AM detector is the simplest approach.
Moreover, it does not presume that there are oscillations of any particular
frequency, and
0 hence is usable during startup conditions. It suffers from a disadvantage
that there is a
leakage of harmonics into the inner loop which will affect the spectrum of the
resultant
oscillations. In addition, the filter adds extra dynamics into the outer loop
such that
compromises need to be made between speed of response and spectral purity. In
particular, an effect of the filter is to constrain the choice of the best Ti.
5 The Fourier series for abs(y) is known to be:
A abs(sin cunt)= ¨2A [2coiit4wõt 6wn t +
it .
3 15 35
As such, the output has to be scaled by 7E/2 to give the correct DC output A,
and the
(even) harmonic terms akcos2kwnt have to be filtered out. As all the filter
needs to do is
to pass the DC component through and reduce all other frequencies, a "brick-
wall" filter
with a cut-off below 2w, is sufficient. However, the dynamics of the filter
will affect the
behavior of the closed-loop. A common choice of filter is in the Butterworth
form. For
example, the third-order low-pass filter with a design break-point frequency
cob is:
1
F(s) =
102 i,s3
1+2S/Wb T_L / UJb -1- / UJb
At the design frequency the response is 3dB down; at 2ob it is -18dB (0.12),
and at 4(4 it
5 is -36dB (0.015) down. Higher-order Butterworth filters have a steeper
roll-off, but most
of their poles are complex and may affect negatively the control-loop's root
locus.
G. Zero Offset Compensation
As noted above, zero offset may be introduced into a sensor voltage signal by
) drift in the pre-amplification circuitry and by the analog-to-digital
converter. Slight
differences in the pre-amplification gains for positive and negative voltages
due to the use
36
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
of differential circuitry may worsen the zero offset effect. The errors vary
between
transmitters, and with transmitter temperature and component wear.
Audio quality (i.e., relatively low cost) analog-to-digital converters may be
employed for economic reasons. These devices are not designed with DC offset
and
amplitude stability as high priorities. Figs. 19A-19D show how offset and
positive and
negative gains vary with chip operating temperature for one such converter
(the AD1879
converter). The repeatability of the illustrated trends is poor, and even
allowing for
temperature compensation based on the trends, residual zero offset and
positive/negative
gain mismatch remain.
0 If phase is calculated using the time difference between zero
crossing points
on the two sensor voltages, DC offset may lead to phase errors. This effect is
illustrated
by Figs 20A-20C. Each graph shows the calculated phase offset as measured by
the
digital transmitter when the true phase offset is zero (i.e., at zero flow).
Fig. 20A shows phase calculated based on whole cycles starting with positive
5 zero-crossings. The mean value is 0.00627 degrees.
Fig. 20B shows phase calculated starting with negative zero-crossings. The
mean value is 0.0109 degrees.
Fig. 20C shows phase calculated every half-cycle. Fig. 20C interleaves the
data from Figs. 20A and 20B. The average phase (-0.00234) is closer to zero
than in
0 Figs. 20A and 20B, but the standard deviation of the signal is about six
times higher.
More sophisticated phase measurement techniques, such as those based on
Fourier methods, are immune to DC offset. However, it is desirable to
eliminate zero
offset even when those techniques are used, since data is processed in whole-
cycle
packets delineated by zero crossing points. This allows simpler analysis of
the effects of,
5 for example, amplitude modulation on apparent phase and frequency. In
addition, gain
mismatch between positive and negative voltages will introduce errors into any
measurement technique.
The zero-crossing technique of phase detection may be used to demonstrate
the impact of zero offset and gain mismatch error, and their consequent
removal. Figs.
0 21A and 21B illustrate the long term drift in phase with zero flow. Each
point represents
an average over one minute of live data. Fig. 21A shows the average phase, and
Fig. 21B
shows the standard deviation in phase. Over several hours, the drift is
significant. Thus,
even if the meter were zeroed every day, which in many applications would be
considered
an excessive maintenance requirement, there would still be considerable phase
drift.
37
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
1. Compensation Technique
A technique for dealing with voltage offset and gain mismatch uses the
computational capabilities of the digital transmitter and does not require a
zero flow
condition. The technique uses a set of calculations each cycle which, when
averaged over
a reasonable period (e.g., 10,000 cycles), and excluding regions of major
change (e.g., set
point change, onset of aeration), converge on the desired zero offset and gain
mismatch
compensations.
Assuming the presence of up to three higher harmonics, the desired waveform
0 for a sensor voltage SV(t) is of the form:
SV(t) = A1 sin(c)t) + A2 sin(2o)t) + A3 sin(30)t) + A4 sin(4cot)
where A1 designates the amplitude of the fundamental frequency component and
A2-A4
designate the amplitudes of the three harmonic components. However, in
practice, the
actual waveform is adulterated with zero offset Zo (which has a value close to
zero) and
5 mismatch between the negative and positive gains Gõ and G. Without any
loss of
generality, it can be assumed that Gp equals one and that Gn is given by:
Gn = 1 + 8G,
where eG represents the gain mismatch.
The technique assumes that the amplitudes A, and the frequency a) are
o constant. This is justified because estimates of Z0 and EG are based on
averages taken
over many cycles (e.g., 10,000 interleaved cycles occurring in about 1 minute
of
operation). When implementing the technique, the controller tests for the
presence of
significant changes in frequency and amplitude to ensure the validity of the
analysis. The
presence of the higher harmonics leads to the use of Fourier techniques for
extracting
5 phase and amplitude information for specific harmonics. This entails
integrating SV(t)
and multiplying by a modulating sine or cosine function.
The zero offset impacts the integral limits, as well as the functional form.
Because there is a zero offset, the starting point for calculation of
amplitude and phase
will not be at the zero phase point of the periodic waveform SV(t). For zero
offset Zo, the
o corresponding phase offset is, approximately,
r
. Zo
9z0 =
38
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
Zo Tzo
(Pzo =- ' tz =
For small phase, Ai with corresponding time delay
The integrals are scaled so that the limiting value (i.e., as Zo and o.)G
approach
zero) equals the amplitude of the relevant harmonic. The first two integrals
of interest
are:
It
--az
20) (1)
TIPS = ¨TE (SV(t) Z0) = sin [o)(t - t0 )] dt, and
tZ0
Tr
2- +tm
2o)
IlNs - ri ¨71 ki EG) (SV(t) Z0) = sin ko(t - t0 )1dt.
TC
¨tzo
CO
These integrals represent what in practice is calculated during a normal
Fourier analysis
of the sensor voltage data. The subscript 1 indicates the first harmonic, N
and P indicate,
0 respectively, the negative or positive half cycle, and s and c indicate,
respectively,
whether a sine or a cosine modulating function has been used.
Strictly speaking, the mid-zero crossing point, and hence the corresponding
integral limits, should be given by 7t/co - tzo, rather than Tc/o) + tzo.
However, the use of
the exact mid-point rather than the exact zero crossing point leads to an
easier analysis,
5 and better numerical behavior (due principally to errors in the location
of the zero
crossing point). The only error introduced by using the exact mid-point is
that a small
section of each of the above integrals is multiplied by the wrong gain (1
instead of 1 + 6G
and vice versa). However, these errors are of order Zo2cG and are considered
negligible.
Using computer algebra and assuming small Zo and CG, first order estimates
o for the integrals may be derived as:
IlPs est
=A1 + ¨4 Z0 1 + 2A2 4A4
______________________________________________ , and
it 3A1 15A1
42A2 4A4
IlNs est =((l+ CG) Al L01-3 ¨Ai + ¨15 ¨Ai =
Useful related functions including the sum, difference, and ratio of the
integrals and their estimates may be determined. The sum of the integrals may
be
expressed as:
5 SUMis (
Jips + 1Ns))
39
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
while the sum of the estimates equals:
4 7 , 2 A2 4 A4
&Mats est =
3 Ai 15 Ai
Similarly, the difference of the integrals may be expressed as:
='lPs ¨/,Nõ
while the difference of the estimates is:
, 2 A2 4 A4
Dtisest = A16GL'Ok+EG) --+ =
3 Ai 15 Ai
Finally, the ratio of the integrals is:
D,.
Ratiois = _______________________________ " '
while the ratio of the estimates is:
0
1 8 15A1 +10A2 + 4./14
Ratiois õt = ________________________ + 4 __________
1
1+ EG 15 RA?
Corresponding cosine integrals are defined as:
= 201 ¨co+rzo (SV(t)+ Zo)cos[co(t ¨t0 )]dt, and
It tzo
5
27E
261(1 + 6G ) 7,0) (SV(t) + Zo) cos[
¨ w(t ¨ tz0 )idt,
it --ktZo
with estimates:
40A2 +16A4 , and
¨ ¨ + _______
/iPc est ZO
157E
40A2 +16A4]
IlNc est = (1+ EG) Z 0
15n
and sums:
Sumic = pc Il NC and
+16A4
Sum1C est = GI 4 A2
Z0 + =
57C
Second harmonic integrals are:
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
12Ps 2wj (SM+ Zi3) sin [26)(t ¨ tz0)] dt, and
tzo
27t
¨ tz0
= ¨2to
Co
I2Ns (1+ EG).1" (SVW + Zo sin [2w(t - tzo )1 dt,
Co
¨Ptzo
with estimates:
A
/ = A + 5 + 9 1 and
2Ps _ est 2 157c 0 ,
8 A
I2Ps _est = 0+8 G ) A2 ________________ Z0 ¨5 +
157C A11'
and sums:
Sum2sI
= -2Ps I2Ns and
8 A3
SUM2p5 est = A2(2 + )15nsGZ0 ¨5+9¨ =
Ai
The integrals can be calculated numerically every cycle. As discussed below,
the equations estimating the values of the integrals in terms of various
amplitudes and the
zero offset and gain values are rearranged to give estimates of the zero
offset and gain
terms based on the calculated integrals.
0
2. Example
The accuracy of the estimation equations may be illustrated with an example.
For each basic integral, three values are provided: the "true" value of the
integral
(calculated within Mathcad using Romberg integration), the value using the
estimation
5 equation, and the value calculated by the digital transmitter
operating in simulation mode,
using Simpson's method with end correction.
Thus, for example, the value for lips calculated according to:
7C
20)
I 1 ps =(SV(t) Z 0 ) sin[o)(t ¨ tz0 )]dt
tzo
is 0.101353, while the estimated value Ps_est) calculated as:
4 , -2A2 4A4
IlPs est = +Z0r 1 + A1 +
7 17-5 --A-1
is 0.101358. The value calculated using the digital transmitter in simulation
mode is
0.101340. These calculations use the parameter values illustrated in Table C.
41
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
Table C
Parameter Value Comment
co 1607c This corresponds to frequency = 80Hz, a typical
value.
A1 0.1 This more typically is 0.3, but it could be smaller
with aeration.
A2 0.01 This more typically is 0.005, but it could be larger
with aeration.
A3 and A4 0.0 The digital Coriolis simulation mode only offers two
harmonics, so these
higher harmonics are ignored. However, they are small (<0.002).
Zo 0.001 Experience
suggests that this is a large value for zero offset.
80 0.001 Experience
suggests this is a large value for gain mismatch.
The exact, estimate and simulation results from using these parameter values
are
illustrated in Table D.
Table D
Integral 'Exact' Value Estimate Digital Coriolis
Simulation
'fps 0.101353 0.101358 0.101340
0.098735 0.098740 0.098751
Iipc 0.007487 0.007488 0.007500
11Nc -0.009496 -0.009498 -0.009531
121's 0.009149 0.009151 0.009118
I2Ns 0.010857 0.010859 0.010885
0 Thus, at least for the particular values selected, the estimates
given by the first
order equations are extremely accurate. As Zo and EG approach zero, the errors
in both
the estimate and the simulation approach zero.
3. Implementation
5 The first order estimates for the integrals define a series of non-
linear
equations in terms of the amplitudes of the harmonics, the zero offset, and
the gain
mismatch. As the equations are non-linear, an exact solution is not readily
available.
However, an approximation followed by corrective iterations provides
reasonable
convergence with limited computational overhead.
42
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
Conduit-specific ratios may be assumed for A1-A4. As such, no attempt is
made to calculate all of the amplitudes A1-A4. Instead, only A1 and A2 are
estimated
using the integral equations defined above. Based on experience of the
relative
amplitudes, A3 may be approximated as A2/2, and A4 may be approximated as
A2/10.
The zero offset compensation technique may be implemented according to the
procedure 2200 illustrated in Fig. 22. During each cycle, the controller
calculates the
integrals I1Ps, IlNs, IINc,
12Põ I2Ns and related functions sumis, ratioiõ sumic and sum25
(step 2205). This requires minimal additional calculation beyond the
conventional
Fourier calculations used to determine frequency, amplitude and phase.
o Every 10,000 cycles, the controller checks on the slope of the sensor
voltage
amplitude A1, using a conventional rate-of-change estimation technique (step
2210). If
the amplitude is constant (step 2215), then the controller proceeds with
calculations for
zero offset and gain mismatch. This check may be extended to test for
frequency
stability.
5 To perform the calculations, the controller generates average values
for the
functions (e.g., sumis) over the last 10,000 cycles. The controller then makes
a first
estimation of zero offset and gain mismatch (step 2225):
Zo = -Sumic/2, and
80 = 1/Ratio15 - 1
0 Using these values, the controller calculates an inverse gain factor
(k) and
amplitude factor (amp_factor) (step 2230):
k = 1.0 / (1.0 + 0.5 * eG), and
amp_factor = 1 + 50/75 * Sum2s/Sumis
The controller uses the inverse gain factor and amplitude factor to make a
first estimation
5 of the amplitudes (step 2235):
A1 = k * [Sum15/2 + 2/7c * Zo * co * ampfactor], and
A2 = k * [Sum25/2 - 4/(3 * n) * Zo * CO
43
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
The controller then improves the estimate by the following calculations,
iterating
as required (step 2240):
xi ¨ Zo,
X2= EG,
EG = [1 + 8/n * xi/Ai * amp_factor]/Ratiois - 1.0,
Zo = - Sumic/2 + x2 * (xi + 2.773/rc*A2)/2,
Ai = k * [Sumis/2 + 2/n * xi * x2 * amp factor],
A2 = k * [Sum2s/2 - 4/(15 * n) * x1 * x2 * (5 - 4.5 * A2)].
The controller uses standard techniques to test for convergence of the values
of Zo and
0 EG. In practice the corrections are small after the first iteration, and
experience suggests
that three iterations are adequate.
Finally, the controller adjusts the raw data to eliminate Zo and EG (step
2245).
The controller then repeats the procedure. Once zero offset and gain mismatch
have been
eliminated from the raw data, the functions (i.e., sumis) used in generating
subsequent
5 values for Zo and 6G. are based on corrected data. Accordingly, these
subsequent values
for Zo and 6G reflect residual zero offset and gain mismatch, and are summed
with
previously generated values to produce the actual zero offset and gain
mismatch. In one
approach to adjusting the raw data, the controller generates adjustment
parameters (e.g.,
S l_off and S2_off) that are used in converting the analog signals from the
sensors to
0 digital data.
Figs. 23A-23C, 24A and 24B show results obtained using the procedure 2200.
The short-term behavior is illustrated in Figs. 23A-23C. This shows
consecutive phase
estimates obtained five minutes after startup to allow time for the procedure
to begin
affecting the output. Phase is shown based on positive zero-crossings,
negative zero-
5 crossings, and both.
The difference between the positive and negative mean values has been reduced
by a factor of 20, with a corresponding reduction in mean zero offset in the
interleaved
data set. The corresponding standard deviation has been reduced by a factor of
approximately 6.
Longer term behavior is shown in Figs. 24A and 24B. The initial large zero
offset
is rapidly corrected, and then the phase offset is kept close to zero over
many hours. The
average phase offset, excluding the first few values, is 6.14e-6, which
strongly suggests
44
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
that the procedure is successful in compensating for changes in voltage offset
and gain
imbalance.
Typical values for Zo and EG for the digital Coriolis meter are Zo = -7.923e-4
and
EG = -1.754e-5 for signal SVI, and Zo = -8.038e-4 and E G = +6.93e-4 for
signal SV2.
H. Dynamic Analysis
In general, conventional measurement calculations for Coriolis meters assume
that
the frequency and amplitude of oscillation on each side of the conduit are
constant, and
that the frequency on each side of the conduit is identical and equal to the
so-called
o resonant frequency. Phases generally are not measured separately for each
side of the
conduit, and the phase difference between the two sides is assumed to be
constant for the
duration of the measurement process. Precise measurements of frequency, phase
and
amplitude every half-cycle using the digital meter demonstrate that these
assumptions are
only valid when parameter values are averaged over a time period on the order
of
5 seconds. Viewed at 100 Hz or higher frequencies, these parameters exhibit
considerable
variation. For example, during normal operation, the frequency and amplitude
values of
SVI may exhibit strong negative correlation with the corresponding SV2 values.
Accordingly, conventional measurement algorithms are subject to noise
attributable to
these dynamic variations. The noise becomes more significant as the
measurement
0 calculation rate increases. Other noise terms may be introduced by
physical factors, such
as flowtube dynamics, dynamic non-linearities (e.g. flowtube stiffness varying
with
amplitude), or the dynamic consequences of the sensor voltages providing
velocity data
rather than absolute position data.
The described techniques exploit the high precision of the digital meter to
monitor
5 and compensate for dynamic conduit behavior to reduce noise so as to
provide more
precise measurements of process variables such as mass flow and density. This
is
achieved by monitoring and compensating for such effects as the rates of
change of
frequency, phase and amplitude, flowtube dynamics, and dynamic physical non-
idealities.
A phase difference calculation which does not assume the same frequency on
each side
D has already been described above. Other compensation techniques are
described below.
Monitoring and compensation for dynamic effects may take place at the
individual
sensor level to provide corrected estimates of phase, frequency, amplitude or
other
parameters. Further compensation may also take place at the conduit level,
where data
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
from both sensors are combined, for example in the calculation of phase
difference and
average frequency. These two levels may be used together to provide
comprehensive
compensation.
Thus, instantaneous mass flow and density measurements by the flowmeter may
be improved by modelling and accounting for dynamic effects of flowmeter
operation. In
general, 80% or more of phase noise in a Coriolis flowmeter may be attributed
to
flowtube dynamics (sometimes referred to as "ringing"), rather than to process
conditions
being measured. The application of a dynamic model can reduce phase noise by a
factor
of 4 to 10, leading to significantly improved flow measurement performance. A
single
o model is effective for all flow rates and amplitudes of oscillation.
Generally,
computational requirements are negligible.
The dynamic analysis may be performed on each of the sensor signals in
isolation
from the other. This avoids, or at least delays, modelling the dynamic
interaction
between the two sides of the conduit, which is likely to be far more complex
than the
5 dynamics at each sensor. Also, analyzing the individual sensor signals is
more likely to
be successful in circumstances such as batch startup and aeration where the
two sides of
the conduit are subject to different forces from the process fluid.
In general, the dynamic analysis considers the impact of time-varying
amplitude,
frequency and phase on the calculated values for these parameters. While the
frequency
0 and amplitude are easily defined for the individual sensor voltages,
phase is
conventionally defined in terms of the difference between the sensor voltages.
However,
when a Fourier analysis is used, phase for the individual sensor may be
defined in terms
of the difference between the midpoint of the cycle and the average 180 phase
point.
Three types of dynamic effects are measurement error and the so-called
5 "feedback" and "velocity" effects. Measurement error results because the
algorithms for
calculating amplitude and phase assume that frequency, amplitude, and phase
are constant
over the time interval of interest. Performance of the measurement algorithms
may be
improved by correcting for variations in these parameters.
The feedback effect results from supplying energy to the conduit to make up
for
energy loss from the conduit so as to maintain a constant amplitude of
oscillation. The
need to add energy to the conduit is only recognized after the amplitude of
oscillation
begins to deviate from a desired setpoint. As a result, the damping term in
the equation of
motion for the oscillating conduit is not zero, and, instead, constantly
dithers around zero.
46
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
Although the natural frequency of the conduit does not change, it is obscured
by shifts in
the zero-crossings (i.e., phase variations) associated with these small
changes in
amplitude.
The velocity effect results because the sensor voltages observe conduit
velocity,
but are analyzed as being representative of conduit position. A consequence of
this is that
the rate of change of amplitude has an impact on the apparent frequency and
phase, even
if the true values of these parameters are constant.
1. Sensor-Level Compensation for Amplitude Modulation
One approach to correcting for dynamic effects monitors the amplitudes of the
sensor signals and makes adjustments based on variations in the amplitudes.
For
purposes of analyzing dynamic effects, it is assumed that estimates of phase,
frequency
and amplitude may be determined for each sensor voltage during each cycle. As
shown
in Fig. 25, calculations are based on complete but overlapping cycles. Each
cycle starts at
5 a zero crossing point, halfway through the previous cycle. Positive
cycles begin with
positive voltages immediately after the initial zero-crossing, while negative
cycles begin
with negative voltages. Thus cycle n is positive, while cycles n-1 and n+1 are
negative.
It is assumed that zero offset correction has been performed so that zero
offset is
negligible. It also is assumed that higher harmonics may be present.
Linear variation in amplitude, frequency, and phase are assumed. Under this
assumption, the average value of each parameter during a cycle equals the
instantaneous
value of the parameter at the mid-point of the cycle. Since the cycles overlap
by 180
degrees, the average value for a cycle equals the starting value for the next
cycle.
For example, cycle n is from time 0 to 27c/co. The average values of
amplitude,
frequency and phase equal the instantaneous values at the mid-point, rc/oD,
which is also
the starting point for cycle n+1, which is from time ir/o) to 37r/o). Of
course, these
timings are approximate, since co also varies with time.
a. Dynamic Effect Compensation Procedure
The controller accounts for dynamic effects according to the procedure 2600
illustrated in Fig. 26. First, the controller produces a frequency estimate
(step 2605) by
using the zero crossings to measure the time between the start and end of the
cycle, as
47
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
described above. Assuming that frequency varies linearly, this estimate equals
the time-
averaged frequency over the period.
The controller then uses the estimated frequency to generate a first estimate
of
amplitude and phase using the Fourier method described above (step 2610). As
noted
above, this method eliminates the effects of higher harmonics.
Phase is interpreted in the context of a single waveform as being the
difference
between the start of the cycle (i.e., the zero-crossing point) and the point
of zero phase for
the component of SV(t) of frequency co, expressed as a phase offset. Since the
phase
offset is an average over the entire waveform, it may be used as the phase
offset from the
midpoint of the cycle. Ideally, with no zero offset and constant amplitude of
oscillation,
the phase offset should be zero every cycle. In practice, however, it shows a
high level of
variation and provides an excellent basis for correcting mass flow to account
for dynamic
changes in amplitude.
The controller then calculates a phase difference (step 2615). Though a number
of
5 definitions of phase difference are possible, the analysis assumes that
the average phase
and frequency of each sensor signal is representative of the entire waveform.
Since these
frequencies are different for SVI and SV2, the corresponding phases are scaled
to the
average frequency. In addition, the phases are shifted to the same starting
point (i.e., the
midpoint of the cycle on SV1). After scaling, they are subtracted to provide
the phase
) difference.
The controller next determines the rate of change of the amplitude for the
cycle n
(step 2620):
roc amp
amp (end of cycle) - amp (start of cycle)
_ .
period of cycle
= (ampn_o ¨ ) freq..
This calculation assumes that the amplitude from cycle n+1 is available when
calculating
5 the rate of change of cycle n. This is possible if the corrections are
made one cycle after
the raw amplitude calculations have been made. The advantage of having an
accurate
estimate of the rate of change, and hence good measurement correction,
outweighs the
delay in the provision of the corrected measurements, which, in one
implementation, is on
the order of 5 milliseconds. The most recently generated information is always
used for
) control of the conduit (i.e., for generation of the drive signal).
48
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
If desired, a revised estimate of the rate of change can be calculated after
amplitude correction has been applied (as described below). This results in
iteration to
convergence for the best values of amplitude and rate of change.
b. Frequency Compensation for Feedback and Velocity Effects
As noted above, the dynamic aspects of the feedback loop introduce time
varying
shifts in the phase due to the small deviations in amplitude about the set-
point. This
results in the measured frequency, which is based on zero-crossings, differing
from the
natural frequency of the conduit. If velocity sensors are used, an additional
shift in phase
) occurs. This additional shift is also associated with changes in the
positional amplitude of
the conduit. A dynamic analysis can monitor and compensate for these effects.
Accordingly, the controller uses the calculated rate of amplitude change to
correct the
frequency estimate (step 2625).
The position of an oscillating conduit in a feedback loop that is employed to
i maintain the amplitude of oscillation of the conduit constant may be
expressed as:
X = A(t) sin(o)ot ¨ 0(0,
where e(t) is the phase delay caused by the feedback effect. The mechanical Q
of the
oscillating conduit is typically on the order of 1000, which implies small
deviations in
amplitude and phase. Under these conditions, 0(0 is given by:
A(t)
0(t)
2w0 A(t)
Since each sensor measures velocity:
SV(t) = k(t) = A(t) sin[w t - 0(t)] + [coo ¨ OW]A(t) cos[coo t ¨0(t)]
\2-1/2
( N / =
(t)
= co0A(t) 1¨ 2 A
¨ + ________ cos(co 0 t ¨ 0(t) ¨
coo woA(t)i
where y(t) is the phase delay caused by the velocity effect:
49
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
y(t) = tan-1 =
0
woAN 1---
'SSince the mechanical Q of the conduit is typically on the order of 1000,
and, hence,
variations in amplitude and phase are small, it is reasonable to assume:
A(t)
-6 and <<1.
0 ( 0 0 A(t)
This means that the expression for SV(t) may be simplified to:
SV(t) coo A(t)cos(o) ot - 0(t) ¨7(t)),
and for the same reasons, the expression for the velocity offset phase delay
may be
simplified to:
)
y(t) At
co 0 A(t)
Summing the feedback and velocity effect phase delays gives the total phase
delay:
A(t) A(t) A(t)
9(0 = (t) + y (t)
ao 0 A(t) w 0 A(t) 2co0A(t)'
and the following expression for SV(t):
SV(t) woA(t)cos[000t ¨ 9(t)].
From this, the actual frequency of oscillation may be distinguished from the
natural frequency of oscillation. Though the former is observed, the latter is
useful for
density calculations. Over any reasonable length of time, and assuming
adequate
amplitude control, the averages of these two frequencies are the same (because
the
average rate of change of amplitude must be zero). However, for improved
instantaneous
density measurement, it is desirable to compensate the actual frequency of
oscillation for
dynamic effects to obtain the natural frequency. This is particularly useful
in dealing
with aerated fluids for which the instantaneous density can vary rapidly with
time.
The apparent frequency observed for cycle n is delineated by zero crossings
occurring at the midpoints of cycles n-1 and n+1. The phase delay due to
velocity change
will have an impact on the apparent start and end of the cycle:
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
obs _ freq. = obs_freq._, + true_ freqn
(qn+i ¨ 9n-1)
2n
= obs_ freq.-1 + true_ freq._, n+1 An-1
271 47t true
freq.A1 zin true_ freq.A., )
r = = \
1 An+1 An-1
= obs freq. +
87c2 An+1 A71 =
Based on this analysis, a correction can be applied using an integrated error
term:
r -
error¨ sumn error sum ¨ ¨1 An+1 An-1 , and
8n2 \, An+1 An-1
est_ freq. = obs _ freq. ¨ error sum. ,
where the value of error sum at startup (i.e., the value at cycle zero) is:
r =
1 A0 A1
error_ sump = ___________________________ + ¨ .
87E2 A0 A1)
Though these equations include a constant term having a value of 1/87E2 ,
actual data has
indicated that a constant term of 1 /87E is more appropriate. This discrepancy
may be due
to unmodelled dynamics that may be resolved through further analysis.
The calculations discussed above assume that the true amplitude of
oscillation, A,
) is available. However, in practice, only the sensor voltage SV is
observed. This sensor
voltage may be expressed as:
SV(t) co 0 A 0 cos(coot - 90)
The amplitude, amp_SV(t), of this expression is:
amp_ SV(t) co0;40 .
The rate of change of this amplitude is:
roc _ amp _SV(t) w0;40
so that the following estimation can be used:
AN roc ¨amp _ SV(t)
AN amp SV(t)
c. Application of Feedback and Velocity Effect Frequency Compensation
Figs. 27A-32B illustrate how application of the procedure 2600 improves the
estimate of the natural frequency, and hence the process density, for real
data from a
meter having a one inch diameter conduit. Each of the figures shows 10,000
samples,
which are collected in just over 1 minute.
51
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
Figs. 27A and 27B show amplitude and frequency data from SVi, taken when
random changes to the amplitude set-point have been applied. Since the conduit
is full of
water and there is no flow, the natural frequency is constant. However, the
observed
frequency varies considerably in response to changes in amplitude. The mean
frequency
value is 81.41 Hz, with a standard deviation of 0.057 Hz.
Figs. 28A and 28B show, respectively, the variation in frequency from the mean
value, and the correction term generated using the procedure 2600. The gross
deviations
are extremely well matched. However, there is additional variance in frequency
which is
not attributable to amplitude variation. Another important feature illustrated
by Fig. 28B
is that the average is close to zero as a result of the proper initialization
of the error term,
as described above.
Figs. 29A and 92B compare the raw frequency data (Fig. 29A) with the results
of
applying the correction function (Fig. 29B). There has been a negligible shift
in the mean
frequency, while the standard deviation has been reduced by a factor of 4.4.
From Fig.
29B, it is apparent that there is residual structure in the corrected
frequency data. It is
expected that further analysis, based on the change in phase across a cycle
and its impact
on the observed frequency, will yield further noise reductions.
Figs. 30A and 30B show the corresponding effect on the average frequency,
which is the mean of the instantaneous sensor voltage frequencies. Since the
mean
frequency is used to calculate the density of the process fluid, the noise
reduction (here by
a factor of 5.2) will be propagated into the calculation of density.
Figs. 31A and 31B illustrate the raw and corrected average frequency for a 2"
diameter conduit subject to a random amplitude set-point. The 2" flowtube
exhibits less
frequency variation that the 1", for both raw and corrected data. The noise
reduction
factor is 4Ø
Figs. 32A and 32B show more typical results with real flow data for the one
inch
flowtube. The random setpoint algorithm has been replaced by the normal
constant
setpoint. As a result, there is less amplitude variation than in the previous
examples,
which results in a smaller noise reduction factor of 1.5.
d. Compensation of Phase Measurement for Amplitude Modulation
Referring again to Fig. 26, the controller next compensates the phase
measurement to account for amplitude modulation assuming the phase calculation
provided above (step 2630). The Fourier calculations of phase described above
assume
52
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
that the amplitude of oscillation is constant throughout the cycle of data on
which the
calculations take place. This section describes a correction which assumes a
linear
variation in amplitude over the cycle of data.
Ignoring higher harmonics, and assuming that any zero offset has been
eliminated,
the expression for the sensor voltage is given by:
SV(t) A1(1+XAt)sin(cot)
where XA is a constant corresponding to the relative change in amplitude with
time. As
discussed above, the integrals II and 12 may be expressed as:
I, = f¨ w SV(t)sin(wt)dt, and
n 0
2o)
'2= ¨ f w SV(t)cos(o)t)dt.
It
Evaluating these integrals results in:
11=A1 1+ ¨ A ,and
/2 = -20) AA.
Substituting these expressions into the calculation for amplitude and
expanding as a series
in XA results in:
1
Amp = A1(1+Lc 2k,A 2
2 =
\ 80)
Assuming X,A is small, and ignoring all terms after the first order term, this
may
be simplified to:
Amp = (1 +X 11 A\ .
j
This equals the amplitude of SV(t) at the midpoint of the cycle (t = rc/ co).
Accordingly,
the amplitude calculation provides the required result without correction.
For the phase calculation, it is assumed that the true phase difference and
frequency are constant, and that there is no voltage offset, which means that
the phase
value should be zero. However, as a result of amplitude modulation, the
correction to be
applied to the raw phase data to compensate for amplitude modulation is:
X
Phase = tan-1 A
2(7a,A +0)))
Assuming that the expression in brackets is small, the inverse tangent
function can be
ignored.
53
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
A more elaborate analysis considers the effects of higher harmonics. Assuming
that the sensor voltage may be expressed as:
SV(t) = (1+ k 40[Ai sin(wt) + A2 sin(20)t) + A3 sin(wt) + A4 sin(46)01
such that all harmonic amplitudes increase at the same relative rate over the
cycle, then
the resulting integrals may be expressed as:
/, = Ai (1 + A ,
CO /
and
¨1
2
________________________ k,k304 +80A2 + 45A3 +32A4)
600)
-
for positive cycles, and
/, --1X -(30A, ¨80A2 +454 ¨3244)
60o)
for negative cycles.
For amplitude, substituting these expressions into the calculations
establishes that
the amplitude calculation is only affected in the second order and higher
terms, so that no
correction is necessary to a first order approximation of the amplitude. For
phase, the
correction term becomes:
¨1, (30A1 +80A2 +45A3 +32A4
60 A _________________
A ( tk A CO)
for positive cycles, and
¨1, (30A1 ¨80A2 +45A3 ¨32A4
0 A __________________
Ai (rcX26, + (0)
for negative cycles. These correction terms assume the availability of the
amplitudes of
the higher harmonics. While these can be calculated using the usual Fourier
technique, it
is also possible to approximate some or all them using assumed ratios between
the
harmonics. For example, for one implementation of a one inch diameter conduit,
typical
amplitude ratios are A1 = 1.0, A2 = 0.01, A3 = 0.005, and A4 = 0.001.
e. Application of Amplitude Modulation Compensation to Phase
Simulations have been carried out using the digital transmitter, including the
simulation of higher harmonics and amplitude modulation. One example uses f =
80 Hz,
A1(t=0) = 0.3, A2 = 0, A3 = 0, A4 = 0, 2 = 1e-5 * 48KHz (sampling rate) ¨
0.47622,
54
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
which corresponds to a high rate of change of amplitude, but with no higher
harmonics.
Theory suggests a phase offset of -0.02706 degrees. In simulation over 1000
cycles, the
average offset is -0.02714 degrees, with a standard deviation of only 2.17e-6.
The
difference between simulation and theory (approx 0.3% of the simulation error)
is
attributable to the model's assumption of a linear variation in amplitude over
each cycle,
while the simulation generates an exponential change in amplitude.
A second example includes a second harmonic, and has the parameters f = 80Hz,
Ai (t=0) = 0.3, A2(t=0) = 0.003, A3 = 0, A4 = 0, XA = -1e-6 * 48KHz (sampling
rate) =
-0.047622. For this example, theory predicts the phase offset to be +2.706e-3,
+/-2.66%
for positive or negative cycles. In simulation, the results are 2.714e-3 +/ -
2.66%, which
again matches well.
Figs. 33A-34B give examples of how this correction improves real flowmeter
data. Fig. 33A shows raw phase data from SVi, collected from a 1" diameter
conduit,
with low flow assumed to be reasonably constant. Fig. 33B shows the correction
factor
calculated using the formula described above, while Fig. 33C shows the
resulting
corrected phase. The most apparent feature is that the correction has
increased the
variance of the phase signal, while still producing an overall reduction in
the phase
difference (i.e., SV2 - SVI) standard deviation by a factor of 1.26, as shown
in Figs. 34A
and 34B. The improved performance results because this correction improves the
correlation between the two phases, leading to reduced variation in the phase
difference.
The technique works equally well in other flow conditions and on other conduit
sizes.
f. Compensation to Phase Measurement for Velocity Effect
The phase measurement calculation is also affected by the velocity effect. A
highly effective and simple correction factor, in radians, is of the form
1
cvk ) = -7tASV(t k )
where ASV(tk) is the relative rate of change of amplitude and may be expressed
as:
SV(t k, )¨ SV(tk-1 ) 1
ASV
(k)- =
tk - tk SV(t k )'
where tk is the completion time for the cycle for which ASV(tk) is being
determined, tk+1
is the completion time for the next cycle, and tk_i is the completion time of
the previous
cycle. ASV is an estimate of the rate of change of SV, scaled by its absolute
value, and
is also referred to as the proportional rate of change of SV.
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
Figs. 35A-35E illustrate this technique. Fig. 35A shows the raw phase data
from a
single sensor (SV1), after having applied the amplitude modulation corrections
described
above. Fig. 35B shows the correction factor in degrees calculated using the
equation
above, while Fig. 35C shows the resulting corrected phase. It should be noted
that the
standard deviation of the corrected phase has actually increased relative to
the raw data.
However, when the corresponding calculations take place on the other sensor
(SV2), there
is an increase in the negative correlation (from -0.8 to -0.9) between the
phases on the two
signals. As a consequence, the phase difference calculations based on the raw
phase
measurements (Fig. 35D) have significantly more noise than the corrected phase
measurements (Fig. 35E).
Comparison of Figs. 35D and 35E shows the benefits of this noise reduction
technique. It is immediately apparent from visual inspection of Fig. 35E that
the process
variable is decreasing, and that there are significant cycles in the
measurement, with the
cycles being attributable, perhaps, to a poorly conditioned pump. None of this
is
discernable from the uncorrected phase difference data of Fig. 35D.
g. Application of Sensor Level Noise Reduction
The combination of phase noise reduction techniques described above results in
substantial improvements in instantaneous phase difference measurement in a
variety of
flow conditions, as illustrated in Figs. 36A-36L. Each graph shows three phase
difference
measurements calculated simultaneously in real time by the digital Coriolis
transmitter
operating on a one inch conduit. The middle band 3600 shows phase data
calculated
using the simple time-difference technique. The outermost band 3605 shows
phase data
calculated using the Fourier-based technique described above.
It is perhaps surprising that the Fourier technique, which uses far more data,
a
more sophisticated analysis, and much more computational effort, results in a
noisier
calculation. This can be attributed to the sensitivity of the Fourier
technique to the
dynamic effects described above. The innermost band of data 3610 shows the
same
Fourier data after the application of the sensor-level noise reduction
techniques. As can
be seen, substantial noise reduction occurs in each case, as indicated by the
standard
deviation values presented on each graph.
Fig. 36A illustrates measurements with no flow, a full conduit, and no pump
noise. Fig. 36B illustrates measurements with no flow, a full conduit, and the
pumps on.
Fig. 36C illustrates measurements with an empty, wet conduit. Fig. 36D
illustrates
56
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
measurements at a low flow rate. Fig. 36E illustrates measurements at a high
flow rate.
Fig. 36F illustrates measurements at a high flow rate and an amplitude of
oscillation of
0.03V. Fig. 36G illustrates measurements at a low flow rate with low aeration.
Fig. 36H
illustrates measurements at a low flow rate with high aeration. Fig. 361
illustrates
measurements at a high flow rate with low aeration. Fig. 36J illustrates
measurements at
a high flow rate with high aeration. Fig. 36K illustrates measurements for an
empty to
high flow rate transition. Fig. 36L illustrates measurements for a high flow
rate to empty
transition.
2. Flowtube Level Dynamic Modelling
A dynamic model may be incorporated in two basic stages. In the first stage,
the
model is created using the techniques of system identification. The flowtube
is
"stimulated" to manifest its dynamics, while the true mass flow and density
values are
kept constant. The response of the flowtube is measured and used in generating
the
dynamic model. In the second stage, the model is applied to normal flow data.
Predictions of the effects of flowtube dynamics are made for both phase and
frequency.
The predictions then are subtracted from the observed data to leave the
residual phase and
frequency, which should be due to the process alone. Each stage is described
in more
detail below.
a. System Identification
System identification begins with a flowtube full of water, with no flow. The
amplitude of oscillation, which normally is kept constant, is allowed to vary
by assigning
a random setpoint between 0.05 V and 0.3 V, where 0.3 V is the usual value.
The
resulting sensor voltages are shown in Fig. 37A, while Figs. 37B and 37C show,
respectively, the corresponding calculated phase and frequency values. These
values are
calculated once per cycle. Both phase and frequency show a high degree of
"structure."
Since the phase and frequency corresponding to mass flow are constant, this
structure is
likely to be related to flowtube dynamics. Observable variables that will
predict this
structure when the true phase and frequency are not known to be constant may
be
expressed as set forth below.
First, as noted above, ASV(tk) may be expressed as:
57
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
ASV(tk ) = SV(tk+i)- Sntk-1 1
t k+i - t k_i AST(t k ).
This expression may be used to determine ASVI and ASV2.
The phase of the flowtube is related to A, which is defined as ASVi - ASV2,
while the frequency is related to At, which is defined as ASVi + ASV2. These
parameters are illustrated in Figs. 37D and 37E. Comparing Fig. 37B to Fig.
37D and
Fig. 37C to Fig. 37E shows the striking relationship between A- and phase and
between
A* and frequency.
Some correction for flowtube dynamics may be obtained by subtracting a
multiple
of the appropriate prediction function from the phase and/or the frequency.
Improved
results may be obtained using a model of the form:
y(k) + aiy(k-1) + + any(k-n) = bou(k) + biu(k-1) + + bn,u(k-m),
where y(k) is the output (i.e., phase or frequency) and u is the prediction
function (i.e., A-
or At). The technique of system identification suggests values for the orders
n and m,
and the coefficients a, and b, of what are in effect polynomials in time. The
value of y(k)
can be calculated every cycle and subtracted from the observed phase or
frequency to get
the residual process value.
It is important to appreciate that, even in the absence of dynamic
corrections, the
digital flowmeter offers very good precision over a long period of time. For
example,
when totalizing a batch of 200kg, the device readily achieves a repeatability
of less that
0.03%. The purpose of the dynamic modelling is to improve the dynamic
precision.
Thus, raw and compensated values should have similar mean values, but
reductions in
"variance" or "standard deviation."
Figs. 38A and 39A show raw and corrected frequency values. The mean values
are similar, but the standard deviation has been reduced by a factor of 3.25.
Though the
gross deviations in frequency have been eliminated, significant "structure"
remains in the
residual noise. This structure appears to be unrelated to the A+ function. The
model
used is a simple first order model, where m = n = 1.
Figs. 38B and 39B show the corresponding phase correction. The mean value is
minimally affected, while the standard deviation is reduced by a factor of
7.9. The model
58
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
orders are n = 2 and m = 10. Some structure appears to remain in the residual
noise. It is
expected that this structure is due to insufficient excitation of the phase
dynamics by
setpoint changes.
More effective phase identification has been achieved through further
simulation
of flowtube dynamics by continuous striking of the flowtube during data
collection (set
point changes are still carried out). Figs. 38C and 39C show the effects of
correction
under these conditions. As shown, the standard deviation is reduced by a
factor of 31.
This more effective model is used in the following discussions.
b. Application to Flow Data
The real test of an identified model is the improvements it provides for new
data.
At the outset, it is useful to note a number of observations. First, the mean
phase,
averaged over, for example, ten seconds or more, is already quite precise. In
the
examples shown, phase values are plotted at 82 Hz or thereabouts. The reported
standard
deviation would be roughly 1/3 of the values shown when averaged to 10Hz, and
1/9
when averages to 1 Hz. For reference, on a one inch flow tube, one degree of
phase
difference corresponds to about 1 kg/s flow rate.
The expected benefit of the technique is that of providing a much better
dynamic
response to true process changes, rather than improving average precision.
Consequently,
in the following examples, where the flow is non-zero, small flow step changes
are
introduced every ten seconds or so, with the expectation that the corrected
phase will
show up the step changes more clearly.
Figs. 38D and 39D show the correction applied to a full flowtube with zero
flow,
just after startup. The ring-down effect characteristic of startup is clearly
evident in the
raw data (Fig. 38D), but this is eliminated by the correction (Fig. 39D),
leading to a
standard deviation reduction of a factor of 23 over the whole data set. Note
that the
corrected measurement closely resembles white noise, suggesting most flowtube
dynamics have been captured.
Figs. 38E and 39E show the resulting correction for a "drained" flowtube.
Noise
is reduced by a factor of 6.5 or so. Note, however, that there appears to be
some residual
structure in the noise.
The effects of the technique on low (Figs. 38F and 39F), medium (Figs. 38G and
390), and high (Figs. 38H and 3911) flow rates are also illustrated, each with
step changes
in flow every ten seconds. In each case, the pattern is the same: the
corrected average
59
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
flows (Figs. 39F-39H) are identical to the raw average flows (Figs. 38F-38H),
but the
dynamic noise is reduced considerably. In Fig. 39H, this leads to the
emergence of the
step changes which previously had been submerged in noise (Fig. 38H).
3. Extensions of Dynamic Monitoring and Compensation Techniques
The previous sections have described a variety of techniques (physical
modelling,
system identification, heuristics) used to monitor and compensate for
different aspects of
dynamic behavior (frequency and phase noise caused by amplitude modulation,
velocity
effect, flowtube dynamics at both the sensor and the flowtube level). By
natural
extension, similar techniques well-known to practitioners of control and/or
instrumentation, including those of artificial intelligence, neural networks,
fuzzy logic,
and genetic algorithms, as well as classical modelling and identification
methods, may be
applied to these and other aspects of the dynamic performance of the meter.
Specifically,
these might include monitoring and compensation for frequency, amplitude
and/or phase
variation at the sensor level, as well as average frequency and phase
difference at the
flowtube level, as these variations occur within each measurement interval, as
well as the
time between measurement intervals (where measurement intervals do not
overlap).
This technique is unusual in providing both reduced noise and improved dynamic
response to process measurement changes. As such, the technique promises to be
highly
valuable in the context of flow measurement.
I. Aeration (Two-Phase Flow)
The digital flowmeter provides improved performance in the presence of
aeration
(also known as two-phase flow) in the conduit. Aeration causes energy losses
in the
conduit that can have a substantial negative impact on the measurements
produced by a
mass flowmeter and can result in stalling of the conduit. Experiments have
shown that
the digital flowmeter has substantially improved performance in the presence
of aeration
relative to traditional, analog flowmeters. This performance improvement may
stem from
the meter's ability to provide a very wide gain range, to employ negative
feedback, to
calculate measurements precisely at very low amplitude levels, and to
compensate for
dynamic effects such as rate of change of amplitude and flowtube dynamics. The
performance improvement also may stem from the meter's use of a precise
digital
amplitude control algorithm.
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
The digital flowmeter detects the onset of aeration when the required driver
gain
rises simultaneously with a drop in apparent fluid density. The digital
flowmeter then
may directly respond to detected aeration. In general, the meter monitors the
presence of
aeration by comparing the observed density of the material flowing through the
conduit
(i.e., the density measurement obtained through normal measurement techniques)
to the
known, non-aerated density of the material. The controller determines the
level of
aeration based on any difference between the observed and actual densities.
The
controller then corrects the mass flow measurement accordingly.
The controller determines the non-aerated density of the material by
monitoring
the density over time periods in which aeration is not present (i.e., periods
in which the
density has a stable value). Alternatively, a control system to which the
controller is
connected may provide the non-aerated density as an initialization parameter.
In one implementation, the controller uses three corrections to account for
the
effects of aeration: bubble effect correction, damping effect correction, and
sensor
imbalance correction. Figs. 40A-40H illustrate the effects of the correction
procedure.
Fig. 40A illustrates the error in the phase measurement as the measured
density
decreases (i.e., as aeration increases) for different mass flow rates, absent
aeration
correction. As shown, the phase error is negative and has a magnitude that
increases with
increasing aeration. Fig. 40B illustrates that the resulting mass flow error
also is
negative. It also is significant to note that the digital flowmeter operates
at high levels of
aeration. By contrast, as indicated by the vertical bar 4000, traditional
analog meters tend
to stall in the presence of low levels of aeration.
A stall occurs when the flowmeter is unable to provide a sufficiently large
driver
gain to allow high drive current at low amplitudes of oscillation. If the
level of damping
requires a higher driver gain than can be delivered by the flowtube in order
to maintain
oscillation at a certain amplitude, then insufficient drive energy is supplied
to the conduit.
This results in a drop in amplitude of oscillation, which in turn leads to
even less drive
energy supplied due to the maximum gain limit. Catastrophic collapse results,
and
flowtube oscillation is not possible until the damping reduces to a level at
which the
corresponding driver gain requirement can be supplied by the flowmeter.
The bubble effect correction is based on the assumption that the mass flow
decreases as the level of aeration, also referred to as the void fraction,
increases. Without
attempting to predict the actual relationship between void fraction and the
bubble effect,
this correction assumes, with good theoretical justification, that the effect
on the observed
61
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
mass flow will be the same as the effect on the observed density. Since the
true fluid
density is known, the bubble effect correction corrects the mass flow rate by
the same
proportion. This correction is a linear adjustment that is the same for all
flow rates. Figs.
40C and 40D illustrate, respectively, the residual phase and mass flow errors
after
correction for the bubble effect. As shown, the residual errors are now
positive and
substantially smaller in magnitude than the original errors.
The damping factor correction accounts for damping of the conduit motion due
to
aeration. In general, the damping factor correction is based on the following
relationship
between the observed phase, (Pobs, and the actual phase, (throe:
( 1,.1
2õ
`1" true
(Pobs = CI true -1 j_
f 2
where 21/4, is a damping coefficient and k is a constant. Fig. 40E illustrates
the damping
correction for different mass flow rates and different levels of aeration.
Fig. 40F
illustrates the residual phase error after correction for damping. As shown,
the phase
error is reduced substantially relative to the phase error remaining after
bubble effect
correction.
The sensor balance correction is based on density differences between
different
ends of the conduit. As shown in Fig. 41, a pressure drop between the inlet
and the outlet
of the conduit results in increasing bubble sizes from the inlet to the
outlet. Since
material flows serially through the two loops of the conduit, the bubbles at
the inlet side
of the conduit (i.e., the side adjacent to the first sensor/driver pair) will
be smaller than
the bubbles at the outlet side of the conduit (i.e., the side adjacent to the
second
sensor/driver pair). This difference in bubble size results in a difference in
mass and
density between the two ends of the conduit. This difference is reflected in
the sensor
signals (5V1 and SV2). Accordingly, the sensor balance correction is based on
a ratio of
the two sensor signals.
Fig. 40G illustrates the sensor balance correction for different mass flow
rates and
different levels of aeration. Fig. 40H illustrates the residual phase error
after applying the
sensor balance correction. At low flow rates and low levels of aeration, the
phase error is
improved relative to the phase error remaining after damping correction.
Other correction factors also may be used. For example, the phase angle of
each
sensor signal may be monitored. In general, the average phase angle for a
signal should
be zero. However, the average phase angle tends to increase with increasing
aeration.
62
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
Accordingly, a correction factor could be generated based on the value of the
average
phase angle. Another correction factor could be based on the temperature of
the conduit.
In general, application of the correction factors tends to keep the mass flow
errors
at one percent or less. Moreover, these correction factors appear to be
applicable over a
wide range of flows and aeration levels.
J. Setpoint Adjustment
The digital flowmeter provides improved control of the setpoint for the
amplitude
of oscillation of the conduit. In an analog meter, feedback control is used to
maintain the
amplitude of oscillation of the conduit at a fixed level corresponding to a
desired peak
sensor voltage (e.g., 0.3 V). A stable amplitude of oscillation leads to
reduced variance in
the frequency and phase measurements.
In general, a large amplitude of oscillation is desirable, since such a large
amplitude provides a large Coriolis signal for measurement purposes. A large
amplitude
of oscillation also results in storage of a higher level of energy in the
conduit, which
provides greater robustness to external vibrations.
Circumstances may arise in which it is not possible to maintain the large
amplitude of oscillation due to limitations in the current that can be
supplied to the
drivers. For example, in one implementation of an analog transmitter, the
current is
limited to 100 mA for safety purposes. This is typically 5-10 times the
current needed to
maintain the desired amplitude of oscillation. However, if the process fluid
provides
significant additional damping (e.g., via two-phase flow), then the optimal
amplitude may
no longer be sustainable.
Similarly, a low-power flowmeter, such as the two-wire meter described below,
may have much less power available to drive the conduit. In addition, the
power level
may vary when the conduit is driven by capacitive discharge.
Referring to Fig. 42, a control procedure 4200 implemented by the controller
of
the digital flowmeter may be used to select the highest sustainable setpoint
given a
maximum available current level. In general, the procedure is performed each
time that a
desired drive current output is selected, which typically is once per cycle,
or once every
half-cycle if interleaved cycles are used.
The controller starts by setting the setpoint to a default value (e.g., 0.3 V)
and
initializing filtered representations of the sensor voltage (filtered SV) and
the drive
current (filtered_DC) (step 4205). Each time that the procedure is performed,
the
63
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
controller updates the filtered values based on current values for the sensor
voltage (SV)
and drive current (DC) (step 4210). For example, the controller may generate a
new
value for filtered_SV as the sum of ninety nine percent of filtered_SV and one
percent of
SV.
Next, the controller determines whether the procedure has been paused to
provide
time for prior setpoint adjustments to take effect (step 4215). Pausing of the
procedure is
indicated by a pause cycle count having a value greater than zero. If the
procedure is
paused, the controller performs no further actions for the cycle and
decrements the pause
cycle count (step 4220).
If the procedure has not been paused, the controller determines whether the
filtered drive current exceeds a threshold level (step 4225). In one
implementation, the
threshold level is ninety five percent of the maximum available current. If
the current
exceeds the threshold, the controller reduces the setpoint (step 4230). To
allow time for
the meter to settle after the setpoint change, the controller then implements
a pause of the
procedure by setting the pause cycle count equal to an appropriate value
(e.g., 100) (step
4235).
If the procedure has not been paused, the controller determines whether the
filtered drive current is less than a threshold level (step 4240) and the
setpoint is less than
a maximum permitted setpoint (step 4245). In one implementation, the threshold
level
equals seventy percent of the maximum available current. If both conditions
are met, the
controller determines a possible new setpoint (step 4250). In one
implementation, the
controller determines the new setpoint as eighty percent of the maximum
available
current multiplied by the ratio of filtered_SV to filtered_DC. To avoid small
changes in
the setpoint (i.e., chattering), the controller then determines whether the
possible new
setpoint exceeds the current setpoint by a sufficient amount (step 4255). In
one
implementation, the possible new setpoint must exceed the current setpoint by
0.02 V and
by ten percent.
If the possible new setpoint is sufficiently large, the controller determines
if it is
greater than the maximum permitted setpoint (step 4260). If so, the controller
sets the
setpoint equal to the maximum permitted setpoint (step 4265). Otherwise, the
controller
sets the setpoint equal to the possible new setpoint (step 4270). The
controller then
implements a pause of the procedure by setting the pause cycle count equal to
an
appropriate value (step 4235).
64
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
Figs. 43A-43C illustrate operation of the setpoint adjustment procedure. As
shown in Fig. 43C, the system starts with a setpoint of 0.3 V. At about eight
seconds of
operation, aeration results in a drop in the apparent density of the material
in the conduit
(Fig. 43A). Increased damping accompanying the aeration results in an increase
in the
drive current (Fig. 43B) and increased noise in the sensor voltage (Fig. 43C).
No changes
are made at this time, since the meter is able to maintain the desired
setpoint.
At about fifteen seconds of operation, aeration increases and the apparent
density
decreases further (Fig. 43A). At this level of aeration, the driver current
(Fig. 43B)
reaches a maximum value that is insufficient to maintain the 0.3 V setpoint.
Accordingly,
the sensor voltage drops to 0.26 V (Fig. 43C), the voltage level that the
maximum driver
current is able to maintain. In response to this condition, the controller
adjusts the
setpoint (at about 28 seconds of operation) to a level (0.23 V) that does not
require
generation of the maximum driver current
At about 38 seconds of operation, the level of aeration decreases and the
apparent
density increases (Fig. 43A). This results in a decrease in the drive current
(Fig. 43B).
At 40 seconds of operation, the controller responds to this condition by
increasing the
setpoint (Fig. 43C). The level of aeration decreases and the apparent density
increases
again at about 48 seconds of operation, and the controller responds by
increasing the
setpoint to 0.3 V.
K. Performance Results
The digital flowmeter has shown remarkable performance improvements relative
to traditional analog flowmeters. In one experiment, the ability of the two
types of meters
to accurately measure a batch of material was examined. In each case, the
batch was fed
through the appropriate flowmeter and into a tank, where the batch was
weighed. For
1200 and 2400 pound batches, the analog meter provided an average offset of
500
pounds, with a repeatability of 200 pounds. By contrast, the digital meter
provided an
average offset of 40 pounds, with a repeatability of two pounds, which clearly
is a
substantial improvement.
In each case, the conduit and surrounding pipework were empty at the start of
the
batch. This is important in many batching applications where it is not
practical to start
the batch with the conduit full. The batches were finished with the flowtube
full. Some
positive offset is expected because the flovvmeter is measuring the material
needed to fill
the pipe before the weighing tank starts to be filled. Delays in starting up,
or offsets
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
caused by aerated flow or low amplitudes of oscillation, are likely to
introduce negative
offsets. For real batching applications, the most important issue is the
repeatability of the
measurement.
The results show that with the analog flowmeter there are large negative
offsets
and repeatability of only 200 pounds. This is attributable to the length of
time taken to
startup after the onset of flow (during which no flow is metered), and
measurement errors
until full amplitude of oscillation is achieved. By comparison, the digital
flowmeter
achieves a positive offset, which is attributable to filling up of the empty
pipe, and a
repeatability of two pounds.
Another experiment compared the general measurement accuracy of the two types
of meters. Fig. 44 illustrates the accuracy and corresponding uncertainty of
the
measurements produced by the two types of meters at different percentages of
the meters'
maximum recommended flow rate. At high flow rates (i.e., at rates of 25% or
more of the
maximum rate), the analog meter produces measurements that correspond to
actual values
to within 0.15% or less, compared to 0.005% or less for the digital meter. At
lower rates,
the offset of the analog meter is on the order of 1.5%, compared to 0.25% for
the digital
meter.
L. Self-Validating Meter
The digital flowmeter may used in a control system that includes self-
validating
sensors. To this end, the digital flowmeter may be implemented as a self-
validating
meter. Self-validating meters and other sensors are described in U.S. Patent
No.
5,570,300, titled "SELF-VALIDATING SENSORS", which is incorporated by
reference.
In general, a self-validating meter provides, based on all information
available to
the meter, a best estimate of the value of a parameter (e.g., mass flow) being
monitored.
Because the best estimate is based, in part, on nonmeasurement data, the best
estimate
does not always conform to the value indicated by the current, possibly
faulty,
measurement data. A self-validating meter also provides information about the
uncertainty and reliability of the best estimate, as well as information about
the
operational status of the sensor. Uncertainty information is derived from
known
uncertainty analyses and is provided even in the absence of faults.
In general, a self-validating meter provides four basic parameters: a
validated
measurement value (VMV), a validated uncertainty (VU), an indication (MV
status) of
the status under which the measurement was generated, and a device status. The
VMV is
66
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
the meter's best estimate of the value of a measured parameter. The VU and the
MV
status are associated with the VMV. The meter produces a separate VMV, VU and
MV
status for each measurement. The device status indicates the operational
status of the
meter.
The meter also may provide other information. For example, upon a request from
a control system, the meter may provide detailed diagnostic information about
the status
of the meter. Also, when a measurement has exceeded, or is about to exceed, a
predetermined limit, the meter can send an alarm signal to the control system.
Different
alarm levels can be used to indicate the severity with which the measurement
has deviated
from the predetermined value.
VMV and VU are numeric values. For example, VMV could be a temperature
measurement valued at 200 degrees and VU, the uncertainty of VMV, could be 9
degrees.
In this case, there is a high probability (typically 95%) that the actual
temperature being
measured falls within an envelope around VMV and designated by VU (i.e., from
191
degrees to 209 degrees).
The controller generates VMV based on underlying data from the sensors. First,
the controller derives a raw measurement value (RMV) that is based on the
signals from
the sensors. In general, when the controller detects no abnormalities, the
controller has
nominal confidence in the RMV and sets the VMV equal to the RMV. When the
controller detects an abnormality in the sensor, the controller does not set
the VMV equal
to the RMV. Instead, the controller sets the VMV equal to a value that the
controller
considers to be a better estimate than the RMV of the actual parameter.
The controller generates the VU based on a raw uncertainty signal (RU) that is
the
result of a dynamic uncertainty analysis of the RMV. The controller performs
this
uncertainty analysis during each sampling period. Uncertainty analysis,
originally
described in "Describing Uncertainties in Single Sample Experiments," S.J.
Kline & F.A.
McClintock, Mech. Eng., 75, 3-8 (1953), has been widely applied and has
achieved the
status of an international standard for calibration. Essentially, an
uncertainty analysis
provides an indication of the "quality" of a measurement. Every measurement
has an
associated error, which, of course, is unknown. However, a reasonable limit on
that error
can often be expressed by a single uncertainty number (ANSI/ASME PTC 19.1-1985
Part
1, Measurement Uncertainty: Instruments and Apparatus).
As described by Kline & McClintock, for any observed measurement M, the
uncertainty in M, wm, can be defined as follows:
67
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
Mr. c {M ¨ Wm ,M wm
where M is true (114,e) with a certain level of confidence (typically 95%).
This
uncertainty is readily expressed in a relative form as a proportion of the
measurement (i.e.
Wm/M).
In general, the VU has a non-zero value even under ideal conditions (i.e., a
faultless sensor operating in a controlled, laboratory environment). This is
because the
measurement produced by a sensor is never completely certain and there is
always some
potential for error. As with the VMV, when the controller detects no
abnormalities, the
controller sets the VU equal to the RU. When the controller detects a fault
that only
partially affects the reliability of the RMV, the controller typically
performs a new
uncertainty analysis that accounts for effects of the fault and sets the VU
equal to the
results of this analysis. The controller sets the VU to a value based on past
performance
when the controller determines that the RMV bears no relation to the actual
measured
value.
To ensure that the control system uses the VMV and the VU properly, the MV
status provides information about how they were calculated. The controller
produces the
VMV and the VU under all conditions--even when the sensors are inoperative.
The
control system needs to know whether VMV and VU are based on "live" or
historical
data. For example, if the control system were using VMV and VU in feedback
control
and the sensors were inoperative, the control system would need to know that
VMV and
VU were based on past performance.
The MV status is based on the expected persistence of any abnormal condition
and on the confidence of the controller in the RMV. The four primary states
for MV
status are generated according to Table 1.
Table 1
Expected Confidence MV Status
Persistence in R1V1V
not applicable nominal CLEAR
not applicable reduced BLURRED
Short zero DAZZLED
Long zero BLIND
68
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
A CLEAR MV status occurs when RMV is within a normal range for given
process conditions. A DAZZLED MV status indicates that RMV is quite abnormal,
but
the abnormality is expected to be of short duration. Typically, the controller
sets the MV
status to DAZZLED when there is a sudden change in the signal from one of the
sensors
and the controller is unable to clearly establish whether this change is due
to an as yet
undiagnosed sensor fault or to an abrupt change in the variable being
measured. A
BLURRED MV status indicates that the RMV is abnormal but reasonably related to
the
parameter being measured. For example, the controller may set the MV status to
BLURRED when the RMV is a noisy signal. A BLIND MV status indicates that the
RMV is completely unreliable and that the fault is expected to persist.
Two additional states for the MV status are UNVALIDATED and SECURE. The
MV status is TIN VALIDATED when the controller is not performing validation of
VMV.
MV status is SECURE when VMV is generated from redundant measurements in which
the controller has nominal confidence.
The device status is a generic, discrete value summarizing the health of the
meter.
It is used primarily by fault detection and maintenance routines of the
control system.
Typically, the device status 32 is in one of six states, each of which
indicates a different
operational status for the meter. These states are: GOOD, TESTING, SUSPECT,
IMPAIRED, BAD, or CRITICAL. A GOOD device status means that the meter is in
nominal condition. A TESTING device status means that the meter is performing
a self
check, and that this self check may be responsible for any temporary reduction
in
measurement quality. A SUSPECT device status means that the meter has produced
an
abnormal response, but the controller has no detailed fault diagnosis. An
IMPAIRED
device status means that the meter is suffering from a diagnosed fault that
has a minor
impact on performance. A BAD device status means that the meter has seriously
malfunctioned and maintenance is required. Finally, a CRITICAL device status
means
that the meter has malfunctioned to the extent that the meter may cause (or
have caused) a
hazard such as a leak, fire, or explosion.
Fig. 45 illustrates a procedure 4500 by which the controller of a self-
validating
meter processes digitized sensor signals to generate drive signals and a
validated mass
flow measurement with an accompanying uncertainty and measurement status.
Initially,
the controller collects data from the sensors (step 4505). Using this data,
the controller
determines the frequency of the sensor signals (step 4510). If the frequency
falls within
an expected range (step 4515), the controller eliminates zero offset from the
sensor
69
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
signals (step 4520), and determines the amplitude (step 4525) and phase (step
4530) of
the sensor signals. The controller uses these calculated values to generate
the drive signal
(step 4535) and to generate a raw mass flow measurement and other measurements
(step
4540).
If the frequency does not fall within an expected range (step 4515), then the
controller implements a stall procedure (step 4545) to determine whether the
conduit has
stalled and to respond accordingly. In the stall procedure, the controller
maximizes the
driver gain and performs a broader search for zero crossings to determine
whether the
conduit is oscillating at all.
If the conduit is not oscillating correctly (i.e., if it is not oscillating,
or if it is
oscillating at an unacceptably high frequency (e.g., at a high harmonic of the
resonant
frequency)) (step 4550), the controller attempts to restart normal oscillation
(step 4555) of
the conduit by, for example, injecting a square wave at the drivers. After
attempting to
restart oscillation, the controller sets the MV status to DAZZLED (step 4560)
and
generates null raw measurement values (step 4565). If the conduit is
oscillating correctly
(step 4550), the controller eliminates zero offset (step 4520) and proceeds as
discussed
above.
After generating raw measurement values (steps 4540 or 4565), the controller
performs diagnostics (step 4570) to determine whether the meter is operating
correctly
(step 4575). (Note that the controller does not necessarily perform these
diagnostics
during every cycle.)
Next, the controller performs an uncertainty analysis (step 4580) to generate
a raw
uncertainty value. Using the raw measurements, the results of the diagnostics,
and other
information, the controller generates the VMV, the VU, the MV status, and the
device
status (step 4585). Thereafter, the controller collects a new set of data and
repeats the
procedure. The steps of the procedure 4500 may be performed serially or in
parallel, and
may be performed in varying order.
In another example, when aeration is detected, the mass flow corrections are
applied as described above, the MV status becomes blurred, and the uncertainty
is
increased to reflect the probable error of the correction technique. For
example, for a
flovvtube operating at 50% flowrate, under normal operating conditions, the
uncertainty
might be of the order of 0.1 - 0.2% of flowrate. If aeration occurs and is
corrected for
using the techniques described above, the uncertainty might be increased to
perhaps 2%
of reading. Uncertainty values should decrease as understanding of the effects
of aeration
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
improves and the ability to compensate for aeration gets better. In batch
situations, where
flow rate uncertainty is variable (e.g. high at start/end if batching from/to
empty, or
during temporary incidents of aeration or cavitation), the uncertainty of the
batch total
will reflect the weighted significance of the periods of high uncertainty
against the rest of
the batch with nominal low uncertainty. This is a highly useful quality metric
in fiscal
and other metering applications.
M. Two Wire Flowmeter
As shown in Fig. 46, the techniques described above may be used to implement a
"two-wire" Coriolis flowmeter 4600 that performs bidirectional communications
on a pair
of wires 4605. A power circuit 4610 receives power for operating a digital
controller
4615 and for powering the driver(s) 4620 to vibrate the conduit 4625. For
example, the
power circuit may include a constant output circuit 4630 that provides
operating power to
the controller and a drive capacitor 4635 that is charged using excess power.
The power
circuit may receive power from the wires 4605 or from a second pair of wires.
The
digital controller receives signals from one or more sensors 4640.
When the drive capacitor is suitably charged, the controller 4615 discharges
the
capacitor 4635 to drive the conduit 4625. For example, the controller may
drive the
conduit once during every 10 cycles. The controller 4615 receives and analyzes
signals
from the sensors 4640 to produce a mass flow measurement that the controller
then
transmits on the wires 4605.
N. Batching From Empty
The digital mass flowmeter 100 provides improved performance in dealing with a
challenging application condition that is referred to as batching from empty.
There are
many processes, particularly in the food and petrochemical industries, where
the high
accuracy and direct mass-flow measurement provided by Coriolis technology is
beneficial
in the metering of batches of material. In many cases, however, ensuring that
the
flowmeter remains full of fluid from the start to the end of the batch is not
practical, and
is highly inefficient. For example, in filling or emptying a tanker, air
entrainment is
difficult to avoid. In food processing, hygiene regulations may require pipes
to be
cleaned out between batches.
71
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
In conventional Coriolis meters, batching from empty may result in large
errors.
For example, hydraulic shock and a high gain requirement may be caused by the
onset of
flow in an empty flowtube, leading to large measurement errors and stalling.
The digital mass flowmeter 100 is robust to the conditions experienced when
batching from empty. More specifically, the amplitude controller has a rapid
response;
the high gain range prevents flowtube stalling; measurement data can be
calculated down
to 0.1 % of the normal amplitude of oscillation; and there is compensation for
the rate of
change of amplitude.
These characteristics are illustrated in Figs. 47A-47C, which show the
response of
the digital mass flowmeter 100 driving a wet and empty 25 mm flowtube during
the first
seconds of the onset of full flow. As shown in Fig. 47A, the drive gain
required to drive
the wet and empty tube prior to the onset of flow (at about 4.0 seconds) has a
value of
approximately 0.1, which is higher than the value of approximately 0.034
needed for a
full flowtube. The onset of flow is characterized by a substantial increase in
gain and a
corresponding drop in amplitude of oscillation. With reference to Fig. 47B, at
about 1.0
second after initiation, the selection of a reduced setpoint assists in the
stabilization of
amplitude while the full flow regime is established. After about 2.75 seconds,
the last of
the entrained air is purged, the conventional setpoint is restored, and the
drive gain
assumes the nominal value of 0.034. The raw and corrected phase difference
behavior is
shown in Fig. 47C.
As shown in Figs. 47A-47C, phase data is given continuously throughout the
transition. In similar circumstances, the analog control system stalls, and is
unable to
provide measurement data until the required drive gain returns to a near-
nominal value
and the lengthy start-up procedure is completed. As also shown, the correction
for the
rate of change of amplitude is clearly beneficial, particularly after 1.0
seconds.
Oscillations in amplitude result in substantial swings in the Fourier-based
and time-based
calculations of phase, but these are substantially reduced in the corrected
phase
measurement. Even in the most difficult part of the transition, from 0.4-1.0
seconds, the
correction provides some noise reduction.
Of course there are still erroneous data in this interval. For example, flow
generating a phase difference in excess of about 5 degrees is physically not
possible.
However, from the perspective of a self validating sensor, such as is
discussed above, this
phase measurement still constitutes raw data that may be corrected. In some
implementations, a higher level validation process may identify the data from
0.4-1.0
72
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
seconds as unrepresentative of the true process value (based on the gain,
amplitude and
other internal parameters), and may generate a DAZZLED mass-flow to suppress
extreme
measurement values.
Referring to Fig. 48A, the response of the digital mass flowmeter 100 to the
onset
of flow results in improved accuracy and repeatability. An exemplary flow rig
4800 is
shown in Fig. 48B. In producing the results shown in Fig. 48A, as in producing
the
results shown in Fig. 44, fluid was pumped through a magnetic flowmeter 4810
and a
Coriolis flowmeter 4820 into a weighing tank 4830, with the Coriolis flowmeter
being
either a digital flowmeter or a traditional analog flowmeter. Valves 4840 and
4860 were
used to ensure that the magnetic flowmeter 4810 was always full, while the
flowtube of
the Coriolis flowmeter 4820 began each batch empty. At the start of the batch,
the
totalizers in the magnetic flowmeter 4810 and in the Coriolis flowmeter 4820
were reset
and flow commenced. At the end of the batch, the shutoff valve 4850 was closed
and the
totals were frozen (hence the Coriolis flowmeter 4820 was full at the end of
the batch).
Three totals were recorded, with one each from the magnetic flowmeter 4810 and
the
Coriolis flowmeter 4820, and one from the weigh scale associated with the
weighing tank
4830. These totals are not expected to agree, as there is a finite time delay
before the
Coriolis flowmeter 4820 and then finally the weighing tank 4830 observes fluid
flow.
Thus, it would be expected that the magnetic flowmeter 4810 would record the
highest
total flow, the Coriolis flowmeter 4820 would record the second highest total,
and the
weighing tank 4830 would record the lowest total.
Fig. 48A shows the results obtained from a series of trial runs using the flow
rig
4800 of Fig 48B, with each trial run transporting about 550 kg of material
through the
flow rig. The monitored value shown is the offset observed between the weigh
scale and
either the magnetic flowmeter 4810 or the Coriolis flowmeter 4820. As
explained above,
positive offsets are expected from both instruments. The magnetic flowmeter
4810
(always full) delivers a consistently positive offset, with a repeatability
(defined here as
the maximum difference in reported value for identical experiments) of 4.0 kg.
The
analog control system associated with the magnetic flowmeter 4810 generates
large
negative offsets, with a mean of -164.2 kg and a repeatability of 87.7 kg.
This poor
performance is attributable to the analog control system's inability to deal
with the onset
of flow and the varying time taken to restart the flowtube. By contrast, the
digital
Coriolis mass flowmeter 4820 shows a positive offset averaging 25.6 kg and a
repeatability of 0.6 kg.
73
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
It would be difficult to assess the true mass-flow through the flowtube, given
its
initially empty state. The reported total mass falls between that of the
magnetic
flowmeter 4810 and the weigh scale, as expected. In an industrial application,
the issue
of repeatability is often of greater importance, as batch recipes are often
adjusted to
accommodate offsets. Of course, the repeatability of the filling process is a
lower bound
on the repeatability of the Coriolis flowmeter total. Similar repeatability
could be
achieved in an arbitrary industrial batch process. Moreover, as shown, the
digital mass
flowmeter 100 provides a substantial performance improvement over its analog
equivalent (magnetic flowmeter 4810) under the same conditions. Again, the
conclusion
drawn is that the digital mass flowmeter 100 in these conditions is not a
significant source
of measurement error.
0. Two-phase Flow
As discussed above with reference to Fig. 40A, two-phase flow, which may
result
from aeration, is another flow condition which presents difficulties for
analog control
systems and analog mass flowmeters. Two-phase flow may be sporadic or
continuous
and results when the material in the flowmeter includes a gas component and a
fluid
component traveling through the flowtube. The underlying mechanisms are very
similar
to the case of batching from empty, in that the dynamics of two-phase gas-
liquid flow
cause high damping. To maintain oscillation, a high drive gain is required.
However, the
maximum drive gain of the analog control system typically is reached at low
levels of gas
fraction in the two-phase material, and, as a result, the flowtube stalls.
The digital mass flowmeter 100 is able to maintain oscillation in the presence
of
two-phase flow. In summary, laboratory experiments conducted thus far have
been
unable to stall a tube of any size with any level of gas phase when controlled
by the
digital controller 105. By contrast, a typical analog control system stalls
with about 2%
gas phase.
Maintaining oscillation is only the first step in obtaining a satisfactory
measurement performance from the flowmeter. As briefly discussed above, a
simple
model, referred to as the "bubble" model, has been developed as one technique
to predict
the mass-flow error.
In the "bubble" or "effective mass" model, a sphere or bubble of low density
gas
is surrounded by fluid of higher density. If both are subject to acceleration
(for example
in a vibrating tube), then the bubble moves within the fluid, causing a drop
in the
74
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
observed inertia of the whole system. Defining the void fraction a as the
proportion of
gas by volume, then the effective mass drops by a proportion R, with
2
R =a
1¨a
Applied to a Coriolis flowmeter, the model predicts that the apparent mass-
flow
will be less than the true mass-flow by the factor R as will, by extension,
the observed
density. Fig. 49 shows the observed mass-flow errors for a series of runs at
different flow
rates, all using a 25 mm flowtube in horizontal alignment, and a mixture of
water and air
at ambient temperature. The x-axis shows the apparent drop in density, rather
than the
void fraction. It is possible in the laboratory to calculate the void
fraction, for example by
measuring the gas pressure and flow rate prior to mixing with the fluid,
together with the
pressure of the two-phase mixture. However, in a plant, only the observed drop
in density
is available, and the true void fraction is not available. Note that with the
analog
flowmeter, air/water mixtures with a density drop of over 5% cause flowtube
stalling so
that no data can be collected.
The dashed line 4910 shows the relationship between mass-flow error and
density
drop as predicted by the bubble model. The experimental data follow a similar
set of
curves, although the model almost always predicts a more negative mass-flow
error. As
discussed above with respect to Fig. 40A, it is possible to develop empirical
corrections to
the mass-flow rate based upon the apparent density as well as several other
internally
observed variables, such as the drive gain and the ratio of sensor voltages.
It is
reasonable to assume that the density of the pure fluid is known or can be
learned. For
example, in many applications the fluid density is relatively constant
(particularly if a
temperature coefficient is accommodated within the controller software).
Fig. 50 shows the corrected mass-flow measurements. The correction is based on
a least squares fit of several internal variables, as well as the bubble model
itself. The
correction process has only limited applicability, and it is less accurate for
lower flow
rates (the biggest errors are for 1.5-1.6 kg/second). In horizontal
orientation, the gas and
liquid phase begin to separate for lower flow rates, and much larger mass-flow
errors are
observed. In these circumstances, the assumptions of the bubble model are no
longer
valid. However, the correction is reasonable for higher flow rates. During on-
line
experimental trials, a similar correction algorithm has restricted mass flow
errors to
within about 2.5% of the mass flow reading.
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
Fig. 51 shows how a self validating digital mass flowmeter 100 responds to the
onset of two-phase flow in its reporting of the mass-flow rate. The lower
waveform 5110
shows an uncorrected mass-flow measurement under a two-phase flow condition,
and the
upper waveform 5120 shows a corrected mass-flow measurement and uncertainty
bound
under the same two-phase flow condition. With single phase flow (up to t = 7
sec.), the
mass-flow measurement is CLEAR and has a small uncertainty of about 0.2% of
the mass
flow reading. With the onset of two-phase flow, a number of processes become
active.
First, the two-phase flow is detected on the basis of the behavior of
internally observed
parameters. Second, a measurement correction process is applied, and the
measurement
status output along with the corrected measurement is set to BLURRED. Third,
the
uncertainty of the mass flow increases with the level of void fraction to a
maximum of
about 2.3% of the mass flow reading. For comparison, the uncorrected mass-flow
measurement 5110 is shown directly below the corrected mass-flow measurement
5120.
The user thus has the option of continuing operation with the reduced quality
of the
corrected mass-flow rate, switching to an alternative measurement if
available, or shutting
down the process.
P. Application of Neural Networks
Another technique for improving the accuracy of the mass flow measurement
during two-phase flow conditions is through the use of a neural network to
predict the
mass-flow error and to generate an error correction factor for correcting any
error in the
mass-flow measurement resulting from two-phase flow effects. The correction
factor is
generated using internally observed parameters as inputs to the digital signal
processor
and the neural network, and has been observed to keep errors to within 2%. The
internally observed parameters may include temperature, pressure, gain, drop
in density,
and apparent flow rate.
Fig. 52 shows a digital controller 5200 that can be substituted for digital
controller
105 or 505 of the digital mass flowmeters 100, 500 of Figs. 1 and 5. In this
implementation of digital controller 5200, process sensors 5204 connected to
the flowtube
generate process signals including one or more sensor signals, a temperature
signal, and
one or more pressure signals (as described above). The analog process signals
are
converted to digital signal data by AID converters 5206 and stored in sensor
and driver
signal data memory buffers 5208 for use by the digital controller 5200. The
drivers 5245
connected to the flowtube generate a drive current signal and may communicate
this
76
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
signal to the AID converters 5206. The drive current signal then is converted
to digital
data and stored in the sensor and driver signal data memory buffers 5208.
Alternatively,
a digital drive gain signal and a digital drive current signal may be
generated by the
amplitude control module 5235 and communicated to the sensor and driver signal
data
memory buffers 5208 for storage and use by the digital controller 5200.
The digital process sensor and driver signal data are further analyzed and
processed by a sensor and driver parameters processing module 5210 that
generates
physical parameters including frequency, phase, current, damping and amplitude
of
oscillation. A raw mass-flow measurement calculation module 5212 generates a
raw
mass-flow measurement signal using the techniques discussed above with respect
to the
flowmeter 500.
A flow condition state machine 5215 receives as input the physical parameters
from the sensor and driver parameters processing module 5210, the raw mass-
flow
measurement signal, and a density measurement 5214 that is calculated as
described
above. The flow condition state machine 5215 then detects a flow condition of
material
traveling through the digital mass flowmeter 100. In particular, the flow
condition state
machine 5215 determines whether the material is in a single-phase flow
condition or a
two-phase flow condition. The flow condition state machine 5215 also inputs
the raw
mass-flow measurement signal to a mass-flow measurement output block 5230.
When a single-phase flow condition is detected, the output block 5230
validates
the raw mass-flow measurement signal and may perform an uncertainty analysis
to
generate an uncertainty parameter associated with the validated mass-flow
measurement.
In particular, when the state machine 5215 detects that a single-phase flow
condition
exists, no correction factor is applied to the raw mass-flow measurement, and
the output
block 5230 validates the mass-flow measurement. If the controller 5200 does
not detect
errors in producing the measurement, the output block 5230 may assign to the
measurement the conventional uncertainty parameter associated with the fault
free
measurement, and may set the status flag associated with the measurement to
CLEAR. If
errors are detected by the controller 5200 in producing the measurement, the
output block
5230 may modify the uncertainty parameter to a greater uncertainty value, and
may set
the status flag to another value such as BLURRED.
When the flow condition state machine 5215 detects that a two-phase flow
condition exists, a two-phase flow error correction module 5220 receives the
raw mass-
flow measurement signal. The two-phase flow error correction module 5220
includes a
77
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
neural network processor for predicting a mass-flow error and for calculating
an error
correction factor. The neural network processor can be implemented in a
software
routine, or alternatively may be implemented as a separate programmed hardware
processor. Operation of the neural network processor is described in greater
detail below.
A neural network coefficients and training module 5225 stores a predetermined
set of neural network coefficients that are used by the neural network
processor. The
neural network coefficients and training module 5225 also may perform an
online training
function using training data so that an updated set of coefficients can be
calculated for use
by the neural network. While the predetermined set of neural network
coefficients are
generated through extensive laboratory testing and experiments based upon
known two-
phase mass-flow rates, the online training function performed by module 5225
may occur
at the initial commissioning stage of the flowmeter, or may occur each time
the
flowmeter is initialized.
The error correction factor generated by the error correction module 5220 is
input
to the mass-flow measurement output block 5230. Using the raw mass-flow
measurement and the error correction factor (if received from the error
correction module
5220 indicating two-phase flow), the mass-flow measurement output block 5230
applies
the error correction factor to the raw mass-flow measurement to generate a
corrected
mass-flow measurement. The measurement output block 5230 then validates the
corrected mass-flow measurement, and may perform an uncertainty analysis to
generate
an uncertainty parameter associated with the validated mass-flow measurement.
The
measurement output block 5230 thus generates a validated mass-flow measurement
signal
that may include an uncertainty and status associated with each validated mass-
flow
measurement, and a device status.
The sensor parameters processing module 5210 also inputs a damping parameter
and an amplitude of oscillation parameter (previously described) to an
amplitude control
module 5235. The amplitude control module 5235 further processes the damping
parameter and the amplitude of oscillation parameter and generates digital
drive signals.
The digital drive signals are converted to analog drive signals by D/A
converters 5240 for
operating the drivers 5245 connected to the flowtube of the digital flowmeter.
In an
alternate implementation, the amplitude control module 5235 may process the
damping
parameter and the amplitude of oscillation parameter and generate analog drive
signals
for operating the drivers 5245 directly.
78
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
Fig. 53 shows a procedure 5250 performed by the digital controller 5200. After
processing begins (step 5251), measurement signals generated by the process
sensors
5204 and drivers 5245 are quantified through an analog to digital conversion
process (as
described above), and the memory buffers 5208 are filled with the digital
sensor and
driver data (step 5252). For every new processing cycle, the sensor and driver
parameters
processing module 5210 retrieves the sensor and driver data from the buffers
5208 and
calculates sensor and driver variables from the sensor data (step 5254). In
particular, the
sensor and driver parameters processing module 5210 calculates sensor
voltages, sensor
frequencies, drive current, and drive gain.
The sensor and driver parameters processing module 5210 then executes a
diagmose_flow_condition processing routine (step 5256) to calculate
statistical values
including the mean, standard deviation, and slope for each of the sensor and
driver
variables. Based upon the statistics calculated for each of the sensor and
driver variables,
the flow condition state machine 5215 detects transitions between one of three
valid flow-
condition states: FLOW CONDITION SHOCK,
FLOW CONDITION HOMOGENEOUS, AND FLOW CONDITION MIXED.
If the state FLOW CONDITION SHOCK is detected (step 5258), the mass-flow
measurement analysis process is not performed due to irregular sensor inputs.
On exit
from this condition, the processing routine starts a new cycle (step 5251).
The processing
routine then searches for a new sinusoidal signal to track within the sensor
signal data and
resumes processing. As part of this tracking process, the processing routine
must find the
beginning and end of the sine wave using the zero crossing tecl-mique
described above. If
the state FLOW_CONDITION SHOCK is not detected, the processing routine
calculates
the raw mass-flow measurement of the material flowing through the flowmeter
100 (step
5260).
If two-phase flow is not detected (i.e., the state
FLOW CONDITION HOMOGENOUS is detected) (step 5270), the material flowing
through the flowmeter 100 is assumed to be a single-phase material. If so, the
validated
mass-flow rate is generated from the raw mass-flow measurement (step 5272) by
the
mass-flow measurement output block 5230. At this point, the validated mass-
flow rate
along with its uncertainty parameter and status flag can be transmitted to
another process
controller. Processing then begins a new cycle (step 5251).
If two-phase flow is detected (i.e., the state FLOW_CONDITION MIXED is
detected) (step 5270), the material flowing through the flowmeter 100 is
assumed to be a
79
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
two-phase material. In this case, the two-phase flow error correction module
5220
predicts the mass-flow error and generates an error correction factor using
the neural
network processor (step 5274). The corrected mass-flow rate is generated by
the mass-
flow measurement output block 5230 using the error correction factor (step
5276). A
validated mass-flow rate then may be generated from the corrected mass-flow
rate. At
this point, the validated mass-flow rate along with its uncertainty parameter
and status
flag can be transmitted to another process controller. Processing then begins
a new cycle
(step 5251).
Referring again to Fig. 52, the neural network processor forming part of the
two-
phase flow error correction module 5220 is a feed-forward neural network that
provides a
non-parametric framework for representing a non-linear functional mapping
between an
input and an output space. The neural network is applied to predict mass flow
errors
during two-phase flow conditions in the digital mass flowmeter. Once the error
is
predicted by the neural network, an error correction factor is applied to the
two-phase
mass flow measurement to correct the errors. Thus, the system allows the
errors to be
predicted on-line by way of the neural network using only internally
observable
parameters derived from the sensor signal, sensor variables, and sensor
statistics.
Of the various neural network models available, the multi-layer perceptron
(MLP)
and the radial basis function (RBF) networks have been used for
implementations of the
digital flowmeter. The MLP with one hidden layer (each unit having a sigmoidal
activation function) can approximate arbitrarily well any continuous mapping.
Therefore,
this type of neural network is suitable to model the non-linear relationship
between the
mass flow error of the flowmeter under two-phase flow and some of the
flowmeter's
internal parameters.
The network weights necessary for accomplishing the desired mapping are
determined during a training or optimization process. During supervised
training, the
neural network is repeatedly presented with the training set (a collection of
input
examples x, and their corresponding desired outputs d,), and the weights are
updated such
that an error function is minimized. For the interpolation problem associated
with the
present technique, a suitable error function is the sum-of-squares error,
which for an MLP
with one output may be represented as:
I=Ee1(02 E(c11 ¨y32
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
where cli is the target corresponding to input xi; yi is the actual neural
network output to xi;
and P is the number of examples in the training set.
An alternative neural network architecture that has been used is the RBF
network.
The RBF methods have their origins in techniques for performing exact
interpolations of
a set of data points in a multi-dimensional space. A RBF network generally has
a simple
architecture of two layers of weights, in which the first layer contains the
parameters of
the basis functions, and the second layer forms linear combinations of the
activation of
the basis functions to generate the outputs. This is achieved by representing
the output of
the network as a linear superposition of basis functions, one for each data
point in the
training set. In this form, training is faster than for a MLP network.
The internal sensor parameters of interest include observed density, damping,
apparent flow rate, and temperature. Each of these parameters is discussed
below.
1. Observed Density
The most widely used metric of two-phase flow is the void (or gas) fraction
defined as the proportion of gas by volume. The equation
Rz= 2a
1¨a
models the mass flow error given the void fraction. For a Coriolis mass
flowmeter, the
reported process fluid density provides an indirect measure of void fraction
assuming the
"true" liquid density in known. This reported process density is subject to
errors similar
to those in the mass-flow measurement in the presence of two-phase flow. These
errors
are highly repeatable and a drop in density is a suitable monotonic but non-
linear
indicator of void fraction, that can be monitored on-line within the
flowmeter. It should
be noted that outside of a laboratory environment, the true void fraction
cannot be
independently assessed, but rather, must be modeled as described above.
Knowledge of the "true" single-phase liquid density can be obtained on-line or
can be configured by the user (possibly including a temperature coefficient).
Both
approaches have been implemented and appear satisfactory.
For the purposes of these descriptions, the drop in density will be used as
the x-
axis parameter in graphs showing two-phase flow behavior. It should be noted
that in the
3D plots of Figs. 54 and 56-57, which capture the results from 134 on-line
experiments, a
full range of density drop points is not possible at high flow rates due to
air pressure
81
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
limitations in the flow system rig. Also, the effects of temperature, though
not shown in
the graphs, have been determined experimentally.
2. Damping
Most Coriolis flowmeters use positive feedback to maintain flowtube
oscillation.
The sensor signals provide the frequency and phase of the flowtube
oscillation, which are
multiplied by a drive gain, Ko to provide the currents supplied to the drivers
5245:
drive signal out (Amps) ID
K= ___________________________________________ =
sensor signal in (Volts) VA
Normally, the drive gain is modified to ensure a constant amplitude of
oscillation,
and is roughly proportional to the damping factor of the flowtube.
One of the most characteristic features of two-phase flow is a rapid rise in
damping. For example, a normally operating 25mm flowtube has VA= 0.3V, D=
10mA,
and hence Ko = 0.033. With two-phase flow, values might be as extreme as VA=
0.03V,
ID = 100mA, and Ko = 3.3, a hundred-fold increase. Fig. 54 shows how damping
varies
with two-phase flow.
3. Apparent Flow Rate and Temperature
As Fig. 49 shows, the mass flow error varies with the true flow rate.
Temperature
variation also has been observed. However, the true mass-flow rate is not
available
within the transmitter (or digital controller) itself when the flowmeter is
subject to two-
phase flow. However, observed (faulty) flow rate, as well as the temperature,
are
candidates as input parameters to the neural network processor.
Q. Network Training and On-Line Correction of Mass Flow Errors
Implementing the neural network analysis for the mass flow measurement error
prediction includes training the neural network processor to recognize the
mass flow error
pattern in training experimental data, testing the performance of the neural
network
processor on a new set of experimental data, and on-line implementation of the
neural
network processor for measurement error prediction and correction.
The prediction quality of the neural network processor depends on the wealth
of
the training data. To collect the neural network data, a series of two-phase
air/water
experiments were performed using an experimental flow rig 5500 shown in Fig.
55. The
flow loop includes a master meter 5510, the self validating SEVA Coriolis
flowmeter
82
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
100, and a diverter 5520 to transfer material from the flowtube to a weigh
scale 5530.
The Coriolis flowmeter 100 has a totalization function which can be triggered
by external
signals. The flow rig control is arranged so that both the flow diverter 5520
(supplying
the weigh scale) and Coriolis totalization are triggered by the master meter
5510 at the
beginning of the experiment, and again after the master meter 5510 has
observed 100 kg
of flow. The weigh scale total is used as the reference to calculate mass flow
error by
comparison with the totalized flow from the digital flowmeter 100, with the
master meter
5510 acting as an additional check. The uncertainty of the experimental rig is
estimated
at about 0.1% on typical batch sizes of 100 kg. For single-phase experiments,
the digital
flowmeter 100 delivers mass flow totals within 0.2% of the weigh scale total.
For two-
phase flow experiments, air is injected into the flow after the master meter
5510 and
before the Coriolis flowmeter 100. At low flow rates, density drops of up to
30% are
achieved. At higher flow rates, at least 15% density drops are achieved.
At the end of each batch, the Coriolis flowmeter 100 reports the batch average
of
each of the following parameters: temperature, damping, density, flow rate,
and the total
(uncorrected) flow. These parameters are thus available for use as the input
data of the
neural network processor.
The output or the target of the neural network is the mass-flow error in
percentage:
mass_en.or% = coriolis - weighscale
________________________________________________ x100
weighscale
Fig. 56 shows how the mass-flow error varies with flow rate and density drop.
Although the general trend follows the bubble model, there are additional
features of
interest. For example, for high flow rates and low density drops, the mass-
flow error
becomes slightly positive (by as much as 1%), while the bubble model only
predicts a
negative error. As is apparent from Fig. 56, in this region of the
experimental space, and
for this flowtube design, some other physical process is taking place to
overcome the
missing mass effects of two-phase flow.
The best results were obtained using only four input parameters to the neural
network: temperature, damping, density drop, and apparent flow rate. Less
satisfactory
perhaps is the result that the best fit was achieved using the neural network
on its own,
rather than as a correction to the bubble model or to a simplified curve fit.
As part of the implementation, a MLP neural network was used for on-line
implementation. Comparisons between RBF and MLP networks with the same data
sets
83
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
and inputs have shown broadly similar performance on the test set. It is thus
reasonable
to assume that the input set delivering the best RBF design also will deliver
a good (if not
the best) MLP design. The MLP neural nets were trained using a scaled
conjugate
gradient algorithm. The facilities of Neural Network Toolbox of the MATLAB
software package were used for neural network training. After further design
options
were explored, the best performance came from an 4-9-1 MLP having as inputs
the
temperature, damping, drop in density, and flow rate, along with the mass flow
error as
output.
Against a validation set, the best neural network provided mass flow rate
predictions to within 2% of the true value. Routines for the detection and
correction for
two-phase flow have been coded and incorporated into the digital Coriolis
transmitter.
Fig. 57 shows the residual mass-flow error when corrected on-line over 134 new
experiments. All errors are with 2%, and most are considerably less. The
random scatter
is due primarily to residual error in the neural network correction algorithm
(as stated
earlier, the uncertainty of the flow rig is 0.1%). Any obvious trend in the
data would
imply scope for further correction. These errors are of course for the average
corrected
mass flow rate (i.e. over a batch).
Fig. 58 shows how the on-line detection and correction for two-phase flow is
reflected in the self-validating interface generated for the mass flow
measurement. In the
graph, the lower, continuous line 5810 is the raw mass flow rate. The upper
line 5820 is a
measurement surrounded by the uncertainty band, and represents the corrected
or
validated mass-flow rate. The dashed line 5830 is the mass flow rate from the
master
meter, which is positioned prior to the air injection point (Fig. 55).
With single-phase flow (up to t = 5 sec.), the mass flow measurement has a
measurement value status of CLEAR and a small uncertainty of about 0.2% of
reading.
Once two-phase flow is detected, the neural network correction is applied
every
interleaved cycle (i.e. at 180 Hz), based on the values of the internal
parameters averaged
(using a moving window) over the last second. During two-phase flow, the
measurement
value status is set to BLURRED, and the uncertainty increases to reflect the
reduced
accuracy of the corrected measurement. The uncorrected measurement (lower dark
line)
shows a large offset error of about 30%.
The master meter reading agrees with a first approximation with the corrected
mass flow measurement. The apparent delay in its response to the incidence of
two-phase
flow is attributable to communication delays in the rig control system, and
its
84
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
squarewave-like response is due to the control system update rate of only once
per
second. Both the raw and corrected measurements from the digital transmitter
show a
higher degree of variation than with single phase flow. The master meter
reading gives a
useful measure of the water phase entering the two-phase zone, and there is a
clear
similarity between the master meter reading and the 'average' corrected
reading.
However, plug flow and air compressibility within the complex 3-D geometry of
the
flowtube may result not only in flow rate variations, but with the
instantaneous mass-flow
into the system being different from the mass-flow out of the system.
Using the self-validating sensor processing approach, the measurement is not
simply labeled good or bad by the sensor. Instead, if a fault occurs, a
correction is
applied as far as possible and the quality of the resulting measurement is
indicated
through the BLURRED status and increased uncertainty. The user thus assesses
application-specific requirements and options in order to decide whether to
continue
operation with the reduced quality of the corrected mass-flow rate, to switch
to an
alternative measurement if available, or to shut down the process. If two-
phase flow is
present only during part of a batch run (e.g. at the beginning or end), then
there will be a
proportionate weighting given to the uncertainty of the total mass of the
batch.
Multiple-Phase Flows
Fig. 59 shows an example process 5900 used to determine a phase-specific
property of a phase included in a multi-phase process fluid. For example, the
process
5900 may be used to determine the mass flow rate and density of each phase of
the multi-
phase process fluid.
As described below, an apparent intermediate value is determined based on,
e.g., a
mass flow rate and density of the multi-phase process fluid (also referred to
as the bulk
mass flow rate and bulk density, respectively) as determined by, for example,
a Coriolis
meter. Although the Coriolis meter continues to operate in the presence of the
multi-
phase process fluid, the presence of the multi-phase fluid effects the motion
of a flowtube
(or conduit) that is part of the Coriolis meter. Thus, the outputs determined
by the meter
may be inaccurate because the meter operates on the assumption that the
process fluid
includes a single phase. These outputs may be referred to as apparent
properties, or raw
properties, of the multi-phase fluid. Thus, in one implementation, the
apparent
intermediate value is determined based on apparent or raw properties of the
multi-phase
fluid. Other implementations may determine the intermediate value based on a
corrected
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
form of the apparent property(ies). To correct for the inaccuracies, the
apparent
intermediate value is input into, e.g., a neural network to produce a
corrected intermediate
value that accounts for the effects of using a multi-phase process fluid. The
corrected
intermediate value is used to determine the mass flow rate and density of each
of the
phases of the multi-phase process fluid. Using an intermediate value rather
than the
apparent mass flow rate and density of the multi-phase process fluid may help
improve
the accuracy of the determination of the mass flow rate and density of each of
the phases
of the multi-phase process fluid.
A multi-phase process fluid is passed through a vibratable flowtube (5905).
Motion is induced in the vibratable flowtube (5910). The vibratable flow tube
may be,
for example, the conduit 120 discussed above with respect to Fig. 1. The multi-
phase
process fluid also may be referred to as a multi-phase flow. The multi-phase
fluid may be
a two-phase fluid, a three-phase fluid, or a fluid that includes more than
three phases. In
general, each phase of the multi-phase fluid may be considered to be
constituents or
components of the multi-phase fluid. For example, a two-phase fluid may
include a non-
gas phase and a gas phase. The non-gas phase may be a liquid, such as oil, and
the gas
phase may be a gas, such as air. A three-phase fluid may include two non-gas
phases and
one gas phase or one non-gas phase and two gas phases. For example, the three-
phase
fluid may include a gas and two liquids such as water and oil. In another
example, the
three-phase fluid may include a gas, a liquid, and a solid (such as sand).
Additionally, the
multi-phase fluid may be a wet gas. While the wet gas may be any of the multi-
phase
fluids described above, wet gas is generally composed of more than 95% gas
phase by
volume. The process 5900 may be applied to any multi-phase fluid.
A first property of the multi-phase fluid may be determined based on the
motion
of the vibratable flowtube (5915). The first property of the multi-phase fluid
may be the
apparent mass flow rate and/or the apparent density of the fluid flowing
through the
vibratable flowtube. Thus, in the example process 5900, the first property may
be the
mass flow rate or the density of the multi-phase fluid. Properties determined
from a
multi-phase fluid may be referred to as apparent, or raw, properties as
compared to true
(or at lest corrected) properties of the multi-phase fluid. The apparent mass
flow rate and
density of the multi-phase fluid are generally inconsistent with the mass flow
rate and
density of each of the individual phases of the multi-phase fluid due to the
effects the
multi-phase fluid has on the motion of the flowtube. For example, if the multi-
phase fluid
has a relatively low gas volume fraction (e.g., the multi-phase fluid includes
more liquid
86
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
than gas), both the apparent density and the apparent mass flow rate of the
multi-phase
fluid obtained from the flowtube tend to be lower than the actual density and
mass
flowrate of the non-gas phase. Although the first property is generally an
apparent
property, in some implementations, the first property may be a corrected or
actual
property. The corrected or actual property may come, for example, from a model
or a
mapping.
As discussed with respect to Fig. 1, the mass flow rate is related to the
motion
induced in the vibratable flowtube. In particular, the mass flow rate is
related to the phase
and frequency behavior of the motion of the flowtube and the temperature of
the
flowtube. Additionally, the density of the fluid is related to the frequency
of motion and
temperature of the flowtube. Thus, because the fluid flowing through the
flowtube
includes more than one phase, the vibratable flowtube provides the mass flow
rate and
density of the multi-phase fluid rather than the mass flow rate and density of
each phase
of the multi-phase fluid. As described in more detail below, the process 5900
may be
used to determine properties of each of the phases of the multi-phase fluid.
In general, to determine the properties of the individual phases in the multi-
phase
fluid, additional information (e.g., the known densities of the materials in
the individual
phases) or additional measurements (e.g., pressure of the multi-phase fluid or
the watercut
of the multi-phase fluid) may be needed at times. However, the properties of
the multi-
phase fluid measured by the meter are typically determined by modification or
correction
to conventional single-phase measurement techniques because of the effects the
multi-
phase flow has on the flowtube as compared to single-phase flow.
Thus, in some implementations, in addition to properties determined based on
the
motion of the conduit, such as the first property discussed above, additional
or "external"
properties of the multi-phase fluid such as temperature, pressure, and
watercut may be
measured and used in the process 5900, e.g., as additional inputs to the
mapping or to
help in determining the flowrates of the individual components of the multi-
phase fluid.
The additional properties may be measured by a device other than the
flowmeter. For
example, the watercut of the multi-phase fluid, which represents the portion
of the multi-
phase fluid that is water, may be determined by a watercut meter. The
additional property
also may include a pressure associated with the flowtube. The pressure
associated with
the flowtube may be, for example, a pressure of the multi-phase process fluid
at an inlet
of the flowtube and/or a differential pressure across the flowtube.
87
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
An apparent intermediate value associated with the multi-phase process fluid
is
determined based on the first property (5920). In some implementations, a
second
property of the multi-phase fluid also may be determined based on the motion
of the
conduit. For example, in such an implementation, the apparent mass flowrate of
the
multi-phase fluid and the apparent density of the multi-phase fluid may be
determined
based on the motion of the conduit, and both of these apparent properties may
be used to
determine one or more apparent, intermediate values (such as liquid volume
fraction and
the volumetric flowrate or the gas Froude number and non-gas Froude number, as
described below). In some implementations, the apparent intermediate values
may be
intermediate values based on one or more corrected or actual properties.
In general, the apparent intermediate value (or values) is a value related to
the
multi-phase fluid that reflects inaccuracies resulting from the inclusion of
more than one
phase in the multi-phase fluid. The apparent intermediate value may be, for
example, a
volume fraction of the multi-phase process fluid. The volume fraction may be a
liquid
volume fraction that specifies the portion of the multi-phase fluid that is a
non-gas. The
volume fraction also may be a gas volume fraction that specifies the portion
of the multi-
phase fluid that is a gas. In general, the volume fraction is a dimensionless
quantity that
may be expressed as a percentage. The gas volume fraction also may be referred
to as a
void fraction. If the multi-phase fluid includes liquids and gases, the liquid
and gas
volume fractions add up to 100%. In other implementations, the apparent
intermediate
value may be a volumetric flow rate of the multi-phase fluid.
In another implementation, the apparent intermediate values may include a non-
gas Froude number and a gas Froude number. Froude numbers are dimensionless
quantities that may represent a resistance of an object moving through a fluid
and which
may be used to characterize multi-phase fluids. In this implementation, the
apparent
intermediate value may be a non-gas Froude number and/or a gas Froude number.
The
apparent gas Froude number may be calculated using the following equation,
where mg
is the apparent gas mass flow rate, pg is an estimate of the density of the
gas phase based
on the ideal gas laws, p1 is an estimate of the density of the liquid in the
non-gas phase of
the multi-phase fluid, A is the cross-sectional area of the flowtube, D is the
diameter of
the flowtube, and g is the acceleration due to gravity:
88
CA 02658810 2009-02-25
WO 2008/027908 PCT/US2007/077018
ma
Fr' = __________________________ g __ =K Va. ___
pgAVg=D\ ¨ pgPg
g
1 ma
where K= ____________________________________________ ,the Apparent Gas
Velocity V' =
Vg=D g
Similarly, the non-gas Froude number (which may be a liquid Froude number)
may be calculated using the following equation, where mia is the apparent
liquid mass
flow rate:
a
Fe mi _________ =K V' I ____ .
piA\Ig=D ¨ pg \ Pr¨ Pg
As discussed in more detail below, the apparent intermediate value is input
into a
mapping that defines a relationship between the apparent intermediate value
and a
corrected intermediate value. The mapping may be, for example, a neural
network, a
polynomial, a function, or any other type of mapping. Prior to inputting the
apparent
intermediate value into the mapping, the apparent intermediate value may be
filtered or
conditioned to reduce measurement and process noise. For example, linear
filters may be
applied to the apparent intermediate value to reduce measurement noise. The
time
constant of the linear filter may be set to a value that reflects the response
time of the
measurement instrumentation (e.g., 1 second) such that the filter remains
sensitive to
actual changes in the fluid flowing through the flowtube (such as slugs of non-
gas fluid)
while also being able to reduce measurement noise.
The development of a mapping for correcting or improving a multiphase
measurement involves the collection of data under experimental conditions,
where the
true or reference measurements are provided by additional calibrated
instrumentation.
Generally, it is not practical to carry out experiments covering all
conceivable multi-
phase conditions, either due to limitations of the test facility, and/or the
cost and time
associated with carrying out possibly thousands of experiments. Additionally,
it is rarely
possible to maintain multiphase flow conditions exactly constant for any
extended period
of time, due to the inherently unstable flow conditions that occur within
multiphase
conditions. Accordingly, it is usually necessary to calculate the average
values of all
relevant parameters, including apparent and true or reference parameter
values, over the
duration of each experiment, which may typically be of 30s to 120s duration.
Thus, the
89
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
mapping may be constructed from experimental data where each data point is
derived
from the average of for example 30s to 120s duration of data.
Difficulties might arise when applying the resulting mapping in the meter
during
multiphase flow in real time, whereby the particular parameter values observed
within the
meter are not included in the mapping provided from the previously collected
experimental data. There are two primary ways in which this may occur. In the
first
instance, although the conditions experienced by the meter, averaged over a
timescale of
about 15 to 120 seconds, do correspond to conditions covered by the mapping,
the
instantaneous parameter values may fall outside of the region, due to
measurement noise
and or instantaneous variations in actual conditions due to the instabilities
inherent in
multiphase flow. As described above, this effect can to some extent be reduced
by time-
averaging or filtering the parameters used as inputs into the mapping
function, though
there is a tradeoff between the noise reduction effects of such filtering and
the
responsiveness of the meter to actual changes in conditions within the
multiphase flow.
Alternatively, averaged parameter values may fall outside of the mapping
because, for
instance, it has not been economically viable to cover all possible multiphase
conditions
during the experimental stage.
It may not be beneficial to apply a mapping function (whether neural net,
polynomial or other function) to data that falls outside of the region for
which the
mapping was intended. Application of the mapping to such data may result in
poor
quality measurements being generated. Accordingly, jacketing procedures may be
applied to ensure that the behavior of the mapping procedure is appropriate
for parameter
values outside the mapped region, irrespective of the reasons for the
parameters falling
outside the mapped region. Data that is included in the region may be referred
to as
suitable data.
Thus, the apparent intermediate value may be "jacketed" prior to inputting the
apparent intermediate value into the mapping. For implementations that include
one
input to the mapping, the region of suitable data may be defined by one or
more limits, a
range, or a threshold. In other implementations, there may be more than one
input to the
mapping. In these implementations, the region of suitable data may be defined
by a series
of lines, curves, or surfaces. Accordingly, as the number of inputs to the
mapping
increases, defining the region of suitable data becomes more complex. Thus, it
may be
desirable to use fewer inputs to the mapping. The gas and liquid Froude
numbers
described above are an example of apparent intermediate values that may be
input into the
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
mapping without additional inputs. Thus, use of the gas and non-gas Froude
numbers
may help to reduce the number of inputs into the mapping, which also may help
reduce
the complexity of the jacketing process. Additionally, using fewer inputs to
the mapping
may result in a simpler mapping, which may help reduce the computational
resources
used by the mapping and help increase the speed of determining corrected
intermediate
values based on the mapping.
An apparent intermediate value having a value that is outside of the defined
region
may be determined to be unsuitable for input to the mapping. In general, rules
are
defined to correct an apparent intermediate value that is determined to be
outside of the
defined region. For example, an apparent intermediate value that is outside of
the defined
region may be ignored by the mapping (e.g., the apparent intermediate value is
not
corrected by the mapping), the apparent intermediate value may not be input to
the
mapping at all, a fixed correction may be applied to the apparent intermediate
value rather
than a correction determined by the mapping, or the correction corresponding
to the
correction that would apply to the value closest to the apparent intermediate
value may be
applied. Other rules for correcting an apparent intermediate value that is
outside of the
defined region may be implemented. In general, the jacketing is specific to a
particular
mapping and is defined for each mapping.
A corrected intermediate value based on a mapping between the apparent
intermediate value and the corrected intermediate value is determined (5925).
The
mapping may be a neural network, a statistical model, a polynomial, a
function, or any
other type of mapping. The neural network or other mapping may be trained with
data
obtained from a multi-phase fluid for which values of the constituent phases
are known.
Similar to the jacketing described above with respect to (5920), the corrected
apparent
value may be jacketed, or otherwise checked, prior to using it in further
processing. A
phase-specific property of the multi-phase process fluid may be determined
based on the
corrected intermediate value (5930). Using one or more of the apparent
intermediate
values discussed above rather than a value directly from the flowtube (e.g.,
mass flow rate
of the multi-phase liquid) may improve the accuracy of the process 5900. The
phase-
specific property may be, for example, a mass flow rate and/or a density of
the non-gas
and gas phases of the multi-phase fluid.
The example process described with respect to Fig. 59 may be implemented in
software or hardware. Figs. 60 and 61 describe an example implementation. With
respect to Figs. 60 and 61, optional components are shown with dashed-lines.
91
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
Specifically, Figs. 60 and 61 demonstrate the application of a digital
flowmeter to a fluid
having multiple phases that are either expected to be encountered frequently
(such as a
batch process described above) or with a fluid flow having a heterogeneous
mixture of
constituents (one or more gas constituents and/or one or more liquid
constituents).
Fig. 60 shows a digital controller 6200 that can be substituted for digital
controller
105 or 505 of the digital mass flowmeters 100, 500 of Figs. 1 and 5. In this
implementation of digital controller 6200, process sensors 6204 connected to
the flowtube
generate process signals including one or more sensor signals, a temperature
signal, and
one or more pressure signals (as described above). The analog process signals
are
converted to digital signal data by AID converters 6206 and stored in sensor
and driver
signal data memory buffers 6208 for use by the digital controller 6200. The
drivers 6245
connected to the flowtube generate a drive current signal and may communicate
this
signal to the AID converters 6206. The drive current signal then is converted
to digital
data and stored in the sensor and driver signal data memory buffers 6208.
Alternatively,
a digital drive gain signal and a digital drive current signal may be
generated by the
amplitude control module 6235 and communicated to the sensor and driver signal
data
memory buffers 6208 for storage and use by the digital controller 6200.
The digital process sensor and driver signal data are further analyzed and
processed by a sensor and driver parameters processing module 6210 that
generates
physical parameters including frequency, phase, current, damping and amplitude
of
oscillation. A raw mass-flow measurement calculation module 6212 generates a
raw
mass-flow measurement signal using the techniques discussed above with respect
to the
flowmeter 500.
Rather than include a dedicated flow condition state machine, such as 5215
discussed in reference to flowmeter 5200, a multiple-phase flow error
correction module
with one or more neural networks receives, as input, the physical parameters
from the
sensor and driver parameters processing module 6210, the raw mass-flow
measurement
signal, and a density measurement 6214 that is calculated as described above.
For
example, if the process fluid involves a known two-phase (e.g., gas and liquid
constituents), three-phase (e.g., gas and two-liquid constituents) or other
multiple-phase
flow (e.g., one or more gas and one or more liquid constituents), the
determination of a
flow condition state may not be necessary. In this example, the process fluid
may be a a
wet-gas that is already known to include a gas volume fraction (gvf) and
liquid volume
fraction (lvf). The wet gas may include, for example, natural gas, a liquid
petroleum
92
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
product, and water. Accordingly, the mass-flow measurement below can
automatically
determine a mass-flow measurement of each phase of a multi-phase process
fluid. A
dedicated neural network for each multiple phase flow condition can be used in
the
multiple phase flow error correction module 6220. Alternatively, or
additionly, a single
neural network that recognizes two-phase and/or three-phase (or more
constituent phases)
flow conditions and applies correction factors based on the actual multiple-
phase flow
condition may be used.
During a multi-phase flow condition, a multiple-phase flow error correction
module 6220 receives the raw (or apparent) mass-flow measurement signal and
the raw
density signal. The apparent mass-flow measurement and density signals reflect
the mass
flow and density of the multi-phase process fluid rather than the mass flow
and density of
each of the phases included in the multi-phase process fluid. The multiple-
phase flow
error correction module 6220 includes a neural network processor for
predicting a mass-
flow error that results from the presence of the multi-phase process fluid.
The neural
network processor can be implemented in a software routine, or alternatively
may be
implemented as a separate programmed hardware processor. Operation of the
neural
network processor is described in greater detail below.
The inputs to the neural network processor may be apparent intermediate values
determined from the raw mass-flow measurement signal and the density
measurement. In
this implementation, the multiple-phase flow error correction module 6220
determines
apparent intermediate values, such as the apparent intermediate values
discussed above
with respect to FIG. 59, from the raw (or apparent) mass flow rate and density
of the
multi-phase process fluid. The apparent intermediate values are input into the
neural
network processor and corrected. The corrected apparent intermediate values
are output
to a mass-flow measurement output block 6230. In other implementations, the
apparent
(or raw) mass-flow measurement and density may be input to the neural network.
A neural network coefficients and training module 6225 stores a predetermined
set or sets of neural network coefficients that are used by the neural network
processor for
each multiple-phase flow condition. The neural network coefficients and
training module
6225 also may perform an online training function using training data so that
an updated
set of coefficients can be calculated for use by the neural network. While the
predetermined set of neural network coefficients are generated through
extensive
laboratory testing and experiments based upon known two-phase, three-phase, or
higher-
phase mass-flow rates, the online training function performed by module 6225
may occur
93
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
at the initial commissioning stage of the flowmeter, or may occur each time
the
flowmeter is initialized.
The corrected intermediate values from the neural network are input to the
mass-flow measurement output block 6230. Using the corrected intermediate
values, the
mass-flow measurement output block 6230 determines the mass flow rate of each
phase
of the multi-phase process fluid. In some implementations, the measurement
output block
6230 validates the mass-flow measurements for the phases and may perform an
uncertainty analysis to generate an uncertainty parameter associated with the
validation.
The sensor parameters processing module 6210 also inputs a damping
0 parameter and an amplitude of oscillation parameter (previously
described) to an
amplitude control module 6235. The amplitude control module 6235 further
processes
the damping parameter and the amplitude of oscillation parameter and generates
digital
drive signals. The digital drive signals are converted to analog drive signals
by D/A
converters 6240 for operating the drivers 6245 connected to the flowtube of
the digital
5 flowmeter. In some implementations, the amplitude control module 6235 may
process
the damping parameter and the amplitude of oscillation parameter and generate
analog
drive signals for operating the drivers 6245 directly.
Fig. 61 shows a procedure 6250 performed by the digital controller 6200. After
processing begins (6251), measurement signals generated by the process sensors
6204
0 and drivers 6245 are quantified through an analog to digital conversion
process (as
described above), and the memory buffers 6208 are filled with the digital
sensor and
driver data (6252). For every new processing cycle, the sensor and driver
parameters
processing module 6210 retrieves the sensor and driver data from the buffers
6208 and
calculates sensor and driver variables from the sensor data (6254). In
particular, the
5 sensor and driver parameters processing module 6210 calculates sensor
voltages, sensor
frequencies, drive current, and drive gain.
The sensor and driver parameters processing module 6210 executes an optional
diagnose_flow_condition processing routine (6256) to calculate statistical
values
including the mean, standard deviation, and slope for each of the sensor and
driver
O variables. The optional diagnose_flow_condition processing routine (-
6256) may be
utilized, for example, to identify a two-phase flow condition and/or to
determine if a
liquid component of the two-phase flow condition includes separate liquid
constituents,
such as oil and water. Based upon the statistics calculated for each of the
sensor and
driver variables, an optional flow condition state machine (6258) may be
utilized to detect
94
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
transitions between one of three valid flow-condition states:
FLOW CONDITION SHOCK, FLOW CONDITION HOMOGENEOUS, AND
FLOW CONDITION MIXED. However, if a process fluid is known to already include
a heteregenous mixture, the process may automatically proceed from step 6254
to a
calculation of raw mass-flow measurement 6260.
If the state FLOW CONDITION SHOCK is detected (6258), the mass-flow
measurement analysis process is not performed due to irregular sensor inputs.
On exit
from this condition, the processing routine starts a new cycle (6251). The
processing
routine then searches for a new sinusoidal signal to track within the sensor
signal data and
resumes processing. As part of this tracking process, the processing routine
must find the
beginning and end of the sine wave using the zero crossing technique described
above. If
the state FLOW_CONDITION SHOCK is not detected, the processing routine
calculates
the raw mass-flow measurement of the material flowing through the flowmeter
100
(6260).
If a multiple-phase flow is already known to exist in the monitored process,
the
material flowing through the flowmeter 100 is assumed to be, for example, a
two-phase or
three-phase material. For example, the material flowing through the flowmeter
100 may
be a multi-phase process fluid, such as a wet gas. In this case, the multiple-
phase flow
error correction module 6220 determines an apparent intermediate value and,
using the
neural network processor(s), corrects the apparent intermediate value using
(6274).
Phase-specific properties of each phase of the multi-phase fluid are
determined by the
mass-flow measurement output block 6230 using the corrected intermediate value
(6276).
Processing then begins a new cycle (6251).
Referring again to Fig. 60, the neural network processor forming part of the
two-phase
flow error correction module 6220 may be a feed-forward neural network that
provides a
non-parametric framework for representing a non-linear functional mapping
between an
input and an output space. Of the various neural network models available, the
multi-
layer perceptron (MLP) and the radial basis function (RBF) networks have been
used for
implementations of the digital flowmeter. The MLP with one hidden layer (each
unit
having a sigmoidal activation function) can approximate arbitrarily well any
continuous
mapping.
In one example, a digital flowmeter 6200 may process a flow known to be a
three-
phase fluid. For example, the three-phase flow may be primarily natural gas,
with a
liquid constituent which includes a mixture of oil and water. In other
examples, the same
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
or a similar process may be applied to a two-phase fluid or a fluid containing
more than
three constituents in the fluid mixture.
Specifically, flowtube operation is maintained in a three-phase flow. The
basic
measurements of frequency and phase, sensor amplitudes and drive gain are
obtained
from the sensor signals and required current. The basic measurements are used,
with any
available external inputs and process or application-specific knowledge to
generate
estimates of the overall fluid and multi-component mass and volumetric flow
rates.
For example, the estimates of the overall fluid and multi-component mass and
volumetric flow rates may be generated as follows. The estimates of frequency,
phase,
and/or amplitudes may be improved using known correlations between the values,
such
as, for example, a rate of change of amplitude correction. The raw estimates
of the mixed
mass flow and density may be produced from best estimates of frequency, phase,
flowtube temperature and calibration constants. A simple linear correction is
applied to
the density measurement for the observed fluid pressure. In some
implementations, the
observed fluid pressure may be obtained from an external input. Because
pressure
expands and stiffens the flowtube, which may cause the raw density to be in
error, a
simple variable offset may work well if gas densities expected in a repeatable
process or
fluid mixture, while a more complete correction may include extra terms for
the variable
fluid density if changes in the liquid and/or gas constituent concentrations
are anticipated
in the process. A transmitter may include configuration parameters defining
the expected
liquid densities (with temperature compensation) and the gas reference
density.
In a three-phase fluid mixture, a fixed water-cut (WC) may be assumed or may
be
measured. The water-cut is the portion of the mixture that is water. The fluid
temperature is measured to calculate an estimate of the true fluid density
(Di) from the
water-cut and the pure oil density (D00 and water density (D water). The
estimate of the
water, =
true fluid density accounts for the known variation of Da and D water with the
fluid
temperature and the fluid pressure.
DI = WC%/100*Dwater + (1-WC%/100)*Doi1
A model (based on, for example, the ideal gas model) for the variation of gas
density (Dg) with observed fluid pressure and fluid temperature, which may be
obtained
from external inputs, is assumed and the raw liquid volume fraction (raw_LVF)
from the
raw mixture density (raw_Dm) is calculated using
96
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
raw_LVF = 100*(raw_Dm-Dg)/(D1-Dg).
The raw volumetric mixture flowrates from the raw mixture are calculated using
raw_mvf = raw_mmf / raw_Dm.
A neural network trained with experimental data is used to generate corrected
estimates of the raw liquid volume fraction and the raw volumetric flowrates
as shown
below. In the equations below, the variable "nnfunction" represents the neural
network.
corrected_LVF = nnfunction ( raw_LVF, raw_mvf, fluid_pressure, flowtube_DP )
corrected_mvf = nnfunction ( raw_LVF, raw_mvf, fluid_pressure, flowtube_DP )
The raw liquid volume fraction (raw LVF) equals 100 ¨ gas volume fraction
(GVF). Additionally, the raw liquid volume fraction is closely related to
density drop.
The raw volumetric flow can be scaled as velocity for example without change
in
approach, the Neural nets could be combined, but could use different inputs.
The liquid and gas flowrates are calculated using the following relationships:
con_liqvf = corr_LVF/100 * corr mvf
con_gasvf = (1-con_LVF/100) * con_mvf
con_liqmf = corr liqvf * D1
con_gasmf = con_gasvf * Dg
A water cut meter can be used to provide measurement used as an additional
input
to the neural network(s), and to help accurately separate the liquid flow into
the
constituent parts. To help accurately separate the liquid flow, the following
relationships
may be used:
con_Watervf = WC%/100*con_liqvf
con_Oilvf = (1-WC%/100)*corr_liqvf
97
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
corr_Watermf = con_Watervf * Dwater
con_Oilmf = con Oilvf * Doil
Alternatively or additionally, apparent gas and non-gas Froude numbers may be
determined, corrected using the neural network, and then used in determining
the mass
flowrates of the constituent components of the multiphase fluid. For example,
the gas
Froude number may be determined based on the following equation where mag is
the
apparent gas mass flow rate, pg is an estimate of the density of the gas phase
of the multi-
phase fluid based on the ideal gas laws, p1 is an estimate of the density of
the liquid in the
) non-gas phase of the multi-phase fluid, A is the cross-sectional area of
the flowtube, D is
the diameter of the flowtube, and g is the acceleration due to gravity. The
apparent gas
mass flow rate is a function of the known or assumed density of the component
fluids in
the multi-phase flow, the apparent density of the multi-phase fluid (the
apparent bulk
density), and the apparent mass flow rate of the multi-phase fluid (apparent
bulk mass
flow rate).
ma
Fr' = _______________
____________________________ = K=V: .! pg
pgAV g = D ¨ pg
1 ma
where K = __________ ,the Apparent Gas Velocity V: =
g = D pgA
Similarly, the apparent non-gas Froude number (which may be the liquid Froude
number) may be calculated using the following equation, where mia is the
apparent
liquid mass flow rate, K is the constant defined above with respect to the gas
Froude
number, and VI" is the apparent liquid velocity determined similarly to the
apparent gas
velocity shown above:
Fria = ________________________________ = K =V ____
piAV g = D ¨ Pr¨ Pg
The apparent gas and non-gas Froude numbers are then corrected using the
neural
network:
98
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
corrected_gas Froude number= nnfunction (apparent gas Froude number, apparent
non-gas Froude number)
corrected_non-gas Froude number = nnfunction (apparent gas Froude number,
apparent non-gas Froude number)
Once the corrected gas and non-gas Froude numbers are determined, the mass
flow rate for the gas and non-gas components of the multi-phase fluid may be
determined.
In particular, once the corrected values of the gas and non-gas Froude numbers
are
obtained, all parameter values for the non-gas and gas components of the multi-
phase
fluid other than the mass flow rate are known. Thus, the corrected mass flow
rate of the
non-gas and gas components of the multi-phase fluid may be determined based on
the
equations above that are used to determine the apparent Froude numbers.
Additionally, as with the implementation using the liquid volume fraction and
the
volumetric flow as inputs to the neural network, a watercut meter may be used
to help
separate the multi-phase fluid into constituent parts. For example, the
watercut meter
may provide a watercut (WC) of the multi-phase fluid that indicates the
portion of the
multi-phase fluid that is water, and the WC may be used to help separate the
multi-phase
fluid into constituent parts using the following equations:
corr_Watervf = WC%/1 00*corr_liqvf
corr_Oilvf = (1-WC%/100)*corr liqvf
con_Watermf = con_Watervf * Dwater
con_Oilmf = corr_Oilvf * Doil.
As discussed with respect to FIG. 59, in certain instances, the neural network
may
produce more accurate corrections of the apparent gas and non-gas Froude
numbers than
other apparent intermediate values. Thus, using the apparent gas and non-gas
Froude
numbers as the inputs to the neural network may result in a more accurate
determination
of the properties of the component flows that make up the multi-phase fluid.
The above description provides an overview of various digital Coriolis mass
flowmeters, describing the background, implementation, examples of its
operation, and a
comparison to prior analog controllers and transmitters. A number of
improvements over
analog controller performance have been achieved, including: high precision
control of
99
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
flowtube operation, even with operation at very low amplitudes; the
maintenance of
flowtube operation even in highly damped conditions; high precision and high
speed
measurement; compensation for dynamic changes in amplitude; compensation for
two-
phase flow; and batching to/from empty. This combination of benefits suggests
that the
digital mass flowmeter represents a significant step forward, not simply a
gradual
evolution from analog technology. The ability to deal with two-phase flow and
external
vibration means that the digital mass flowmeter 100 gives improved performance
in
conventional Coriolis applications while expanding the range of applications
to which the
flow technology can be applied. The digital platform also is a useful and
flexible vehicle
for carrying out research into Coriolis metering in that it offers high
precision, computing
power, and data rates.
An additional application of the digital flowmeter 6200 to a three-phase
fluid, e.g.,
a wet-gas having a gas (methane) and liquid component (oil and water) is
described and
shown in connection with FIGS. 62-72. Fig. 62 is a graphical view of a test
matrix for
wellheads tested based on actual testing at various well pressures and gas
velocities. Fig.
62 is a graphical view of raw density errors at various liquid void fraction
percentages
and of wells at various velocities and pressures. Fig. 64 is a graphical view
of raw mass
flow errors at various liquid void fraction percentages and of wells at
various velocities
and pressures. Fig. 65 is a graphical view of raw liquid void fraction errors
of wells at
various velocities and pressures. Fig. 66 is a graphical view of raw
volumetric flow
errors for wells at various velocities and pressures. Fig. 67 is a graphical
view of
corrected liquid void fractions of wells at various velocities and pressures.
Fig. 68 is a
graphical view of corrected mixture volumetric flow of wells at various
velocities and
pressures. Fig. 69 is a graphical view of corrected gas mass flow of wells at
various
velocities and pressures. Fig. 70 is a graphical view of corrected gas
cumulative
probability of the digital flowmeter tested. Fig. 71 is a graphical view of
corrected liquid
mass flow error of wells at various velocities and pressures. Fig. 72 is a
graphical view of
corrected gas cumulative probability of the digital flowmeter tested.
With respect to Figs. 62-72, test trials on wet gas, containing water and air,
covered test well parameters throughout a wide rang of the meter. The Liquid
Volumetric
Flowrate Fraction (LVF = 100% - GVF) points covered, included: 0.0, 0.2, 0.4,
0.6, 0.8,
1.0, 1.5, 2.0, 3.0, 4.0, 5.0%. With respect to the mass flow and density
errors detected,
errors for gas/liquid mixture mass flow and density assumed no liquid hold-up,
a static
mixer to control a stable flow profile in the field, and positive density
errors due to liquid
100
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
hold-up in meter. The positive density errors due to liquid hold-up in the
meter are at the
highest at low flows and at low gas densities. The negative mass flow errors
were similar
to Coriolis meter two-phase responses.
Applicable modeling strategies use apparent mass flow and apparent density to
apply correction factors or curve-fitting of collected data to create real
measurements,
such as a real-time density measurement. However, the wide range of gas
densities, such
as 175-900psi, favors additional approaches as well. For example, alternative
parameters
have been identified, including model parameters that include two principle
parameters
for errors. Specifically, mixture volumetric flow (assuming no slip between
phases) ¨
based on mass/density ratio and Liquid Volume Fraction (LVF) i.e. 100% - GVF.
Corrections for each are provided, using their raw values and additional
pressure data
(only). Given corrected values of LVF and volumetric flow, the mass flow rates
of the
gas and liquid components can be calculated as follows:
ml = pl . LVF/100% . Volflow
mg = pg . (1 - LVF/100%) . Volflow
The resulting errors are shown in Figs. 69-72. The model covers a wide range
of
conditions, including various pressures, and flowrates. A more restricted set
of
conditions may yield improved results, such as, for example, a higher pressure
resulting
in smaller raw errors, a "natural" operating range for the meter, very high
pressure drops
with high LVF and velocity, and/or a need to review meter sizing with wet gas.
The model may be expanded or modified so that basic pressure "corrections,"
which may include fitting raw data through curve-fitting to directly output
actual
measurements, e.g., without a true correction factor, and are applied to
density before we
apply the neural net. Current inputs are fluid specific, e.g., volumetric flow
depends on
actual fluid density. The inputs may be made less dimensional, for example, by
converting volumetric flow to velocity, then express the velocity as a
percentage of
maximum velocity that may be accommodated by the conduit, followed by
normalizing
the data to determine constituents. Operating pressures may include a 60bar
flow
pressure, with a 2-3 bar differential, and will support higher operating
pressures in the
range of 150psi-1000psi. The detailed model calibration trials referenced in
Figs. 62-72
101
CA 02658810 2009-02-25
WO 2008/027908
PCT/US2007/077018
utilize natural gas in the range of approximately 375 psi from a wellhead. The
flowtube
sizing may also be determined based on the pressure drop.
R. Source Code Listing
The following source code, which is hereby incorporated into this application,
is
used to implement the mass-flow rate processing routine in accordance with one
implementation of the flowmeter. It will be appreciated that it is possible to
implement
the mass-flow rate processing routine using different computer code without
departing
from the scope of the described techniques. Thus, neither the foregoing
description nor
the following source code listing is intended to limit the described
techniques.
Source code listing
void calculate_massflow(meas_data_type *p,
meas_data_type *op,
int validating)
double Tz, zl, z2, z3, z4, z5, z6, z7, t, x, dd, m, g,
flow error;
double noneu_mass_flow, phase_bias_unc, phase_prec_unc,
this density;
int reset, freeze;
/* calculate non- engineering units mass flow */
if (amp_svl < le-6)
noneu_mass_flow = 0.0;
else
noneu_mass_flow = tan(my_pi * p->phase_diff/180);
/* convert to engineering units */
Tz = p->temperature_value - 20;
p->massflow_value = flow_factor * 16.0 * (FC1 * Tz + FC3 * Tz*Tz
+ FC2) * noneu_mass_flow /p->v_freq;
/* apply two-phase flow correction if necessary */
if (validating && do_two_phase_correction) {
/* call neural net for calculation of mass flow correction */
= vmv_temp_stats.mean; // mean VMV temperature
= RMV_dens_stats.mean; // mean RMV density
102
CA 02658810 2013-10-02
dd = (TX_true_density - x) / TX_true_density * 100.0;
m = RMV_mass_stats.mean; // mean RMV mass flow;
g = gain_stats.mean; // mean gain;
nn_predict (t, dd, m, g, &flow_error);
p->massflow_value = 100.0*m/(100.0 + flow_error);
S. Notice of Copyright
A portion of the disclosure of this patent document contains material which is
subject to copyright protection. The copyright owner has no objection to the
facsimile
reproduction by any one of the patent document or the patent disclosure, as it
appears in the
Patent and Trademark Office patent file or records, but otherwise reserves all
copyright
rights whatsoever.
Modifications will be apparent to those skilled in the art.
103
