A two-phase flow density dynamic compensation method and system for a multi-parameter flowmeter

Abstract

The present application relates to multiphase flow metering technical field, specifically disclose a kind of two-phase flow density dynamic compensation method of multi-parameter flowmeter, upstream pressure fluctuation and downstream differential pressure fluctuation are synchronously collected on throttle element by sliding window, to distinguish transmission function with both and extract phase lag angle at crossover frequency;The difference between current phase lag angle and reference value is input as phase shift amount into erosion phase-outflow coefficient correlation model, obtain outflow coefficient decay rate and get real outflow coefficient by sliding weighted average;The real outflow coefficient is substituted into momentum conservation equation, and the cross-sectional gas rate is used as constraint, and the real mixed density is obtained by solving mixed density inverse iteration;Real mixed density is used to calculate mass flow, and the current outflow coefficient and the current phase lag angle are updated as reference value of subsequent window, the present application realizes that density compensation reference is on-line self-migration in erosion state, and compensation effectiveness can be maintained for a long time without shutdown calibration.

Description

Technical Field

[0001] This invention relates to the field of multiphase flow metering technology, specifically to a method and system for dynamic compensation of two-phase flow density in a multi-parameter flowmeter. Background Technology

[0002] In industrial processes, gas-liquid two-phase flow is widely present in pipeline transportation, wellhead metering, and process control. With the development of Industrial Internet of Things (IIoT) technology, the application of smart sensors in multi-parameter flow meters is becoming increasingly common. Current mainstream multi-parameter flow meters typically integrate smart differential pressure transmitters, temperature sensing units, and pressure sensing units; some high-end models are also equipped with smart phase content sensing elements based on sound velocity or capacitance principles. These smart sensors integrate signal conditioning, analog-to-digital conversion, and digital communication functions, enabling them to output high-frequency response raw measurement signals. This provides the necessary hardware foundation for the density compensation method based on upstream and downstream pressure pulsation signal analysis in this invention.

[0003] The existing technology has the following shortcomings: When the throttling edge of the throttling element gradually becomes passive due to long-term erosion of the sand-containing gas-liquid two-phase flow, its discharge coefficient undergoes irreversible decay. The existing density compensation method lacks real-time sensing means for the erosion process and cannot effectively distinguish between the hardware degradation component caused by erosion and the component caused by fluid density change in the differential pressure signal. At the same time, the reference value of the pre-stored correction spectrum becomes systematically invalid as the physical profile of the element changes, and this reference cannot be migrated and updated online. As a result, the density compensation model loses its reference basis and gradually fails during continuous operation, thus restricting the long-term metering stability of the multi-parameter flow meter under the condition of two-phase flow containing solid particles. Summary of the Invention

[0004] The purpose of this invention is to provide a method and system for dynamic compensation of two-phase flow density in a multi-parameter flow meter, so as to solve the problems mentioned above.

[0005] The objective of this invention can be achieved through the following technical solutions: A method for dynamic compensation of two-phase flow density in a multi-parameter flow meter includes the following steps: S1: Collect the pressure pulsation signal on the upstream side and the differential pressure pulsation signal on the downstream side of the throttling element in the multi-parameter flow meter to form two sets of time-series data within the current sliding window; S2, taking upstream pressure pulsation as excitation input and differential pressure pulsation as response output, identifies the amplitude and phase frequency characteristics of non-parametric dynamic transfer characteristics within the current window, and extracts the phase hysteresis angle at the cross frequency; S3. Compare the current phase lag angle with the factory-calibrated reference phase lag angle to obtain the phase offset. Based on the erosion phase-out coefficient correlation model, convert the phase offset into the attenuation rate of the current outflow coefficient, and then obtain the current true outflow coefficient. S4. Substitute the current true outflow coefficient into the momentum conservation equation of the throttling element, and use the synchronously measured cross-sectional gas content as a constraint to perform reverse iterative solution of the mixing density until convergence, and obtain the current true mixing density. S5 replaces the original apparent density with the current true mixing density for mass flow calculation, updates the current true outflow coefficient to the historical baseline, and sets the current phase lag angle as the new baseline phase lag angle for use in the phase offset calculation of subsequent windows.

[0006] As a further aspect of the present invention: S2 specifically includes: Construct a cross-correlation sequence using the upstream pressure pulsation sequence and the downstream differential pressure pulsation sequence within the current window; The cross-correlation sequence is fitted with a least-squares polynomial in the complex frequency domain to obtain the amplitude frequency estimation curve and phase frequency estimation curve of the current window. Locate the frequency corresponding to an amplitude ratio of 1 on the amplitude-frequency estimation curve, and extract the phase value at the corresponding frequency as the phase lag angle on the phase-frequency estimation curve.

[0007] As a further aspect of the present invention: the least squares polynomial fitting in the complex frequency domain specifically includes: Set the highest order of the fitted polynomial based on the current window length; Using the minimum weighted complex residual sum of squares at discrete frequency points of the cross-correlation sequence as the criterion, the coefficients of the corresponding order of the polynomial are fitted in progressively increasing order, and the constraint of minimizing the sum of the magnitudes of the coefficients of each order is applied during the fitting process to suppress ill-conditioned solutions. After each order of fitting is completed, the corresponding residual sum of squares is calculated. The order corresponding to the first occurrence of the minimum value of the residual sum of squares is selected as the optimal order, and the amplitude frequency estimation curve and phase frequency estimation curve under the optimal order are used as the output of the current window.

[0008] As a further aspect of the present invention: the method for constructing the erosion phase-outflow coefficient correlation model is as follows: Before leaving the factory, select standard throttling elements from the same batch and measure their reference phase hysteresis angle and reference outflow coefficient under non-erosion conditions. The grinding method gradually increases the erosion degree of the throttling element. After each erosion stage, the change in phase lag angle and the change in discharge coefficient are measured simultaneously, and the discrete corresponding points of multiple stages are recorded. Using the change in phase lag angle as the input independent variable and the change in outflow coefficient as the output dependent variable, a continuous mapping curve is constructed by piecewise conformal interpolation, and the curve is pre-stored as an erosion phase-outflow coefficient correlation model.

[0009] As a further aspect of the present invention: the calculation process of the current actual outflow coefficient is as follows: The factory-calibrated reference outflow coefficient is gradually corrected by the attenuation rate of the current window output to obtain the undetermined outflow coefficient; A sliding weighted average is applied to the sequence of undetermined outflow coefficients obtained under multiple consecutive windows, and the weighted average is calculated by taking the current window with the largest weight and the historical windows with decreasing weights in successive steps. Set the weighted average value as the actual outflow coefficient for the current window.

[0010] As a further aspect of the present invention: obtaining the current true mixing density specifically includes: The currently measured apparent density is used as the initial value for the mixing density during iteration; Substitute the current actual outflow coefficient and the initial value of the iteration into the momentum conservation equation to calculate the current differential pressure estimate, and compare it with the measured differential pressure value to obtain the differential pressure residual. Using the cross-sectional gas content as a physical boundary constraint, the iterative value of the mixing density is gradually adjusted according to the sign and amplitude of the differential pressure residual. The iteration stops when the absolute value of the difference between two adjacent differential pressure residuals is less than a preset threshold, and the current iteration value is output as the true mixing density of the current window.

[0011] As a further aspect of the present invention: the calculation of the current differential pressure estimate specifically includes: Based on the current fluid temperature and pressure state within the window, the coefficients in the momentum conservation equation are corrected online. Substitute the current actual outflow coefficient and the initial value of the iteration into the momentum conservation equation after state correction to calculate the initial differential pressure estimate. The weighted cumulative amount of differential pressure residuals during historical iterations is used as a feedback correction term to dynamically compensate the initial differential pressure estimate and output the current differential pressure estimate.

[0012] As a further aspect of the present invention: S5 specifically includes: The actual mixed density output in the current window is used as the density value in the current mass flow calculation; After completing the mass flow calculation for the current window, the current actual outflow coefficient is stored as the historical baseline outflow coefficient for the next window; The current phase lag angle is stored as the reference phase lag angle for the next window, overwriting the reference value of the previous window, so that the phase offset calculation of subsequent windows is based on the latest reference.

[0013] A two-phase flow density dynamic compensation system for a multi-parameter flow meter includes: The signal acquisition module acquires the pressure pulsation signal upstream of the throttling element and the differential pressure pulsation signal downstream of the multi-parameter flow meter, forming two sets of time-series data within the current sliding window; The transmission feature extraction module uses upstream pressure pulsation as excitation input and differential pressure pulsation as response output to identify the amplitude and phase frequency features of non-parametric dynamic transmission characteristics within the current window, and extracts the phase hysteresis angle at the cross frequency. The erosion correction and coefficient calculation module compares the current phase lag angle with the factory-calibrated reference phase lag angle to obtain the phase offset. Based on the erosion phase-outflow coefficient correlation model, the phase offset is converted into the attenuation rate of the current outflow coefficient, and then the current true outflow coefficient is obtained. The density inverse solution module substitutes the current true outflow coefficient into the momentum conservation equation of the throttling element, and uses the synchronously measured cross-sectional gas content as a constraint to perform inverse iterative solution of the mixing density until convergence, thus obtaining the current true mixing density. The benchmark update and metering application module replaces the original apparent density with the current true mixed density for mass flow calculation, updates the current true outflow coefficient to the historical benchmark, and sets the current phase lag angle as the new benchmark phase lag angle for use in the phase offset calculation of subsequent windows.

[0014] The beneficial effects of this invention are: (1) This invention extracts the phase lag angle at the cross frequency in the dynamic transmission characteristics of the throttling element as the erosion characterization quantity. It utilizes the physical property that the phase shift caused by erosion passivation is independent of the phase content and flow velocity changes of the fluid to separate the hardware degradation caused by erosion from the two-phase flow density change in the signal feature domain. This avoids the density compensation failure caused by the inability to distinguish the source of differential pressure amplitude change in traditional methods. It provides an online self-sensing compensation basis that does not rely on external calibration for harsh two-phase flow conditions containing solid particles.

[0015] (2) This invention constructs an erosion phase-outflow coefficient correlation model and uses a sliding weighted average to update the real outflow coefficient. At the same time, the real mixed density obtained in each window is used for mass flow calculation. The current outflow coefficient and the current phase lag angle are then rolled to update the reference value of the subsequent window, forming a closed-loop update link where the reference value evolves over time. This allows the density compensation reference of the flow meter to automatically migrate with the gradual degradation of the geometric profile of the throttling element, and maintains the compensation effectiveness for a long time without stopping the machine for disassembly and calibration. Attached Figure Description

[0016] The invention will now be further described with reference to the accompanying drawings.

[0017] Figure 1 This is a flowchart of the method of the present invention; Figure 2 This is a flowchart illustrating the construction process of the erosion phase-outflow coefficient correlation model in this invention. Figure 3 This is a system block diagram of the present invention. Detailed Implementation

[0018] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0019] Please see Figure 1 As shown, this invention provides a method for dynamic compensation of two-phase flow density in a multi-parameter flow meter, comprising the following steps: S1: Collect the pressure pulsation signal on the upstream side and the differential pressure pulsation signal on the downstream side of the throttling element in the multi-parameter flow meter to form two sets of time-series data within the current sliding window; S2, taking upstream pressure pulsation as excitation input and differential pressure pulsation as response output, identifies the amplitude and phase frequency characteristics of non-parametric dynamic transfer characteristics within the current window, and extracts the phase hysteresis angle at the cross frequency; S3. Compare the current phase lag angle with the factory-calibrated reference phase lag angle to obtain the phase offset. Based on the erosion phase-out coefficient correlation model, convert the phase offset into the attenuation rate of the current outflow coefficient, and then obtain the current true outflow coefficient. S4. Substitute the current true outflow coefficient into the momentum conservation equation of the throttling element, and use the synchronously measured cross-sectional gas content as a constraint to perform reverse iterative solution of the mixing density until convergence, and obtain the current true mixing density. S5 replaces the original apparent density with the current true mixing density for mass flow calculation, updates the current true outflow coefficient to the historical baseline, and sets the current phase lag angle as the new baseline phase lag angle for use in the phase offset calculation of subsequent windows.

[0020] In S1, the pressure pulsation signal upstream of the throttling element and the differential pressure pulsation signal downstream of the multi-parameter flowmeter are collected to form two sets of time-series data within the current sliding window, specifically including: During the data acquisition process, a first pressure sensing element is installed upstream of the throttling element of the multi-parameter flow meter. The first pressure sensing element is a piezoresistive pressure sensor, whose measuring end extends into the fluid being measured to sense the instantaneous pressure fluctuation of the upstream fluid in real time and output a first voltage signal corresponding to the instantaneous pressure. A second pressure sensing element is installed downstream of the throttling element of the multi-parameter flow meter. This second pressure sensing element is also a piezoresistive pressure sensor. Its measuring end extends into the fluid on the downstream side to sense the instantaneous pressure fluctuation of the fluid on the downstream side in real time and output a second voltage signal corresponding to the instantaneous pressure. The first voltage signal and the second voltage signal are respectively connected to the high-pressure input terminal and the low-pressure input terminal of the same differential pressure transmitter. The differential pressure transmitter contains a signal conditioning circuit and an analog-to-digital conversion circuit, which are used to convert the real-time pressure difference between the high-pressure end and the low-pressure end into a digital differential pressure instantaneous value sequence. At the same time, the first voltage signal is output as an upstream instantaneous pressure value sequence after analog-to-digital conversion.

[0021] The instantaneous pressure value sequence and differential pressure instantaneous value sequence on the upstream side are synchronously acquired at a fixed sampling frequency. This sampling frequency is not less than five times the natural frequency of the throttling element to ensure that the two sequences are strictly aligned on the time axis. Each sampling moment corresponds to an instantaneous upstream pressure value and an instantaneous differential pressure value. After collecting a preset number of sampling points consecutively, the latest preset number of sampling points constitute the current sliding window. The value of the preset length is not less than the minimum number of sampling points required for the transition process of the throttling element to the step excitation response. Once the current sliding window is formed, the instantaneous upstream pressure values ​​corresponding to all sampling times within the window are arranged in chronological order to form the upstream pressure pulsation time series data for the current window. The instantaneous differential pressure values ​​corresponding to the same time within the window are also arranged in chronological order to form the differential pressure pulsation time series data for the current window. Whenever a new sampling point is acquired, the sliding window advances by one sampling interval, discards the earliest sampling point data in the window and incorporates the latest acquired sampling point data, forming two sets of time-series data for the next sliding window. This enables the sliding window to continuously scroll and update on the time axis for feature identification in subsequent windows.

[0022] In S2, using upstream pressure pulsation as the excitation input and differential pressure pulsation as the response output, the amplitude and phase frequency characteristics of the non-parametric dynamic transfer characteristics are identified within the current window, and the phase lag angle at the cross-frequency is extracted, specifically including: Within the current window, retrieve the upstream pressure pulsation time series data and the differential pressure pulsation time series data. The two sets of time series data have the same sequence length, which is recorded as the first length.

[0023] For each delay value within the range from negative first length minus 1 to positive first length minus 1, multiply each value in the upstream pressure pulsation time series data from the first value to the (first length minus the absolute value of the current delay) value with the corresponding value in the differential pressure pulsation time series data from (the absolute value of the current delay plus 1) to the first length value digit by digit, and then sum them up. Divide the sum by (first length minus the absolute value of the current delay) to obtain the cross-correlation value corresponding to the current delay value. After traversing all delay values, construct the cross-correlation sequence of the current window with each delay value as the x-axis and the cross-correlation value corresponding to each delay value as the y-axis.

[0024] Transform the cross-correlation sequence of the current window from the time domain to the frequency domain to obtain the complex values ​​of the cross-correlation sequence at each discrete frequency point; set the highest order of the fitting polynomial according to the first length of the current window, where the highest order is the first length divided by 30 and rounded down, and the highest order does not exceed 20. Using the complex values ​​at each discrete frequency point as the fitting target, the polynomial coefficients corresponding to each order from the 1st to the highest order are fitted in an order-by-order incremental manner. The fitting criterion is to minimize the weighted complex residual square sum of the complex values ​​at each discrete frequency point. The weight of each discrete frequency point is set in an inversely proportional manner from low to high frequency. In the fitting process of each order, the constraint condition of minimizing the sum of the magnitudes of the coefficients of each order is applied to the polynomial coefficients to suppress ill-conditioned solutions caused by uneven distribution of discrete frequency points.

[0025] After completing the polynomial fitting for each order, the difference between the complex values ​​at each discrete frequency point and the complex values ​​of the cross-correlation sequence at the corresponding discrete frequency point is calculated. The squared magnitude of each difference is multiplied by the weight of its corresponding discrete frequency point and then summed to obtain the residual sum of squares for that order. Starting from order 1, the residual sums of squares for each order are compared sequentially. When the residual sum of squares for a certain order is less than the residual sum of squares for the previous order and simultaneously less than the residual sum of squares for the next order, that order is determined to be the optimal order. If no such minimum value appears in the entire order range, the order corresponding to the smallest residual sum of squares is determined to be the optimal order.

[0026] The coefficients of the polynomial at the optimal order are used to construct the numerator and denominator polynomials of the transfer function for the current window. With frequency as the independent variable, the ratio of the amplitude of the numerator polynomial at each frequency point to the amplitude of the denominator polynomial at the same frequency point is calculated to obtain the amplitude-frequency estimate at each frequency point, and the amplitude-frequency estimate curve is plotted. The difference between the phase value of the numerator polynomial at each frequency point and the phase value of the denominator polynomial at the same frequency point is calculated to obtain the phase-frequency estimate at each frequency point, and the phase-frequency estimate curve is plotted.

[0027] On the amplitude-frequency estimation curve, the frequency is scanned point by point from low frequency to high frequency. When the amplitude-frequency estimation value at a certain frequency point is equal to the value 1, the frequency point is determined as the crossover frequency. If the amplitude-frequency estimation value is not exactly equal to 1, the frequency point closest to 1 before the amplitude-frequency estimation value is less than 1 and its two adjacent frequency points are selected. Local parabolic interpolation is performed using the amplitude-frequency estimation values ​​and frequency values ​​of these three points. The frequency corresponding to the amplitude value equal to 1 after interpolation is determined as the crossover frequency. The phase value corresponding to the crossover frequency is taken from the phase-frequency estimation curve. This phase value is output as the phase lag angle of the current window for use in subsequent steps.

[0028] Please see Figure 2As shown, in S3, the current phase lag angle is compared with the factory-calibrated reference phase lag angle to obtain the phase offset. Based on the erosion phase-outflow coefficient correlation model, the phase offset is converted into the attenuation rate of the current outflow coefficient, and then the current true outflow coefficient is calculated. Specifically, this includes: During the calibration stage before leaving the factory, the throttling element processed in the same batch as the flow meter to be shipped is selected as the standard part. The standard throttling element is installed in the calibration pipeline, and the upstream pressure pulsation and downstream differential pressure pulsation are measured under pure water single-phase flow conditions. The reference phase lag angle of the standard throttling element under non-erosion state is obtained according to the S1 to S2 method. At the same time, the reference outflow coefficient under this state is calibrated by weighing method.

[0029] Diamond polishing paste was used to grind the throttling edge of the standard throttling element step by step. After each grinding stage, a profilometer was used to measure the radius of the arc of the throttling edge and record the grinding amount to simulate the geometric passivation caused by sand erosion in the field. After each grinding stage, the standard throttling element was reinstalled in the calibration pipeline, and the phase lag angle and discharge coefficient were measured again under the same flow conditions. The difference between the phase lag angle after each grinding stage and the reference phase lag angle was recorded as the phase offset, and the attenuation ratio of the discharge coefficient after each grinding stage relative to the reference discharge coefficient was recorded as the attenuation rate. Multiple sets of discrete corresponding points were recorded step by step.

[0030] Using the phase offset as input and the attenuation rate as output, the above-mentioned multiple sets of discrete corresponding points are arranged in ascending order of phase offset. Between two adjacent discrete points, a cubic conformal interpolation method is used to construct an interpolation polynomial. This interpolation polynomial ensures that the first derivative at the discrete point is continuous and its monotonicity remains unchanged, thus forming a continuous mapping curve from the phase offset to the attenuation rate. This continuous mapping curve is stored in the flowmeter's storage unit in the form of discrete nodes and interpolation polynomial coefficients, serving as an erosion phase-outflow coefficient correlation model.

[0031] During online operation, after obtaining the current phase lag angle by pressing S2 within a certain window, the factory-calibrated reference phase lag angle is retrieved from the storage unit. The difference between the current phase lag angle and the reference phase lag angle is calculated, and the absolute value of this difference is taken as the current phase offset. The current phase offset is input into the aforementioned erosion phase-outflow coefficient correlation model. The interpolation interval where the current phase offset is located is found on the continuous mapping curve. The corresponding attenuation rate is calculated by substituting it into the cubic conformal interpolation polynomial of that interval. The calculation method is as follows: Let the current phase offset be... The two endpoints of the interpolation interval in which it lies are and The corresponding attenuation rates are respectively and ,and Then the attenuation rate The calculation formula is: ; in, The attenuation rate corresponding to the current phase offset. This is the attenuation rate calibration value corresponding to the left endpoint of the interpolation interval. This is the attenuation rate calibration value corresponding to the right endpoint of the interpolation interval. This is the current phase offset. This is the calibration value for the phase offset corresponding to the left endpoint of the interpolation interval. The current phase offset is the calibration value corresponding to the right endpoint of the interpolation interval. When the current phase offset is less than the minimum phase offset in the calibration data, the attenuation rate corresponding to the minimum phase offset is taken as the current attenuation rate. When the current phase offset is greater than the maximum phase offset in the calibration data, the attenuation rate corresponding to the maximum phase offset is taken as the current attenuation rate.

[0032] The reference outflow coefficient read from the storage unit is multiplied by the difference between the value 1 and the aforementioned attenuation rate to obtain the unsmoothed outflow coefficient of the current window. The validity of the unsmoothed outflow coefficient is checked. If the unsmoothed outflow coefficient is less than 0.7 times or greater than 1.0 times the reference outflow coefficient, it is determined to be an invalid value. The unsmoothed outflow coefficient of the current window is discarded and the real outflow coefficient of the previous window is used as the real outflow coefficient of the current window.

[0033] For multiple consecutive valid windows, the undetermined outflow coefficients corresponding to each window are cached to form a sequence of undetermined outflow coefficients. A sliding weighted average is applied to this sequence, where the weight coefficient for the current window is 0.5, the weight coefficient for the previous window is 0.3, the weight coefficient for the second-to-last window is 0.15, and the weight coefficients for the third-to-last and earlier windows are set in an exponentially decreasing manner with an exponential decay factor of 0.5. The undetermined outflow coefficients of each window are multiplied by their corresponding weight coefficients and summed. The sum is then divided by the sum of all weight coefficients to obtain the weighted average. During the calculation, if there are fewer than four valid windows, only the existing valid windows are weighted, and the weight coefficients of each window are scaled proportionally to make the total weight equal to 1. This weighted average is then set as the actual outflow coefficient output for the current window for use in subsequent steps.

[0034] In S4, the current true outflow coefficient is substituted into the momentum conservation equation of the throttling element, and the mixing density is solved in reverse iteration until convergence, constrained by the synchronously measured cross-sectional gas content, to obtain the current true mixing density. Specifically, this includes: Within the current window, the apparent density value of the two-phase flow given by the densitometer or empirical correlation at the previous moment is extracted and used as the initial value for the mixing density iteration. This initial value is only used to start the reverse iteration process of this window, and its accuracy does not affect the final convergence result of the iteration.

[0035] Obtain the cross-sectional gas content measurement value synchronously measured in the current window, and set the cross-sectional gas content measurement value as the physical boundary constraint for the reverse iteration of the mixing density. The meaning of this constraint is: the mixing density value after each iteration adjustment must meet the mixing density value range determined by the cross-sectional gas content, pure gas phase density, and pure liquid phase density.

[0036] Within the current window, the current actual outflow coefficient and the trial value of the mixing density used in the current iteration are substituted into the momentum conservation equation of the throttling element. This equation describes the quantitative relationship between mass flow rate, outflow coefficient, throttling ratio, expandable coefficient, throttling area, mixing density and differential pressure. After substituting, a differential pressure calculation value is obtained. The differential pressure calculation value is subtracted from the measured differential pressure value in the current window to obtain the differential pressure residual of the current iteration step.

[0037] In the first iteration step, the apparent density is used as the initial value for iteration and substituted into the momentum conservation equation to calculate the initial differential pressure estimate, and subsequent corrections are made to obtain the current differential pressure estimate. Before calculating the initial differential pressure estimate, the coefficients in the momentum conservation equation are corrected online based on the measured fluid temperature and fluid pressure values ​​within the current window. The correction methods are as follows: temperature correction is performed on the fluid density and fluid viscosity based on the current temperature value, and the correction coefficient is obtained by consulting the fluid property database; pressure correction is performed on the expandability coefficient based on the current pressure value, and the correction coefficient is obtained by consulting the fluid property database.

[0038] Substitute the current actual outflow coefficient and the initial iteration value into the momentum conservation equation after state correction, and multiply them in the order of the current values ​​of the parameters on the right side of the equation to calculate the initial differential pressure estimate.

[0039] After the initial differential pressure estimate is calculated, the weighted cumulative amount of differential pressure residuals from each iteration step in the historical iteration process is used as a feedback correction term to dynamically compensate the initial differential pressure estimate. The historical iteration steps are all the iteration steps executed within the current window. The weighting method is as follows: the historical differential pressure residuals closer to the current iteration step are assigned a larger weight coefficient, and the historical differential pressure residuals farther from the current iteration step are assigned a smaller weight coefficient. The weight coefficients decrease according to the exponential decay law, and the exponential decay factor is set to 0.6. All historical differential pressure residuals are multiplied by their corresponding weight coefficients and then summed to obtain the weighted cumulative amount. The initial differential pressure estimate is subtracted from this weighted cumulative amount to obtain the current differential pressure estimate for the current iteration step.

[0040] Subtract the current differential pressure estimate from the measured differential pressure value in the current window to obtain the differential pressure residual for the current iteration step; determine whether the absolute value of the difference between the differential pressure residuals of two adjacent iteration steps is less than a preset convergence threshold, which is one-thousandth of the measured differential pressure value in the current window; if it is less than the preset convergence threshold, stop the iteration and output the trial value of the mixing density used in the current iteration step as the true mixing density of the current window.

[0041] If the absolute value of the difference between the differential pressure residuals of two adjacent iterations is greater than or equal to the preset convergence threshold, the iterative value of the mixing density is adjusted according to the sign and amplitude of the current differential pressure residual. The adjustment method is as follows: if the differential pressure residual is positive, it indicates that the current trial value of the mixing density is too large, so the mixing density is adjusted in the decreasing direction; if the differential pressure residual is negative, it indicates that the current trial value of the mixing density is too small, so the mixing density is adjusted in the increasing direction. The adjustment step size is the ratio of the absolute value of the differential pressure residual to the current differential pressure estimate multiplied by one-tenth of the current trial value of the mixing density. The adjusted mixing density value also needs to be boundary-checked with the above-mentioned cross-sectional gas content constraint. If the adjusted mixing density value exceeds the physical value range determined by the cross-sectional gas content, the boundary value of the value range is taken as the new iterative value of the mixing density, and then the next iteration step is entered. The above differential pressure estimate calculation and residual comparison process is repeated until the convergence condition is met.

[0042] In S5, the current true mixing density replaces the original apparent density for mass flow calculation, and the current true outflow coefficient is updated to the historical baseline. Simultaneously, the current phase lag angle is set as the new baseline phase lag angle for subsequent window phase offset calculations. Specifically, this includes: After obtaining the true mixed density through S4 in the current window, the true mixed density is directly assigned to the mass flow calculation channel of the flow meter, replacing the original apparent density value in the channel that was directly measured by the density meter or estimated by empirical correlation, and used as the density value in the mass flow calculation of the current window.

[0043] The multi-parameter flow meter synchronously obtains the volumetric flow rate measurement value of the current window through the velocity sensing element, multiplies the volumetric flow rate measurement value with the above-mentioned actual mixing density, and obtains the mass flow rate value of the current window; the mass flow rate value is output as the final measurement result of the current window to the flow meter's display unit or data bus for the host computer or field control system to read and use.

[0044] After the mass flow rate calculation for the current window is completed, the actual outflow coefficient of the current window obtained in S3 is read from the working register and written to the historical reference storage area specially allocated in the storage unit to overwrite the original historical reference outflow coefficient in the storage area. After the write operation is completed, the outflow coefficient stored in the historical reference storage area is updated to the actual outflow coefficient of the current window, which can be retrieved and used by S3 of the next window when calculating the phase offset and attenuation rate.

[0045] Simultaneously, the current window phase lag angle extracted in S2 is read from the working register and written to the specially allocated reference phase storage area in the storage unit to overwrite the original factory-calibrated reference phase lag angle or the previously updated reference phase lag angle in the storage area; after the writing operation is completed, the phase lag angle stored in the reference phase storage area is updated to the phase lag angle of the current window.

[0046] When the next sliding window starts, the phase offset calculation in S3 will directly retrieve the latest stored reference phase lag angle as the comparison benchmark, instead of retrieving the factory-calibrated reference phase lag angle. This ensures that the phase offset calculation of subsequent windows is based on the phase lag angle of the latest window. At the same time, the latest historical reference outflow coefficient retrieved when calculating the current true outflow coefficient in S3 is also updated to the latest written current window true outflow coefficient.

[0047] For a continuously rolling sliding window sequence, starting from the second window, the S3 of each window is calculated from the baseline phase lag angle updated in the previous window and the historical baseline outflow coefficient, thus forming a closed-loop link in which the baseline values ​​of the phase lag angle and the outflow coefficient are gradually updated over time.

[0048] In the first window after the flow meter is first powered on or reset, since there is no updated reference value from the previous window, the reference phase lag angle is taken from the factory calibration value, and the historical reference outflow coefficient is taken from the initial factory calibration outflow coefficient; from the second window onwards, all the above reference values ​​are switched to the updated values ​​from the previous window.

[0049] When the S4 reverse iteration fails to converge within a certain window, causing the window to be unable to output an effective true mixing density, the S5 baseline value update operation is not performed within that window. The baseline phase lag angle and historical baseline outflow coefficient in the storage unit remain unchanged to maintain the continuity of the baseline value and avoid introducing erroneous baselines due to invalid windows.

[0050] Both the historical reference storage area and the reference phase storage area in the storage unit adopt a single-value overwrite writing method, that is, only the latest written value is retained each time it is updated, and the entire reference value of the storage window is not accumulated, so as to control the storage space occupation and ensure the reading efficiency. Each write operation is executed only after all calculations in the current window have been completed and the iteration convergence has been confirmed to be valid, ensuring that the reference values ​​written to the storage unit are all valid results of the valid window.

[0051] Please see Figure 3 As shown, a two-phase flow density dynamic compensation system for a multi-parameter flow meter includes: The signal acquisition module acquires the pressure pulsation signal upstream of the throttling element and the differential pressure pulsation signal downstream of the multi-parameter flow meter, forming two sets of time-series data within the current sliding window; The transmission feature extraction module uses upstream pressure pulsation as excitation input and differential pressure pulsation as response output to identify the amplitude and phase frequency features of non-parametric dynamic transmission characteristics within the current window, and extracts the phase hysteresis angle at the cross frequency. The erosion correction and coefficient calculation module compares the current phase lag angle with the factory-calibrated reference phase lag angle to obtain the phase offset. Based on the erosion phase-outflow coefficient correlation model, the phase offset is converted into the attenuation rate of the current outflow coefficient, and then the current true outflow coefficient is obtained. The density inverse solution module substitutes the current true outflow coefficient into the momentum conservation equation of the throttling element, and uses the synchronously measured cross-sectional gas content as a constraint to perform inverse iterative solution of the mixing density until convergence, thus obtaining the current true mixing density. The benchmark update and metering application module replaces the original apparent density with the current true mixed density for mass flow calculation, updates the current true outflow coefficient to the historical benchmark, and sets the current phase lag angle as the new benchmark phase lag angle for use in the phase offset calculation of subsequent windows.

[0052] The working principle of this invention is as follows: Upstream pressure pulsations and downstream differential pressure pulsations of the throttling element are simultaneously acquired via a sliding window. The upstream pressure pulsations serve as the excitation, and the differential pressure pulsations as the response. Within the window, the amplitude and phase frequency characteristics of the transfer function are identified, and the phase lag angle at the cross-frequency is extracted. The difference between the current phase lag angle and the reference phase lag angle is used as the phase offset input to a pre-constructed erosion phase-outflow coefficient correlation model. The outflow coefficient attenuation rate is calculated, and after correcting the factory calibration reference value, a sliding weighted average is used to obtain the current true outflow coefficient. This true outflow coefficient is substituted into the momentum conservation equation, and the synchronously measured cross-sectional gas content is used as the physical boundary constraint. The mixing density is then iteratively solved in reverse until convergence to obtain the current true mixing density. This true mixing density replaces the original apparent density for mass flow rate calculation. Simultaneously, the current true outflow coefficient and the current phase lag angle are updated to the historical reference and the new reference phase lag angle, respectively, for use in subsequent sliding window phase offset calculations.

[0053] The foregoing has provided a detailed description of one embodiment of the present invention, but this description is merely a preferred embodiment and should not be construed as limiting the scope of the invention. All equivalent variations and modifications made within the scope of the claims of this invention should still fall within the patent coverage of this invention.

Claims

1. A method for dynamic compensation of two-phase flow density in a multi-parameter flowmeter, characterized in that, Includes the following steps: S1: Collect the pressure pulsation signal on the upstream side and the differential pressure pulsation signal on the downstream side of the throttling element in the multi-parameter flow meter to form two sets of time-series data within the current sliding window; S2, taking upstream pressure pulsation as excitation input and differential pressure pulsation as response output, identifies the amplitude and phase frequency characteristics of non-parametric dynamic transfer characteristics within the current window, and extracts the phase hysteresis angle at the cross frequency; S3. Compare the current phase lag angle with the factory-calibrated reference phase lag angle to obtain the phase offset. Based on the erosion phase-out coefficient correlation model, convert the phase offset into the attenuation rate of the current outflow coefficient, and then obtain the current true outflow coefficient. S4. Substitute the current true outflow coefficient into the momentum conservation equation of the throttling element, and use the synchronously measured cross-sectional gas content as a constraint to perform reverse iterative solution of the mixing density until convergence, and obtain the current true mixing density. S5 replaces the original apparent density with the current true mixing density for mass flow calculation, updates the current true outflow coefficient to the historical baseline, and sets the current phase lag angle as the new baseline phase lag angle for use in the phase offset calculation of subsequent windows.

2. The two-phase flow density dynamic compensation method for a multi-parameter flowmeter according to claim 1, characterized in that, S2 specifically includes: Construct a cross-correlation sequence using the upstream pressure pulsation sequence and the downstream differential pressure pulsation sequence within the current window; The cross-correlation sequence is fitted with a least-squares polynomial in the complex frequency domain to obtain the amplitude frequency estimation curve and phase frequency estimation curve of the current window; Locate the frequency corresponding to an amplitude ratio of 1 on the amplitude-frequency estimation curve, and extract the phase value at the corresponding frequency as the phase lag angle on the phase-frequency estimation curve.

3. The two-phase flow density dynamic compensation method for a multi-parameter flowmeter according to claim 2, characterized in that, The least squares polynomial fitting in the complex frequency domain specifically includes: Set the highest order of the fitted polynomial based on the current window length; Using the minimum weighted complex residual sum of squares at discrete frequency points of the cross-correlation sequence as the criterion, the coefficients of the corresponding order of the polynomial are fitted in progressively increasing order, and the constraint of minimizing the sum of the magnitudes of the coefficients of each order is applied during the fitting process to suppress ill-conditioned solutions. After each order of fitting is completed, the corresponding residual sum of squares is calculated. The order corresponding to the first occurrence of the minimum value of the residual sum of squares is selected as the optimal order, and the amplitude frequency estimation curve and phase frequency estimation curve under the optimal order are used as the output of the current window.

4. The two-phase flow density dynamic compensation method for a multi-parameter flowmeter according to claim 1, characterized in that, The method for constructing the erosion phase-outflow coefficient correlation model is as follows: Before leaving the factory, select standard throttling elements from the same batch and measure their reference phase lag angle and reference outflow coefficient under non-erosion conditions. The grinding method gradually increases the erosion degree of the throttling element. After each erosion stage, the change in phase lag angle and the change in discharge coefficient are measured simultaneously, and the corresponding points of multiple discrete stages are recorded. Using the change in phase lag angle as the input independent variable and the change in outflow coefficient as the output dependent variable, a continuous mapping curve is constructed by piecewise conformal interpolation, and it is pre-stored as an erosion phase-outflow coefficient correlation model.

5. The two-phase flow density dynamic compensation method for a multi-parameter flowmeter according to claim 1, characterized in that, The calculation process for the current actual outflow coefficient is as follows: The factory-calibrated reference outflow coefficient is gradually corrected by the attenuation rate of the current window output to obtain the undetermined outflow coefficient; A sliding weighted average is applied to the sequence of undetermined outflow coefficients obtained under multiple consecutive windows, and the weighted average is calculated by taking the current window with the largest weight and the historical windows with decreasing weights in successive steps. Set the weighted average value as the actual outflow coefficient for the current window.

6. The two-phase flow density dynamic compensation method for a multi-parameter flowmeter according to claim 1, characterized in that, Obtaining the current true mixing density specifically includes: The currently measured apparent density is used as the initial value for the mixing density during iteration; Substitute the current actual outflow coefficient and the initial value of the iteration into the momentum conservation equation to calculate the current differential pressure estimate, and compare it with the measured differential pressure value to obtain the differential pressure residual. Using the cross-sectional gas content as a physical boundary constraint, the iterative value of the mixing density is gradually adjusted according to the sign and amplitude of the differential pressure residual. The iteration stops when the absolute value of the difference between two adjacent differential pressure residuals is less than a preset threshold, and the current iteration value is output as the true mixing density of the current window.

7. The two-phase flow density dynamic compensation method for a multi-parameter flowmeter according to claim 6, characterized in that, The calculation of the current differential pressure estimate specifically includes: Based on the current fluid temperature and pressure state within the window, the coefficients in the momentum conservation equation are corrected online. Substitute the current actual outflow coefficient and the initial value of the iteration into the momentum conservation equation after state correction to calculate the initial differential pressure estimate. The weighted cumulative amount of differential pressure residuals during historical iterations is used as a feedback correction term to dynamically compensate the initial differential pressure estimate and output the current differential pressure estimate.

8. The two-phase flow density dynamic compensation method for a multi-parameter flowmeter according to claim 1, characterized in that, S5 specifically includes: The actual mixed density output in the current window is used as the density value in the current mass flow calculation; After completing the mass flow calculation for the current window, the current actual outflow coefficient is stored as the historical baseline outflow coefficient for the next window; The current phase lag angle is stored as the reference phase lag angle for the next window, overwriting the reference value of the previous window, so that the phase offset calculation of subsequent windows is based on the latest reference.

9. A two-phase flow density dynamic compensation system for a multi-parameter flowmeter, characterized in that, A method for dynamic compensation of two-phase flow density in a multi-parameter flowmeter according to any one of claims 1-9 includes: The signal acquisition module acquires the pressure pulsation signal upstream of the throttling element and the differential pressure pulsation signal downstream of the multi-parameter flow meter, forming two sets of time-series data within the current sliding window; The transmission feature extraction module uses upstream pressure pulsation as excitation input and differential pressure pulsation as response output to identify the amplitude and phase frequency features of non-parametric dynamic transmission characteristics within the current window, and extracts the phase hysteresis angle at the cross frequency. The erosion correction and coefficient calculation module compares the current phase lag angle with the factory-calibrated reference phase lag angle to obtain the phase offset. Based on the erosion phase-outflow coefficient correlation model, the phase offset is converted into the attenuation rate of the current outflow coefficient, and then the current true outflow coefficient is obtained. The density inverse solution module substitutes the current true outflow coefficient into the momentum conservation equation of the throttling element, and uses the synchronously measured cross-sectional gas content as a constraint to perform inverse iterative solution of the mixing density until convergence, thus obtaining the current true mixing density. The benchmark update and metering application module replaces the original apparent density with the current true mixed density for mass flow calculation, updates the current true outflow coefficient to the historical benchmark, and sets the current phase lag angle as the new benchmark phase lag angle for use in the phase offset calculation of subsequent windows.

Images