Abstract:
The occurrence and magnitude of a leak in a pipeline are determined by sensing a product velocity, pressure and temperature at an upstream point of flow of a product of known composition in said pipeline and sensing product velocity, pressure and temperature at a downstream point in said pipeline. 
     The ambient temperature of the material surrounding said pipeline is sensed at spaced locations along said pipeline between said upstream and downstream points. 
     The above parameters are employed to determine from four equations (conservation of mass, conservation of momentum, conservation of energy, and the equation of state) a single valued function of time representative of the product imbalance of product flowing from said upstream point to said downstream point.

Description:
TECHNICAL FIELD 
     This invention relates to operation of fluid carrying pipelines and more particularly to the rapid detection of the occurrance or existence and magnitude of a leak therein. 
     BACKGROUND ART 
     In the operation of pipelines through which fluids such as crude oil, refined products, natural gas and the like flow under pressure of pumping systems, it is highly desirable to be able to portray on real-time basis any loss of product that might occur through a rupture in the pipeline or through some other cause. 
     Heretofore there have been many attempts to detect such events in order to direct efforts to locate and repair pipeline leaks. 
     Heretofore attempts have been made to determine loss, utilizing pressure and flow velocity measured at two points along a pipeline. Such attempts are described in the following publications: 
     &#34;Network System Transient Calculations by Implicit Method&#34;, M. A. Stoner, Journal of the Society of Petroleum Engineers, December 1971, pp. 352-365; and 
     &#34;Some Applications of Transient Flow Simulation&#34;, V. L. Streeter and E. B. Wylie, Journal of the Society of Petroleum Engineers, December 1978, pp. 357-364. 
     DISCLOSURE OF THE INVENTION 
     In accordance with the present invention, variable pipeline parameters are measured at two spaced points along the pipeline. Upstream, the variables are: (1) product velocity, (2) product pressure, (3) product temperature, and (4) product composition. At the downstream point, the variables are (1) product velocity, (2) pressure and (3) temperature. Further measurements are made at various points along the pipeline and between the upstream and downstream points of the temperature of the earth at the level at which the pipeline is buried, or in the case of surface pipelines, the temperature is determined on a continuing basis at spaced points. 
     Based upon the above parameters, a function is generated represented by the solution of four simultaneous equations. 
     The equations are: 
     1. conservation of mass equation, 
     2. conservation of momentum equation, 
     3. conservation of energy equation, and 
     4. the equation of state. 
     Any product imbalance is registered on a real-time basis and as a function of time. The plotted or registered function is established through optimizing the solutions in accordance with the finite difference explicit method which is well known. Integration steps are made as long as possible in terms of the nature of three variables: (1) product velocity, (2) pipeline size, and (3) the pipeline wall thickness. 
     More particularly, in accordance with the present invention, the occurrence and magnitude of a pipeline leak is determined by: 
     (a) sensing product velocity, pressure, and temperature of flow of product of known composition in said pipeline at an upstream point, 
     (b) sensing product velocity, pressure, and temperature at a downstream point in said pipeline, 
     (c) determining ground temperature at spaced locations between said upstream and downstream points, 
     (d) utilizing the measured parameters for the simultaneous solution of: (1) the conservation of mass equation, (2) the conservation of momentum equation, (3) conservation of energy equation, and (4) the equation of state for said pipeline with said product therein, and 
     (e) registering the product imbalance function resulting from repetitive solutions to said equations as function of time. The imbalance is equal to the product fed into the pipeline minus the product delivered from the pipeline minus the product packed into the pipeline. The imbalance appearing on a real-time chart will provide indicia as to the time occurrance and magnitude of a leak or rupture. 
    
    
     BRIEF DESCRIPTION OF DRAWINGS 
     FIG. 1 illustrates a pipeline system embodying the present invention. 
    
    
     DETAILED DESCRIPTION 
     Referring now to FIG. 1, a pipeline 10 is buried beneath the earth&#39;s surface 11 and extends from station A to station B. Stations A and B may be spaced apart many miles. Ordinarily they would be spaced in accordance with the locations of pumping stations. As shown in FIG. 1 a measuring unit 12 is located at station A and a measuring unit 13 is located at station B. 
     The present invention is directed to the determination of the occurrence and magnitude of a leak such as might occur at any point along line 10, such as at point 14. The product may be gas or liquids. It is desirable, as soon as possible after occurrence, that the existence of a leak is determined so that its location can then be established and corrective action taken. 
     In accordance with the present invention, a unique set of parameters of the line 10 and of the product in line 10 are utilized repetitively in order to provide a substantially continuous time function as to the magnitude of any product imbalance between the quantity leaving the station 12 and quantity arriving at pump station 13. This is accomplished in accordance with the present invention by utilizing the following measurements. 
     At station A an analyzer 20 is provided to indicate on the output channel 21 the density of the product in line 10. A sensor 22 is provided to indicate on channel 23 the velocity of the product in line 10 at station A. The sensor 24 provides an output on channel 25 indicative of the pressure in the product at station A and sensor 26 provides an input on channel 27 representative of the temperature of the product. A further sensor 28 provides an input on channel 29 indicating the temperature of the ground in the region of the pipeline 10. Input channels 21, 23, 25, 27 and 29 are connected to a processor 30. 
     Along line 10, ground temperature sensors 31 and 32 are located and are connected to processor 30. Sensors 31 and 32 are located at intermediate points along line 10. 
     At station B, a ground temperature sensor 33 is located. Also at station B is a sensor 34 which provides on channel 35 a signal representing the velocity of the product at station B. Sensor 36 provides on channel 37 a signal representing the pressure in the product at station B. Sensor 38 provides on channel 39 a signal representative of the temperature of the product at station B. 
     With the above variables applied to processor 30, four simultaneous equations are repeatedly solved to yield a single function solution which is plotted as a time dependent trace 41 on chart recorder 40. The trace 41 is the imbalance between the quantity of product at station A and the quantity of product delivered at station B. 
     If the nature of the product is known, the use of analyzer 20 is not necessary and a signal on channel 21 may be provided, on a suitable compatible scale, proportional to the density of the product. 
     Many types of transient flow conditions are possible in a system depending mainly upon its geometry and upon the cause of the transient. The unsteady flow can be of very extended duration as in starting up a long transmission system; its duration may be in the order of one day as in the supply to industrial, commercial or residential areas, or it may be short term, perhaps one hour in the case of unexpected mechanical failures or power shutdowns. Even faster transients become significant upon operation of relief valve mechanisms, or more importantly breaks or leaks on line 10. 
     CONTROLLING RELATIONSHIPS 
     Several physical equations used for steady and unsteady state analysis are set forth below in differential equation form: 
     
         ______________________________________(i)  Conservation of Mass Equation ##STR1##(ii) Conservation of Momentum Equation ##STR2##(iii)Conservation of Energy Equation ##STR3## K(T - T.sub.g)πD = ##STR4##(iv) The Equation of State ##STR5## ##STR6## ##STR7##______________________________________Where:     For Propane*     For N-Butane*______________________________________A.sub.0    =          18,634.7    =      32,544.7B.sub.0    =          .964762     =      1.56588C.sub.0    =          7.96178 × 10.sup.9                      =      1.37436 × 10.sup.10D.sub.0    =          4.53708 × 10.sup.11                      =      3.33159 × 10.sup.11E.sub.0    =          2.56053 × 10.sup.13                      =      2.30902 × 10.sup.12a   =          4.00664 × 10.sup.4                      =      7.11818 × 10.sup.4b   =          5.46248     =      9.14006c   =          2.74461 × 10.sup.10                      =      7.00044 × 10.sup.10d   =          1.50520 × 10.sup.7                      =      3.64238 × 10.sup.7α    =          2.01402     =      4.00985γ    =          4.56182     =      7.54122ρ    =  dρ/dtρ    = densityv   = velocityE   = wall elasticityμ    = poisson&#39;s ratio = 1.33ξ.sub.1    = lateral stress tensor = 5/4ξ.sub.2    = hoop stress tensor = 7/3D   = densityg   = gravityh   = elevation changet   = timeA   = areaT   = temperatureT.sub.g    = ground temperatureR   = gas constantM   = molecular weightx   = distance along the pipeline______________________________________ *From Properties of Light Petroleum Fluids, Kenneth D. Starling, Gulf Publishing Co., 1972. 
    
     Within the equation of momentum, the frictional factor and elevation head terms are produced. Within the mass equation, pipewall deformation is accounted for. In the given form, the equations operate on four (4) variables, T, P, v, and ρ. 
     The mass and momentum equation contain no property of the product other than density. Thus the operation can be employed for any product, given that the density ρ can be produced for the product from an adequate equation of state. 
     An accurate equation of state ranging through the critical regions is employed so that the Joule-Thompson effect arises from a direct calculation of the energy flow. Heat conduction to and from the ground is included. Boundary conditions are pressure and temperature upstream and mass flow downstream, all of which can be varied linearly over selected time intervals. In actual operations, they are provided by the real-time data from instruments. 
     Frictional factor in the above equations is computed on an interval-by-interval basis, using the known Colebrook-White relation discussed in the Proceedings of the Royal Society in London, Serial A, Vol. 161, p. 367 (1973). Viscosity is computed using correlation such as the Dean and Stiel correlation discussed in &#34;Viscosity of Non-Polar Gas Mixtures at Moderate and High Pressure&#34;, Journal of American Institute of Chemical Engineers, Vol. 11, p. 525, 1965. 
     More than one product in line 10, with a diffusive interface between them may be accommodated. The interface is treated as a sharp transition from one product to the other. An auxiliary routine determines the growth and shape of the interface due to turbulent diffusion. The sharp interface approximation is accurate if the length of the interface is small compared to the selected integration interval along the pipeline. 
     Useful operations are accurate, fast, flexible, and take every significant parameter into account. Such desirable properties are not necessarily compatible with each other or with a processor that may be available. The present operation reconciles these requirements as far as possible. Ultimately, decisions must be made as to what is most needed. Several methods, such as characteristics, exist for the numerical integration of simultaneous differential equations. The operation herein preferably uses a first order finite difference approach, with careful attention to the protocol of incrementing the various dependent variables. Finite difference integration is fast, simple and retains all the essential elements in a discernible form. It also tends to be very stable. 
     Accuracy and speed depend upon the size and distance steps. The nature of flow in pipelines is such that surprisingly long (miles) distance integration increments for integration give accurate results. This is important, because the integration time increment must be significantly less than the time required for an acoustic wave to travel down a distance increment or else instabilities become unmanageable. As a result, the time step employed is much less than would be needed for accuracy alone. 
     A fluid flows down a pipeline with its acceleration determined by the balance between pressure gradient and friction. Friction is usually expressed as some frictional factor multiplied by a function of velocity. The frictional factor itself is a slowly varying function of the fluid properties in the pipe. The work done on the fluid by the pressure gradient force goes into heat due to the friction. Since this work comes from enthalpy of the fluid, the flow would be isothermal, except for conduction to the ground and the Joule-Thompson effect. A pipeline may be thought of as a long porous plug. Because many products are pumped at near critical regions, an adequate equation of state, both to relate state variables and permit determination of thermodynamic functions at high speeds, is most difficult. 
     Leaks have been heretofore detected externally by observing some effect of the escaping product, or internally by observation of the product loss. Although chemical detection is extraordinarily sensitive, the cost of continual sampling along a pipeline has been limited generally to the observation of dead vegetation. Loss of product can be observed by monitoring mass flow at two stations. However, without an accurate transient dependent operation to represent the fluctuating line pack, hours are required to detect 10% leaks and days are required to detect 1% leaks. The present transient dependent operation in real-time, with boundary conditions supplied and updated from pipeline instrumentation and supervisory control systems, permits the identification of leaks within the resolution of the instruments (typically 1%) in the time required for the effects of the leak to propagate to the instruments. 
     Transient responsive operations are herein made available with greater flexibility so that most pipeline operations can be simulated with speed to permit operational use in real-time. The product properties which are important in lighter products, particularly in the transition through the critical region, is handled by an equation of state covering a wide range. The operation, once tuned to the pipeline with adequate instrumentation, permits rapid identification of leaks and improper conditions. In one embodiment of the computer program for use in the system of FIG. 1 was as set out in Table I. 
     
                                           TABLE I__________________________________________________________________________C  DYNAME, FR.C  ** THIS SUBROUTINE PERFORMS THE DYNAMIC GAS MODEL ONC  ** A SINGLE PIPELINE SEGMENTC  ** &lt;THIS VERSION DOES NOT SUPPORT INTERFACES&gt;C  **C  **C  ** THE VARIABLES UPON ENTRY:C  ** PRODUCT (PRODUCT ID) - THE ID OF THE PRODUCT ON THE HEAD OFC  **         THE SEGMENTC  ** FRFACT  (FRICTION FACTOR) - THE FRICTION FACTOR OF THE PRODUCTC  **         IN THE HEAD OF THE SEGMENTC  **C--   SUBROUTINE DYNAM (SEGNUM)   COMMON/PIPE1/PRS(40),TEM(40),DEN(40),VEL(40)   COMMON/PIPE2/MIL(40),GTM(40),HLC(40),GRV(40)   COMMON/CONFG/HEAD(10),TAIL(10),SUCT(10),TLOC(10)   COMMON/SEGDX/LEN(10)   COMMON/CONS1/FOURODIAM,TWODI   COMMON/DUNKP/ITNUM,PRODUCT,FRFACT,KOKMI,MOLWTI,KMIMI   COMMON/PARMS/DT,NUMOFSEGS   DIMENSION DELTAV(40),DELTAD(40),DELTAT(40)   REAL LEN,MIL,MOLWT,MOLWTI,KOKMI,KMIMI   INTEGER HEAD,TAIL,SUCT,TLOC   INTEGER START,STOP,SEGNUM,PRODUCTC  **   PARAMETER R = 49706.5  ; = 32.17 * 10.73 * 144.0   PARAMETER G = 32.17C  **   DELTAXI = 1.0/LEN(SEGNUM)                     ; GET SEGMENT&#39;S INTERVAL LEN.C  **C  ** SET THE ABSOLUTE INTEGRATION LIMITS FOR THIS SEGMENTC  **10 START = HEAD(SEGNUM)   ; GET THE FIRST INTERVAL   STOP = TAIL(SEGNUM)    ; GET THE LAST INTERVAL NUMBERC  **C  ** NOW ESTABLISH THE VELOCITY PROFILE USING A FORWARD STEP OFC  ** THE CONVERSATION OF MOMENTUM EQUATIONC  **J = 0                     ; DELTAV ARRAY INDEX   DO 100 I = START,STOP - 1   J = J + 1   DELTAV(J) = ((VEL(I)*(VEL(I) - VEL(I + 1)) + (PRS(1) - PRS(I +   1))/DEN(I) &amp; + (GRV(I) - GRV(I + 1)))*DELTAXI - FRFACT*VEL(I)*ABS(VEL(I))*TWOD1)   *DT   VEL(I) = VEL(I) + DELTAV(J)100   CONTINUEDELTAV(J + 1) = 0.0       ; A BOUNDARY CONDITIONC  **C  ** NOW USE THE GAS ENERGY EQUATIONS TO DETERMINE THE INTERVALC  ** TEMPERATURE PROFILES.C  **J = 1                     ; DELTA ARRAYS INDEXC  **C  ** COMPUTE THE (I - 1) TERM TO START THINGS OFFC  **   I = START   H = KOKMI*R*TEM(1)*MOLWTI   TERM2 = DEN(I)*VEL*(0.5*VEL(I)*VEL(1) + H &amp; + GRV(I) + G*ENTHD(DEN(I),TEM(I),PRODUCT)*MOLWT1)C  **C  ** THIS IS THE MAIN INTEGRATION LOOPC  **   DO 250 I = START + 1,STOPJ = J + 1                 ; DELTA ARRAYS INDEXC  **C  ** COMPUTE THE NEW DENSITY FOR THIS INTERVALC  **   DELTAD(J) = (DEN(I - 1)*VEL(I - 1) - DEN(I)*VEL(I))*DELTAXI*DTC  **C  ** COMPUTE THE VARIOUS TERMS OF THE ENERGY EQUATION.C  ** NOTE THE THE (I - 1) TERM HAS BEEN SAVED FROM THE PREVIOUS STEP.C  **TERM1 = TERM2             ; FROM THE PREVIOUS STEP   H = KOKMI*R*TEM(I)*MOLWTI   TERM2 = DEN(I)*VEL(I)*(0.5*VEL(I)*VEL(I) + H &amp; + GRV(I) + G*ENTHD(DEN(I),TEM(I),PRODUCT)*MOWLTI)   TERM3 = HLC(I)*(GTM(I) - TEM(I))*FOURODIAM   U = H - PRS(I)/DEN(I)   TERM4 = 0.5*DEN(I)*VEL(I)*DELTAV(J) + (0.5*VEL(I)*VEL(I) &amp; + U + G*ENGD(DEN(I),TEM(I),PRODUCT)*MOLWTI)*DELTAD(J)   TERM5 = (U/TEM(I) + G*SPHD(DEN(1),TEM(I),PRODUCT)*MOLWTI)*DEN(I)   DELTAT(J) = (((TERM1 - TERM2)*DELTAXI + TERM3)*DT - TERM4)/TERM5C  **250   CONTINUEC  **C  ** NOW IT IS TIME TO ADJUST THE DENSITIES AND TEMPERATURES TO REFLECTC  ** THE INCREMENTS JUST COMPUTED. (NOTE THAT PRESSURE IS COMPUTED FROMC  ** THE EQUATION OF STATE)C  **   J = 1   DO 450 I = START + 1,STOP   J = J + 1   DEN(I) = DEN(I) + DELTAD(J)   TEM(I) = TEM(I) + DELTAT(J)   PRS(I) = STARL(DEN(I),TEM(I))450   CONTINUEC  **   RETURNC  **   END__________________________________________________________________________